Текст подпрограммы и версий
amc1r_c.zip , amc1c_c.zip 		Тексты тестовых примеров
tamc1r_c.zip , tamc1c_c.zip

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

Назначение

Вычисление произведения вещественной циркулянтной матрицы на вещественный вектор или произведения двух вещественных циркулянтных матриц на основе алгоритма Быстрого Преобразования Фурье. Подпрограмма может использоваться также для вычисления вещественных периодических сверток.

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

Задана циркулянтная матрица А размерности N*N: 

                   |  a0     aN-1  aN-2   ...   a2  a1    |
                   |  a1     a0     aN-1   ...   a3  a2    |
                   |  a2     a1      a0     ...   a4   a3    |
(1)    A   =   |  . . . . . . . . . . . . . . . . . . . . . . .            ,
                   |  aN-2  aN-3  aN-4   ...  a0  aN-1  |
                   |  aN-1  aN-2  aN-3   ...  a1   a0    |

которая определяется N элементами первого столбца  a0, a1, ..., aN-1, переставляемыми циклически. Требуется вычислить вектор  y = Ax = (y0, y1,..., yN - 1), где вектор  x = (x0, x1,..., xN - 1) также задан. Эта задача сводится к вычислению периодической (круговой) свертки двух последовательностей  as,  xs,  s = 0, 1, ..., N - 1 одинаковой длины N [1], т.е. 

                            N-1 
            ys      =      ∑     as -j x j ,       s = 0, 1, ... , N-1 ,
                             j=0

где ряд as продолжен периодически с периодом N для отрицательных значений индекса  s:  a - s = aN - s ,  s = 1, 2, ..., N - 1.

Пусть F - оператор Дискретного Преобразования Фурье, задаваемый формулой (например, для ряда  as): 

                          N-1
Bm  =  F(as)  =    ∑     as exp (-2 π i m s / N) ,     m = 0, 1,..., N-1 ,   i = (-1)1/2;
                           s=0

аналогично Vm = F(xs), Wm = F(ys),  s = 0, 1, ..., N - 1. На основании теоремы о дискретной круговой свертке [1] 

            Wm = Bm* Vm ,    m = 0, 1, ..., N-1 .

Искомый вектор y находится с помощью Обратного Дискретного Преобразования Фурье F -1 от произведения ДПФ рядов  as и  xs : 

  ys  =  F -1(Wm)  =
                                         N-1
                          =  1/N      ∑    Bm* Vm exp (2 π i m s / N) ,    s = 0, 1,..., N-1,
                                         m=0

Для эффективного вычисления ДПФ и ОДПФ применяется быстрое преобразование Фурье [2], за счет которого число операций в pеализованном алгоритме пропорционально N log2 N и требуется память длины 2N.

Описанный метод можно также применять для быстрого умножения двух циркулянтных матриц: А и матрицы C вида (1), определяемой вектором  x = (x0, x1,..., xN - 1), при этом вектор  y = Ax = (y0, y1,..., yN - 1), определяет циркулянтную матрицу вида (1) произведения A C.

1.  Б.Голд, Ч.Pейдер. Цифровая обработка сигналов. М.,"Советское радио", 1973.
2.  В.А.Морозов, Н.Н.Кирсанова, А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов, Сб. "Численный анализ на ФОPТPАНе", вып.15, Изд - во МГУ, М., 1976.

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

    int amc1r_c (integer *n, real *b, real *x)

   Параметры
n - размерность исходной матрицы A, равная целой степени числа два, n ≥ 4 (тип: целый);
b - вещественный одномерный массив длины n, содержащий заданные определяющие элементы матрицы A:  a0,  a1, ..., an - 1,
x - вещественный одномерный массив длины n, содержащий заданные компоненты вектора x = (x0, x1,..., xn - 1) на входе и вычисленные компоненты вектора  y = (y0, y1,..., yn - 1) на выходе. 
   Версии
amc1c_c - вычисление произведения комплексной циркулянтной матрицы A на комплексный вектор x или произведения двух комплексных циркулянтных матриц. Для подпрограммы amc1c_c параметры b, x имеют тип complex. 
   Вызываемые подпрограммы
 ftf1c_c -
 ftf2c_c   		подпрограммы вычисления Дискретного или Обратного Дискретного Преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье (отличаются друг от друга способом представления исходного ряда); используются в amc1r_c и amc1c_c соответственно. 
   Замечания по использованию
  1. 

После выхода из подпрограммы amc1r_c и amc1c_c значения определяющих элементов исходной матрицы A в массиве b не сохраняются.

  2. 

Если циркулянтная матрица A имеет вид 

              | a0        a1        a2  ...  aN-1   |
              | aN-1     a0        a1  ...  aN-2   |
              | aN-2     aN-1     a0  ...  aN-3   |
              | . . . . . . . . . . . . . . . . . . . . .
              | a1         a2        a3  ...  a0      |

и задается n элементами первой строки  a0,  a1, ..., an - 1, то до обращения к подпрограммам amc1r_c и amc1c_c эти элементы нужно упорядочить в массиве b следующим образом:  a0, an - 1, an - 2, ..., a1.

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


1.
int main(void)
{
    /* Initialized data */
    static float b[8] = { 1.f,8.f,7.f,6.f,5.f,4.f,3.f,2.f };
    static float x[8] = { 1.f,-1.f,-2.f,3.f,0.f,2.f,-3.f,0.f };

    /* Local variables */
    extern int amc1r_c(int *, float *, float *);
    static int i__, n;

    n = 8;
    amc1r_c(&n, b, x);

    for (i__ = 1; i__ <= 8; ++i__) {
         printf("\n  %16.7e \n", x[i__ - 1]);
    }
/* L1: */
    return 0;
} /* main */


Результат:    y  =  (-4., 4., -4., -20., 4., 4., 20., -4.)

2.
int main(void)
{
    /* Initialized data */
    static complex b[4] = { {1.f,1.f},{4.f,-1.f},{3.f,2.f},{2.f,0.f} };
    static complex x[4] = { {1.f,-1.f},{-1.f,0.f},{-2.f,1.f},{3.f,-2.f} };

    /* Local variables */
    extern int amc1c_c(int *, complex *, complex *);
    static int i__, n;

    n = 4;
    amc1c_c(&n, b, x);

    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n  %16.7e %16.7e \n", x[i__ - 1].r, x[i__ - 1].i);
    }
/* L1: */
    return 0;
} /* main */

Результат:  y  =  ( (2., -12.), (11., -4.), (4., -5.), (-3., 3.) )