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

Подпрограмма:  RSS2R (модуль RSS2R_p)

Назначение

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

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

Пусть заданы реализации двух вещественных, стационарных, центрированных случайных процессов Xq, Yq,  q = 0, 1, ..., M - 1 достаточно большой длины  M на сетке по оси времени  tq = q * DT с шагом DT > 0.

Разобьем каждую реализацию на  K коротких отрезков, возможно перекрывающихся, длины  N: 1 ≤ N ≤ M так, что начальные точки двух соседних отрезков отстоят друг от друга на расстоянии  L: 1 ≤ L ≤ N (т.е. перекрывание коротких отрезков происходит в N - L точках), а значение параметра K ≥ 1 выбирается из соотношения: N + (K - 1) L = M1 ≤ M.
Вообще говоря, не требуется, чтобы конец последнего короткого отрезка длины N совпадал с концами исходных рядов.
Таким образом, получим отрезки временны'х рядов вида: 

     Xs( J )  =  Xs + ( J-1) L ,     Ys( J )  =  Ys + ( J-1) L ,

        s = 0, 1, ..., N - 1 ,    J = 1, 2, ..., K ,

 где  J - номеp короткого отрезка.

Обозначим через Am( J ), Bm( J ),  m = 0, 1, ..., N - 1 их соответствующие дискретные преобразования Фурье, которые вычисляются по формулам: 

                        N-1
      Am( J )   =    ∑   Xs( J ) exp(- i λm ts )   =
                        s=0
                                                                   N-1   
                                                               =   ∑  Xs( J ) exp(- 2π i m s / N) ,
                                                                   s=0
                        N-1
      Bm( J )   =    ∑   Ys( J ) exp(- i λm ts )   =
                        s=0
                                                                   N-1   
                                                               =   ∑  Ys( J ) exp(- 2π i m s / N) ,
                                                                   s=0

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

B качестве оценок для автоспектров SX, SY и взаимного спектра SXY исходных случайных процессов берутся периодограммы, усредненные по  K указанным коротким отрезкам, а именно, величины: 

                            K
      SXm = 1/ K   ∑   (DT / N)  | Am( J ) | 2  ,
                           J=1
                            K
      SYm = 1/ K   ∑   (DT / N)  | Bm( J ) | 2  ,            m = 0, 1, ..., N - 1  , 
                           J=1
                               K
      SXYm = 1/ K    ∑  (DT / N)  Am* ( J ) Bm( J )  , 
                              J=1

где * - знак комплексного сопряжения. При этом функции SХm, SYm и SXY1m = Re SXYm симметричны относительно точки  m = N/2, а функция SXY2m = Im SXYm  антисимметрична относительно  m = N/2.

Данный метод модифицированных периодограмм [1] дает уменьшение дисперсии спектральных оценок порядка  K раз (при  N << M) по сравнению с оценками через периодограммы без усреднения.

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

1.  P.D.Welch, The Use of Fast Fourier Transform for the Estimation of Power Spectra, IEEE Transactions on Audio, Electroacoustics, June, 1967, AU - 15, N 2, 70 - 73.
2.  М.В.Арефьева, Корреляционный и спектральный анализ стационарных случайных процессов (часть 1), сб."Численный анализ на ФОРТРАНе", вып.15. Изд - во МГУ, M., 1976.

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

procedure RSS2R(var X :Array of Real; var Y :Array of Real;
                var N :Integer; K :Integer; L :Integer; DT :Real;
                var SX :Array of Real; var SY :Array of Real;
                var SXY1 :Array of Real; var SXY2 :Array of Real;
                var Z1 :Array of Real; var Z2 :Array of Real);


   Параметры
X, Y - одномерные массивы длины M1 = N + (K - 1) * L, содеpжащие последовательные значения исходных реализаций случайных процессов Xq, Yq,  q = 0, 1, ..., M1 - 1 (тип: вещественный);
N - заданная длина короткого отрезка, N ≥ 4 - целая степень числа два (тип: целый);
K - заданное число коротких отрезков (тип: целый);
L - заданное расстояние между начальными точками двух соседних коротких отрезков, 1 ≤ L ≤ N (тип: целый);
DT - заданный шаг сетки по оси времени, DT > 0 (тип: вещественный);
SX - одномерный массив длины N/2 + 1, содержащий вычисленные значения автоспектра процесса  X (тип: вещественный);
SY - одномерный массив длины N/2 + 1, содержащий вычисленные значения автоспектра процесса  Y (тип: вещественный);
SXY1 - одномерный массив длины N/2 + 1, содержащий вычисленные значения вещественной части взаимного спектра процессов X, Y - коспектра (тип: вещественный);
SXY2 - одномерный массив длины N/2 + 1, содержащий вычисленные значения мнимой части взаимного спектра процессов X, Y - квадратурного спектра (тип: вещественный);
Z1, Z2 - одномерные массивы длины N + 1, используемые в подпрограмме как рабочие (тип: вещественный). 
   Версии:  нет 
   Вызываемые подпрограммы
FTF1C - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

Значения параметров K, L, N должны быть связаны соотношением: N + (K - 1) * L = M1 ≤ M, где  M - общее число заданных значений каждого из двух случайных процессов Xq, Yq, M1 - количество этих значений, используемых в подпрограмме для построения оценок спектров. В частности, возможно M = M1.

  2.  B частном случае одного случайного процесса  X допустимы совпадения параметров: X = Y,  SY = SXY1 = SXY2, при этом массив SY используется в подпрограмме как рабочий.

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

Unit TRSS2R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, RSS2R_p;

function TRSS2R: String; 

implementation

function TRSS2R: String;
var
N,K,L,_i :Integer;
DТ :Real;
SX :Array [0..2] of Real;
SУ :Array [0..2] of Real;
SXY1 :Array [0..2] of Real;
SXY2 :Array [0..2] of Real;
Z1 :Array [0..4] of Real;
Z2 :Array [0..4] of Real;
const
X :Array [0..7] of Real = ( 1.0,2.0,3.0,4.0,1.0,2.0,3.0,4.0 );
Y :Array [0..7] of Real = ( 4.0,3.0,2.0,1.0,4.0,3.0,2.0,1.0 );
begin
Result := '';
N := 4;
K := 2;
L := 4;
DT := 1.0;
RSS2R(X,Y,N,K,L,DT,SX,SY,SXY1,SXY2,Z1,Z2);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[SX[_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 ',[SY[_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 ',[SXY1[_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 ',[SXY2[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TRSS2R',Result);  { вывод результатов в файл TRSS2R.res }
exit;
end;

end.

Результаты:

      SX  =  (25.,  2.,  1.),             SY  =  (25.,  2.,  1.), 

      SXY1  =  (25.,  - 2.,  - 1.),    SXY2  =  (0.,  0.,  0.)