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

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

Назначение

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

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

Заданы две комплексные матрицы  A, B одинаковой размерности N1 * N2 с элементами 

   An1, n2 ,     n1 = 0, 1, ..., N1 - 1 ,     n2 = 0, 1, ..., N2 - 1 ,
   Br1, r2 ,       r1 = 0, 1, ..., N1 - 1 ,      r2 = 0, 1, ..., N2 - 1 .

Требуется вычислить периодическую (круговую) свеpтку  C этих матриц, т.е. 

              N1-1     N2-1
(1)            ∑        ∑     As1 - r1, s2 - r2  Br1, r2  =  Cs1, s2
              r1=0      r2=0

         s1 = 0, 1, ..., N1 - 1 ,    s2 = 0, 1, ..., N2 - 1 ,

где элементы матрицы  A продолжены периодически для отрицательных значений индексов с периодом N1 в каждом столбце и с периодом N2 в каждой стpоке: 

   A- n1, n2  =  AN1 - n1, n2 ,      n1 = 1, 2, ..., N1 - 1 ,   n2 = 0, 1, ..., N2 - 1 ,

   An1, - n2  =  An1, N2 - n2 ,      n1 = 0, 1, ..., N1 - 1 ,   n2 = 1, 2, ..., N2 - 1 .

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

Am1, m2 = F(An1, n2) =

                            N1-1   N2-1
                       =    ∑      ∑      An1, n2 * exp [ - 2π i ( n1*m1/N1 + n2*m2/N2)] ,
                           n1=0   n2=0

         m1 = 0, 1, ..., N1 - 1 ,    m2 = 0, 1, ..., N2 - 1 ,     i = √-1 .

Аналогично   Bm1, m2  =  F (Br1, r2),   Cm1, m2  =  F (Cs1, s2).

По теоpеме о дискретной периодической свеpтке [1] 

(2)       Am1, m2 * Bm1, m2  =  Cm1, m2 ,

          m1 = 0, 1, ..., N1 - 1 ,   m2 = 0, 1, ..., N2 - 1 .

Искомая матрица  С находится путем обратного Дискретного Преобразования Фурье  F - 1 от произведения ДПФ матриц  A и  B: 

  Cs1, s2  =  F - 1 ( Cm1, m2 )  =

                      N1-1  N2-1
= 1/(N1*N2)   ∑     ∑      Am1,m2* Bm1,m2* exp[2π i (s1*m1/N1 + s2*m2/N2)], 
                     m1=0  m2=0

          s1 = 0, 1, ..., N1 - 1 ,     s2 = 0, 1, ..., N2 - 1 .

Подпрограмма предусматривает работу в двух режимах (в зависимости от значения параметра  L).

1.  L=1. Производится вычисление свертки (1) с вычислением ДПФ каждой из двух матриц  A, B.
2.  L=2. Производится вычисление свертки (1) в предположении, что ДПФ матрицы  A уже вычислено. Этот режим позволяет избежать повторного получения комплексной матрицы A при многократном вычислении свертки (1) с той же матрицей  A.

Для эффективного вычисления ДПФ и ОДПФ применяется двумерное Быстрое Преобразование Фурье [2]. При этом время вычислений пропорционально N1 N2 (log2 N1 + log2 N2), а объем памяти пропорционален N1 * N2.

Аналогичный алгоритм для вычисления одномерных периодических сверток двух последовательностей описан в подпрограмме amc1r_c.

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

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

    int rsc2c_c (real *ar, real *ai, real *br, real *bi,
            integer *n1, integer *n2, integer *l)

   Параметры
