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

Назначение

Вычисление прямой и обратной дискретной периодической (круговой) свертки двух одномерных массивов вещественных чисел.

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

Пусть известны  N значений  s_j = s ( t_j ) функции сигнала  s ( t ) на равномерной сетке  t_j = j Δt,  j = 0, 1, 2, ..., N - 1 (где  Δt - шаг сетки, а  N равняется целой степени двух). Предполагается, что сигнал  s_j является периодическим с периодом  N (отсюда, в частности, следует, что значения  s_j известны и для отрицательных значений  j, т.е.  s_- j = s_N - j,  j = 1, 2, ..., N - 1). Значения  s_j задаются в вещественном одномерном массиве ARRAY длины  N:  ARRAY ( J ) = s_j ( J = 1, 2, ..., N;  j = 0, 1, 2, ..., N - 1).

Пусть известны также  M значений  r_k = r ( t_k ) импульсной функции отклика  r ( t ) на той же равномерной сетке  t_k = k Δt,  k = - (M - 1)/2, ..., 0, ..., (M - 1)/2 (где  Δt - тот же шаг сетки, что и для функции сигнала, а  M - нечетное натуральное число такое, что M < N). Значения  r_k задаются в  M первых элементах одномерного вещественного массива RESP длины  N следующим образом:

       RESP(1)  =  r0
      RESP( I )  =  rk ,     I = 2, 3, ..., (M+1)/2 ;   k = 1, 2, ..., (M-1)/2
      RESP( I )  =  rk ,     I = (M+3)/2, ..., M ;      k = (M-1)/2, ..., - 1 

Подпрограмма rsa1r_c имеет два режима работы, задаваемых при обращении к ней значением параметра IREG. В первом режиме ( IREG = 1) выполняется прямая дискретная периодическая свертка массивов ARRAY и RESP (L = 1, 2, ..., N):

                               N
             (r*s)L   =   ∑   ARRAY(L - J) × RESP( J )
                              j =1 

где  * - знак свертки, а  × - знак умножения.
Предварительно в массив RESP добавляются нули следующим образом:

   RESP( I )  =  RESP( I )                 ,   I = 1, 2, ..., (M+1)/2 
   RESP( I )  =  0.0                           ,   I = (M+3)/2, ..., N - (M-1)/2 
   RESP(N+1-I )  =  RESP(M+1-I ) ,   I = 1, 2, ..., (M-1)/2 

Для вычисления указанной свертки используется быстрое дискретное преобразование Фурье в соответствии с теоремой о дискретной периодической (круговой) свертке. Полученные значения (r * s)_L помещаются в первых  N компонентах массива ANS.

Во втором режиме ( IREG = - 1) выполняется восстановление значений  s_j функции сигнала по известным значениям  r_k импульсной функции ответа, задаваемым в массиве RESP, и значениям соответствующей свертки, задаваемым в массиве ARRAY. Полученные значения  s_j помещаются в первых  N элементах массива ANS.

Р.Отнес, Л.Эноксон. Прикладной анализ временных рядов. Изд - во "Мир", 1982.

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

    int rsa1r_c (real *array, integer *n, real *resp, integer *m,
        complex *ans, complex *fft, integer *ireg, integer *iflag)

   Параметры 

   Версии 

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

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

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

int main(void)
{
    /* Initialized data */
    static float array[16] = { 0.f,0.f,.5f,0.f,0.f,-.5f,0.f,0.f,.5f,.5f,.5f,
                                .5f,.5f,0.f,0.f,0.f };
    /* Local variables */
    static float resp[16];
    extern int rsa1r_c(float *, int *, float *, int *, float *, complex *,
                       int *, int *);
    static int m, n, iflag, i;
    static complex fft[16];
    static float ans[32];

    n = 16;
    m = 3;
    resp[0] = .7f;
    resp[1] = .4f;
    resp[2] = .4f;
    rsa1r_c(array, &n, resp, &m, ans, fft, &c__1, &iflag);

    printf("\n %5i \n", iflag);
    for (i = 0; i <= 28; i += 4) {
         printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                ans[i], ans[i+1], ans[i+2], ans[i+3]);
    }
    return 0;
} /* main */


Результаты:

       ans = (- 0.596046e - 7, 0.2, 0.35, 0.2, - 0.2, - 0.35, - 0.2, 
                       0.2, 0.55, 0.75, 0.75, 0.75, 0.55, 0.2, - 0.745058e - 8,  
                       0.521541e - 07)