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

Назначение

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

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

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

Разобьем каждую реализацию на  K коротких отрезков, возможно перекрывающихся, длины  N: 1 ≤ N ≤ M так, что начальные точки двух соседних отрезков отстоят друг от друга на расстоянии  L: 1 ≤ L ≤ N (т.е. перекрывание коротких отрезков происходит в N - L точках), а значение параметра K ≥ 1 выбирается из соотношения: N + (K - 1) L = M₁ ≤ 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 короткого отрезка. 

Обозначим через A_m^( J ), B_m^( 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, SY_m и SXY1_m = Re SXY_m симметричны относительно точки  m = N/2, а функция SXY2_m = Im SXY_m  антисимметрична относительно  m = N/2.

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

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

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

    int rss2r_c (real *x, real *y, integer *n, integer *k,
            integer *l, real *dt, real *sx, real *sy, real *sxy1, real *sxy2,
            real *z1, real *z2)

   Параметры 

   Версии:  нет 

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

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

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

int main(void)
{
    /* Initialized data */
    static float x[8] = { 1.f,2.f,3.f,4.f,1.f,2.f,3.f,4.f };
    static float y[8] = { 4.f,3.f,2.f,1.f,4.f,3.f,2.f,1.f };

    /* Local variables */
    static float sxy1[3], sxy2[3];
    static int k, l;
    extern int rss2r_c(float *, float *, int *, int *, int *, float *,
                       float *, float *, float *, float *, float *, float *);
    static int n;
    static float z1[5], z2[5], dt, sx[3], sy[3];

    n = 4;
    k = 2;
    l = 4;
    dt = 1.f;
    rss2r_c(x, y, &n, &k, &l, &dt, sx, sy, sxy1, sxy2, z1, z2);

    printf("\n %16.7e %16.7e %16.7e \n", sx[0], sx[1], sx[2]);
    printf("\n %16.7e %16.7e %16.7e \n", sy[0], sy[1], sy[2]);
    printf("\n %16.7e %16.7e %16.7e \n", sxy1[0], sxy1[1], sxy1[2]);
    printf("\n %16.7e %16.7e %16.7e \n", sxy2[0], sxy2[1], sxy2[2]);
    return 0;
} /* main */


Результаты:

      sx  =  (25.,  2.,  1.),             sy  =  (25.,  2.,  1.), 

      sxy1  =  (25.,  - 2.,  - 1.),    sxy2  =  (0.,  0.,  0.)