БЧА НИВЦ МГУ. FTA6R_C. Ускоренные дискретные преобразования Фурье

Текст подпрограммы и версий
fta6r_c.zip , fta6d_c.zip

Тексты тестовых примеров
tfta6r_c.zip , tfta6d_c.zip

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

Назначение

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

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

Пусть известны значения h_{k1,k2,...,kNDIM} = h (t_1,k1, t_2,k2, ..., t_NDIM,kNDIM) комплексно - значной функции h (t₁, t₂, ..., t_NDIM) от NDIM вещественных переменных на многомерной сетке размерности NDIM, узлы которой по каждой размерности (т.е. по каждой переменной) определяются следующим образом:

t_1,k1 = k1 δt₁ ( k1 = 0, 1, 2, ..., N₁ - 1; δt₁ - шаг сетки по первой размерности; N₁ равняется целой степени двух и определяет число узлов сетки по первой размерности)

t_2,k2 = k2 δt₂ ( k2 = 0, 1, 2, ..., N₂ - 1; δt₂ - шаг сетки по второй размерности; N₂ равняется целой степени двух и определяет число узлов сетки по второй размерности )
   .
   .
   .

t_NDIM,kNDIM = kNDIM δt_NDIM ( kNDIM = 0, 1, 2, ..., N_NDIM ; δt_NDIM - шаг сетки по последней размерности NDIM; N_NDIM равняется целой степени двух и определяет число узлов сетки по последней размерности NDIM )

Эти комплексные значения h_{k1,k2,...,kNDIM} задаются в виде одномерного вещественного массива ARRAY длины 2 * N₁ * N₂ * ... * N_NDIM (здесь * означает знак умножения) по следующему правилу (I)
(Re означает вещественную часть, а Im - мнимую часть комплексного числа):

   Re(h_{k1,k2,...,kNDIM}) → ARRAY( 1 + 2(k1 + N₁*k2 + N₁*N₂*k3 + ... +
                                                     + N₁*N₂*...*N_NDIM-1*kNDIM) )

   Im(h_{k1,k2,...,kNDIM}) → ARRAY( 2 + 2(k1 + N₁*k2 + N₁*N₂*k3 + ... +
                                                     + N₁*N₂*...*N_NDIM-1*kNDIM) )

Таким образом, исходные комплексные значения представляются в виде вещественного одномерного массива по правилам Фортрана, определенным для представления комплексных многомерных массивов в виде вещественных одномерных массивов (см. также "Замечания по использованию").

Информация о количестве узлов NDIM - мерной сетки по каждой размерности указывается в компонентах целого массива MESH длины NDIM следующим образом:

               MESH  =  (N₁, N₂, ..., N_NDIM)

Подпрограмма fta6r_c имеет два режима работы, задаваемых при обращении к ней значением параметра IREG. В первом режиме (IREG = 1) выполняется NDIM - мерное прямое быстрое дискретное преобразование Фурье заданного указанным выше образом массива ARRAY, состоящее в получении N₁ * N₂ * ... * N комплексных чисел

H_{n1,n2,...,nNDIM} ( n1 = 0, 1, 2, ..., N₁ - 1, n2 = 0, 1, 2, ..., N₂ - 1, ..., nNDIM = 0, 1, 2, ..., N_NDIM - 1 ) по формуле:

H_{n1,n2,...,nNDIM}  =

                         N_NDIM-1        N₁-1
                     =       ∑          ...    ∑       exp(2πi * kNDIM * nNDIM / N_NDIM) *
                           ^kNDIM=0       ^k1=0

                                                * ... * exp(2πi * k1 * n1 / N₁) h_{k1,k2,...,kNDIM}

Вычисленные значения H_{n1,n2,...,nNDIM} располагаются в том же массиве ARRAY по правилу (I).

Во втором режиме (IREG = - 1) выполняется NDIM - мерное обратное быстрое дискретное преобразование Фурье, состоящее в восстановлении из значений H_{n1,n2,...,nNDIM}, заданных по правилу (I) в одномерном вещественном массиве ARRAY длины 2 * N₁ * N₂ * ... * N_NDIM, значений h_{k1,k2,...,kNDIM} по формуле

h_{k1,k2,...,kNDIM}  =  1 / (N₁*N₂*...*N_NDIM) *

                          N_NDIM-1      N₁-1
                      *      ∑          ...   ∑      exp(- 2πi * nNDIM * kNDIM / N_NDIM) *
                         ^nNDIM=0        ⁿ¹⁼⁰

                                               * ... * exp(- 2πi * n1 * k1 / N₁) H_{n1,n2,...,nNDIM}

Вычисленные значения h_{k1,k2,...,kNDIM} располагаются в том же массиве ARRAY по правилу (I).

Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков. Численные методы. Изд - во "Наука", 1987.

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

    int fta6r_c (real *array, integer *ndim, integer *mesh,
            integer *ireg)

   Параметры

array -	вещественный одномерный массив длины 2 * n₁ * n₂ * ... * n_ndim, содержащий на входе h_{k1,k2,...,kndim} (ireg = 1) или h_{n1,n2,...,nndim} (ireg = - 1), а на выходе значения h_{n1,n2,...,nndim} (ireg = 1) или h_{k1,k2,...,kndim} (ireg = - 1);
ndim -	число размерностей заданной многомерной сетки (тип: целый);
mesh -	целый массив, i - ая компоненты которого содержит число узлов заданной ndim - мерной сетки по i - му направлению ( i = 1, 2, 3, ..., ndim), т.е. mesh = (n₁, n₂, ..., n_ndim); каждое значение n_i должно быть целой степенью двух;
ireg -	задает режим работы подпрограммы (тип: целый); при этом:

ireg= 1 -	когда выполняется прямое преобразование;
ireg= -1 -	когда выполняется обратное преобразование

   Версии

fta6d_c -

выполнение прямого и обратного быстрого дискретного преобразования Фурье многомерного массива комплексных чисел, длина каждой размерности которого равняется целой степени двух, в режиме удвоенной точности; при этом параметр array должен иметь тип double

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

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

В подпрограммах fta6r_c и fta6d_c проверка того, что каждое из значений n₁, n₂, ..., n_ndim должно быть степенью двух, не производится.

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

int main(void)
{
    int i__1, i__2;

    /* Local variables */
    static int ndim, mesh[2];
    extern int fta6r_c(float *, int *, int *, int *);
    static int i;
    static float array[32];
    static int k1, k2, n1, n2;
    static int c__1 = 1;

    ndim = 2;
    n1 = 4;
    n2 = 4;
    mesh[0] = n1;
    mesh[1] = n2;
    i__1 = n1;
    for (k2 = 1; k2 <= i__1; ++k2) {
        i__2 = n2;
        for (k1 = 1; k1 <= i__2; ++k1) {
            i = (k1 - 1 + n1 * (k2 - 1) << 1) + 1;
            array[i - 1] = k1 * 4.f + k2 - 4.f;
            array[i] = 0.f;
/* l1: */
        }
    }
    for (i = 0; i <= 28; i+= 4) {
        printf("\n  %15.6e %15.6e %15.6e %15.6e \n",
               array[i], array[i+1], array[i+2], array[i+3]);
    }

    fta6r_c(array, &ndim, mesh, &c__1);

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


Результаты:

    array = (136, 0, - 32, - 32, - 32, 0, - 32, 32, - 8, - 8, 0, 0, 0, 0,
                       0, 0, - 8, 0, 0, 0, 0, 0, 0, 0, - 8, 8, 0, 0, 0, 0, 0, 0)