Текст подпрограммы и версий
rsa1r_p.zip, rsa1e_p.zip 		Тексты тестовых примеров
trsa1r_p.zip, trsa1e_p.zip

Подпрограмма:  RSA1R (модуль RSA1R_p)

Назначение

Вычисление прямой и обратной дискретной периодической (круговой) свертки двух одномерных массивов вещественных чисел.

Математическое описание

Пусть известны  N значений  sj = s ( tj ) функции сигнала  s ( t ) на равномерной сетке  tj = j Δt,  j = 0, 1, 2, ..., N - 1 (где  Δt - шаг сетки, а  N равняется целой степени двух). Предполагается, что сигнал  sj является периодическим с периодом  N (отсюда, в частности, следует, что значения  sj известны и для отрицательных значений  j, т.е.  s- j = sN - j,  j = 1, 2, ..., N - 1). Значения  sj задаются в вещественном одномерном массиве ARRAY_ длины  N:  ARRAY_ ( J ) = sj ( J = 1, 2, ..., N;  j = 0, 1, 2, ..., N - 1).

Пусть известны также  M значений  rk = r ( tk ) импульсной функции отклика  r ( t ) на той же равномерной сетке  tk = k Δt,  k = - (M - 1)/2, ..., 0, ..., (M - 1)/2 (где  Δt - тот же шаг сетки, что и для функции сигнала, а  M - нечетное натуральное число такое, что M < N). Значения  rk задаются в  M первых элементах одномерного вещественного массива RESP длины  N следующим образом: 

       RESP(1)  =  r0
      RESP( I )  =  rk ,     I = 2, 3, ..., (M+1)/2 ;   k = 1, 2, ..., (M-1)/2
      RESP( I )  =  rk ,     I = (M+3)/2, ..., M ;      k = (M-1)/2, ..., - 1

Подпрограмма RSA1R имеет два режима работы, задаваемых при обращении к ней значением параметра IREG. В первом режиме ( IREG = 1) выполняется прямая дискретная периодическая свертка массивов ARRAY_ и RESP (L = 1, 2, ..., N): 

                               N
             (r*s)L   =   ∑   ARRAY_(L - J) × RESP( J )
                              j =1

где  * - знак свертки, а  × - знак умножения.
Предварительно в массив RESP добавляются нули следующим образом: 

   RESP( I )  =  RESP( I )                 ,   I = 1, 2, ..., (M+1)/2 
   RESP( I )  =  0.0                           ,   I = (M+3)/2, ..., N - (M-1)/2 
   RESP(N+1-I )  =  RESP(M+1-I ) ,   I = 1, 2, ..., (M-1)/2

Для вычисления указанной свертки используется быстрое дискретное преобразование Фурье в соответствии с теоремой о дискретной периодической (круговой) свертке. Полученные значения (r * s)L помещаются в первых  N компонентах массива ANS.

Во втором режиме ( IREG = - 1) выполняется восстановление значений  sj функции сигнала по известным значениям  rk импульсной функции ответа, задаваемым в массиве RESP, и значениям соответствующей свертки, задаваемым в массиве ARRAY_. Полученные значения  sj помещаются в первых  N элементах массива ANS.

Р.Отнес, Л.Эноксон. Прикладной анализ временных рядов. Изд - во "Мир", 1982.

Использование 

procedure RSA1R(var ARRAY_ :Array of Real; var N :Integer;
                var RESP :Array of Real; M :Integer;
                var ANS :Array of Complex; var FFT :Array of Complex;
                IREG :Integer; var IFLAG :Integer);


   Параметры
ARRAY_ - вещественный одномерный массив длины  N, содержащий либо значения  sj ( IREG = 1), либо значения свертки ( IREG = - 1);
N - длина массива ARRAY, равная целой степени двух (тип: целый);
RESP - вещественный одномерный массив длины  N, в первых  M компонентах которого задаются значения  rk;
M - количество заданных значений  rk (тип: целый);
ANS - вещественный одномерный массив длины 2N, в первых  N элементах которого в результате работы подпрограммы размещаются либо компоненты вектора свертки (IREG = 1), либо восстановленные значения  sj (IREG = - 1);
FFT - комплексный одномерный массив длины  N, используемый в подпрограмме в качестве рабочего;
IREG - задает режим работы подпрограммы (тип: целый); при этом:
IREG=  1 - когда вычисляется свертка;
IREG= -1 - когда восстанавливается функция сигнала  sj
IFLAG - целая переменная, служащая для сообщения о характере окончания работы подпрограммы во втором режиме ( IREG = - 1); при этом:
IFLAG=0 - когда восстановлены все значения  sj;
IFLAG=I - когда  I - ая компонента преобразования Фурье импульсной функции отклика  rk равна нулю; это означает, что значение  sI не может быть восстановлено из - за потери информации в векторе свертки; значения sj в этом случае подпрограммой не восстанавливаются. 
   Версии
RSA1E - вычисление прямой и обратной дискретной периодической (круговой) свертки двух одномерных массивов вещественных чисел в режиме расширенной (Extended) точности; при этом параметры ARRAY_, RESP и ANS должны иметь тип Extended, параметр FFT - тип Complex * 16 (это означает, что данная версия существует только для тех компиляторов, входной язык которых допускает комплексные данные расширенной (Extended) точности). 
   Вызываемые подпрограммы
       FTA2R -
       FTA2E   		одновременное выполнение прямого быстрого дискретного преобразования Фурье двух одномерных массивов вещественных чисел, каждый из которых имеет длину  N, равной целой степени двух, в режимах одинарной и расширенной (Extended) точности; используются в подпрограммах RSA1R и RSA1E соответственно.
       FTA3R -
       FTA3E   		выполнение прямого и обратного быстрых преобразований Фурье одномерного массива вещественных чисел длины 2N, где  N равняется целой степени двух, в режимах одинарной и расширенной (Extended) точности; используются в подпрограммах RSA1R и RSA1E соответственно. 
   Замечания по использованию
  В подпрограммах RSA1R и RSA1E проверка того, что значение  N должно быть целой степенью двух, а  M должно быть нечетным числом, меньшим  N, не производится.

Пример использования 

Unit TRSA1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, RSA1R_p;

function TRSA1R: String;

implementation

function TRSA1R: String;
var
N,M,_i,IFLAG :Integer;
RESP :Array [0..15] of Real;
ANS :Array [0..31] of Real;
FFT :Array [0..15] of Complex;
const
ARRАY_ :Array [0..15] of Real = ( 0.0,0.0,0.5,0.0,0.0,-0.5,0.0,0.0,0.5,0.5,0.5,
0.5,0.5,0.0,0.0,0.0 );
begin
Result := '';  { результат функции }
N := 16;
M := 3;
RESP[0] := 0.7;
RESP[1] := 0.4;
RESP[2] := 0.4;
RSA1R(ARRAY_,N,RESP,M,ANS,FFT,1,IFLAG);
Result := Result + Format('%5d ',[IFLAG]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 31 do
 begin
  Result := Result + Format('%20.16f ',[ANS[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TRSA1R',Result);  { вывод результатов в файл TRSA1R.res }
exit;
end;

end.

Результаты:

       ANS = (- 0.596046E - 7, 0.2, 0.35, 0.2, - 0.2, - 0.35, - 0.2, 
                       0.2, 0.55, 0.75, 0.75, 0.75, 0.55, 0.2, - 0.745058E - 8,  
                       0.521541E - 07)