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

Назначение

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

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

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

Пусть известны также  M значений  r_k = r ( t_k ) импульсной функции отклика  r ( t ) на той же равномерной сетке  t_k = k Δt,  k = - (M - 1)/2, ..., 0, ..., (M - 1)/2 (где  Δt - тот же шаг сетки, что и для функции сигнала, а  M - нечетное натуральное число такое, что M < N). Значения  r_k задаются в  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) выполняется восстановление значений  s_j функции сигнала по известным значениям  r_k импульсной функции ответа, задаваемым в массиве RESP, и значениям соответствующей свертки, задаваемым в массиве ARRAY_. Полученные значения  s_j помещаются в первых  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);

   Параметры 

   Версии 

   Вызываемые подпрограммы 

   Замечания по использованию  

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

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)