Текст подпрограммы и версий ( Фортран )
rss2r.zip 		Тексты тестовых примеров ( Фортран )
trss2r.zip
Текст подпрограммы и версий ( Си )
rss2r_c.zip 		Тексты тестовых примеров ( Си )
trss2r_c.zip
Текст подпрограммы и версий ( Паскаль )
rss2r_p.zip 		Тексты тестовых примеров ( Паскаль )
trss2r_p.zip

Подпрограмма:  RSS2R

Назначение

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

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

Пусть заданы реализации двух вещественных, стационарных, центрированных случайных процессов 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.

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

    SUBROUTINE  RSS2R (X, Y, N, K, L, DT, SX, SY, SXY1, SXY2, S1,
                                             Z2) 
   Параметры
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 используется в подпрограмме как рабочий.

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

       DIMENSION  X(8), Y(8), SX(3), SY(3), SXY1(3), SXY2(3), Z1(5), 
      *                      Z2(5)
       DATA  X /1., 2., 3., 4., 1., 2., 3., 4./
       DATA  Y /4., 3., 2., 1., 4., 3., 2., 1./
       N = 4
       K = 2
       L = 4
       DT = 1.
       CALL  RSS2R (X, Y, N, K, L, DT, SX, SY, SXY1, SXY2, Z1, Z2)

Результаты:

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

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