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

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

Назначение

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

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

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

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

    int rss1r_c (real *x1, real *y1, integer *n, real *dt,
            real *px, real *py, real *pxy1, real *pxy2)

   Параметры
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_c - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

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

  2. 

B частном случае построения автопериодограммы одного случайного процесса  X в массив y1 следует заслать нули, при этом допустимы совпадения параметров: px = py = pxy1 = pxy2.

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

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

int main(void)
{
    /* Initialized data */
    static float x1[5] = { 1.f,2.f,3.f,4.f,0.f };
    static float y1[5] = { 4.f,3.f,2.f,1.f,0.f };

    /* Local variables */
    static float pxy1[3], pxy2[3];
    extern int rss1r_c(float *, float *, int *, float *, float *,
                       float *, float *, float *);
    static int n;
    static float dt, px[3], py[3];

    n = 4;
    dt = 1.f;
    rss1r_c(x1, y1, &n, &dt, px, py, pxy1, pxy2);

    printf("\n %16.7e %16.7e %16.7e \n", px[0], px[1], px[2]);
    printf("\n %16.7e %16.7e %16.7e \n", py[0], py[1], py[2]);
    printf("\n %16.7e %16.7e %16.7e \n", pxy1[0], pxy1[1], pxy1[2]);
    printf("\n %16.7e %16.7e %16.7e \n", pxy2[0], pxy2[1], pxy2[2]);
    return 0;
} /* main */


Результаты:

       px  =  (25.,   2.,   1.) ,           py  =  (25.,  2.,  1.) , 

       pxy1  =  (25., - 2., - 1.) ,     pxy2  =  ( 0.,  0.,  0.)