Текст подпрограммы и версий
rss1r_p.zip 		Тексты тестовых примеров
trss1r_p.zip

Подпрограмма:  RSS1R (модуль RSS1R_p)

Назначение

Построение автопериодограмм и взаимных периодограмм двух вещественных случайных процессов по одному отрезку их реализаций.

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

Пусть заданы реализации двух вещественных, стационарных, центрированных случайных процессов  Xq,  Yq,  q = 0, 1, ..., M - 1 на сетке по оси времени  tq = q * DT с шагом DT > 0. Рассмотрим отрезки этих реализаций длины  N:  1 ≤ N ≤ M, левые границы которых  a, b:  0 ≤ a ≤ M - N,  0 ≤ b ≤ M - N, вообще говоря, различны, т.е. отрезки вида:  Xa + s,  Yb + s,  s = 0, 1, ..., N - 1, и вычислим их дискретные преобразования Фурье  Am, Bm: 

                  N-1
      Am  =    ∑   Xa + s exp(- i λm ts )  =
                  s=0
                                                              N-1
                                                          =   ∑   Xa + s exp(- 2π i m s / N) ,
                                                              s=0
                  N-1
      Bm  =    ∑   Yb + s exp(- i λm ts )  =
                  s=0
                                                              N-1
                                                          =   ∑   Yb + s exp(- 2π i m s / N) ,
                                                              s=0

где   λm = m * Δλ,   m = 0, 1, ..., N - 1,   Δλ = 2π/(N * DT),   i = √-1.

Соответствующие автопериодограммы и взаимные периодограммы строятся по формулам [1]: 

    PXm    =  (DT / N) | Am | 2 ,            PYm  =  (DT / N)  | Bm | 2 ,

    PXYm  =  (DT / N)  A*m Bm ,             m  =  0, 1, ..., N -1 ,

где * - знак комплексного сопряжения, при этом функции РХm, PYm и PXY1m = Re PXYm симметричны относительно точки  m = N/2, а функция PXY2m = Im PXYm  антисимметрична относительно  m = N/2.

Периодограммы, полученные по одному отрезку реализаций, являются при N → +∞ асимптотически несмещенными оценками спектров случайных процессов, но эти оценки асимптотически несостоятельны.

Полное описание реализованного алгоритма содержится в статье [2] (подпрограмма PERSEG).

1.  Дж.Бендат, А.Пирсол, Измерение и анализ случайных процессов, Изд - во "Мир", M., 1974.
2.  М.В.Арефьева, Корреляционный и спектральный анализ стационарных случайных процессов (часть 1), сб. "Численный анализ на ФОРТРАНе", вып.15. Изд - во МГУ, M., 1976.

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

procedure RSS1R(var X1 :Array of Real; var Y1 :Array of Real;
                var N :Integer; DT :Real; var PX :Array of Real;
                var PY :Array of Real; var PXY1 :Array of Real;
                var PXY2 :Array of Real);


   Параметры
X1, Y1 - одномерные мссивы длины N + 1, первые  N элементов которых содержат значения исходных реализаций случайных процессов Xa + s,  Yb + s, s = 0, 1, ..., N - 1 (тип: вещественный);
N - количество заданных значений каждого из двух исходных рядов, N ≥ 4 - целая степень числа два (тип: целый);
DT - заданный шаг сетки по оси времени, DT > 0 (тип: вещественный);
PX - одномерный массив длины N/2 + 1, содержащий вычисленные значения автопериодограммы процесса  X (тип: вещественный);
PY - одномерный массив длины N/2 + 1, содержащий вычисленные значения автопериодограимы процесса  Y (тип: вещественный);
PXY1 - одномерный массив длины N/2 + 1, содержащий вычисленные значения вещественной части взаимной периодограммы процессов X, Y (тип: вещественный);
PXY2 - одномерный массив длины N/2 + 1, содержащий вычисленные значения мнимой части взаимной периодограммы процессов X, Y (тип: вещественный). 
   Версии:  нет 
   Вызываемые подпрограммы
FTF1C - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

После выхода из подпрограммы RSS1R значения исходных случайных процессов в массивах X1, Y1 не сохраняются.

  2. 

B частном случае построения автопериодограммы одного случайного процесса  X в массив Y1 следует заслать нули, при этом допустимы совпадения параметров: PX = PY = PXY1 = PXY2.

  3.  B частном случае построения только автопериодограмм двух случайных процессов X, Y без вычисления их взаимной периодограммы допустимы совпадения параметров: P = PXY1 = PXY2.

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

Unit TRSS1R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, RSS1R_p;

function TRSS1R: String; 

implementation

function TRSS1R: String;
var
N,_i :Integer;
DT :Real;
РХ :Array [0..2] of Real;
PY :Array [0..2] of Real;
PXY1 :Array [0..2] of Real;
PXY2 :Array [0..2] of Real;
const
X1 :Array [0..4] of Real = ( 1.0,2.0,3.0,4.0,0.0 );
Y1 :Array [0..4] of Real = ( 4.0,3.0,2.0,1.0,0.0 );
begin
Result := '';
N := 4;
DT := 1.0;
RSS1R(X1,Y1,N,DT,PX,PY,PXY1,PXY2);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[PX[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[PY[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[PXY1[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[PXY2[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TRSS1R',Result);  { вывод результатов в файл TRSS1R.res }
exit;
end;

end.

Результаты:

       PX  =  (25.,   2.,   1) ,           PY  =  (25.,  2.,  1.) , 

       PXY1  =  (25., - 2., - 1.) ,     PXY2  =  ( 0.,  0.,  0.)