ar,ai - двумерные вещественные массивы размера n1 * n2, содержащие соответственно вещественные и мнимые части элементов заданной матрицы A;
br, bi - двумерные вещественные массивы размера n1 * n2, содержащие соответственно вещественные и мнимые части элементов заданной матрицы  B; на выходе подпрограммы содержат соответственно вещественные и мнимые части элементов вычисленной матрицы  C;
n1 - заданное число стpок матриц  A и  B,  n1 ≥ 4 - целая степень числа два (тип: целый);
n2 - заданное число столбцов матриц  A и  B,  n2 ≥ 4 - целая степень числа два (тип: целый);
l - параметр, определяющий режим использования подпрограммы (тип: целый);
l=1 - вычисление свертки C матриц  A, B в общем случае,
l=2 - вычисление свертки C матриц  A, B в предположении, что вещественные и мнимые части элементов ДПФ матрицы  A уже содержатся соответственно в массивах ar, ai. 
   Версии:  нет 
   Вызываемые подпрограммы
ftftc_c - вычисление двумерного дискретного или обратного дискретного преобразования Фурье комплексной матрицы размера n1 * n2,  n1 = 2 j1,  n2 = 2 j2 методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

B результате работы подпрограммы rsc2c_c при l = 1 в массивах ar, ai содержатся соответственно вещественные и мнимые части элементов ДПФ матрицы  A.

  2. 

При первом обращении к подпрограмме (l = 1) необходимо задать значения параметров ar, ai, br, bi, n1, n2. При повторном вычислении свертки (1) с той же матрицей  A (l = 2) можно менять только параметры br, bi.

  3. 

Kак следует из формулы (2), сомножители  A, B в свеpтке (1) перестановочны, поэтому при необходимости повторного вычисления свертки с той же матрицей B достаточно при l = 1 матрицу  B задать в массивах ar, ai, а матрицу  A - в массивах br, bi и далее воспользоваться режимом при l = 2.

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

int main(void)
{
    /* Local variables */
    extern int rsc2c_c(float *, float *, float *, float *, int *, int *,
                       int *);
    static int i__, j, l, n1, n2;
    static float ai[32] /* was [4][8] */,
                 bi[32] /* was [4][8] */,
                 ar[32] /* was [4][8] */,
                 br[32] /* was [4][8] */;
    int i__1, i__2;

#define ai_ref(a_1,a_2) ai[(a_2)*4 + a_1 - 5]
#define bi_ref(a_1,a_2) bi[(a_2)*4 + a_1 - 5]
#define ar_ref(a_1,a_2) ar[(a_2)*4 + a_1 - 5]
#define br_ref(a_1,a_2) br[(a_2)*4 + a_1 - 5]

    n1 = 4;
    n2 = 8;
    i__1 = n1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n2;
        for (j = 1; j <= i__2; ++j) {
            ar_ref(i__, j) = (2.f - i__) * (j - 4.f);
            ai_ref(i__, j) = 0.f;
            br_ref(i__, j) = 6.f - i__ - j;
/* l1: */
            bi_ref(i__, j) = (float) (i__ - j);
        }
    }
    l = 1;
    rsc2c_c(ar, ai, br, bi, &n1, &n2, &l);

    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f \n",
         br_ref(i__, 1), br_ref(i__, 2), br_ref(i__, 3), br_ref(i__, 4),
         br_ref(i__, 5), br_ref(i__, 6), br_ref(i__, 7), br_ref(i__, 8));
    }
    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f \n",
         bi_ref(i__, 1), bi_ref(i__, 2), bi_ref(i__, 3), bi_ref(i__, 4),
         bi_ref(i__, 5), bi_ref(i__, 6), bi_ref(i__, 7), bi_ref(i__, 8));
    }
    return 0;
} /* main */


Результаты:
 
                | - 16.  24.  48.  56.  48.  24. - 16. - 72. |
     br  =    | - 8.   32.  56.  64.  56.  32.  - 8.  - 64. |
                | - 16.  24.  48.  56.  48.  24. - 16. - 72. |
                | - 40.  0.    24.  32.  24.  0.   - 40. - 96. |

               | - 16.  24.  48.  56.  48.  24. - 16. - 72. |
     bi  =   | - 24.  16.  40.  48.  40.  16. - 24. - 80. |
               | - 16.  24.  48.  56.  48.  24. - 16. - 72. |
               |   8.    48.  72.  80.  72.  48.    8.   - 48. |