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

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

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

Подпрограмма: FTA6R (модуль FTA6R_p)

Назначение

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

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

Пусть известны значения 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 имеет два режима работы, задаваемых при обращении к ней значением параметра 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.

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

procedure FTA6R(var ARRAY_ :Array of Real; NDIM :Integer;
                var MESH :Array of Integer; IREG :Integer);

   Параметры

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 -	когда выполняется обратное преобразование.

   Версии

FTA6E -

выполнение прямого и обратного быстрого дискретного преобразования Фурье многомерного массива комплексных чисел, длина каждой размерности которого равняется целой степени двух, в режиме расширенной (Extended) точности; при этом параметр ARRAY_ должен иметь тип Extended.

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

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

	1.	В подпрограммах FTA6R и FTA6E проверка того, что каждое из значений N₁, N₂, ..., N_NDIM должно быть степенью двух, не производится.
	2.	В правиле (I) предполагается, что k_i меняются от 0 до N_i - 1 ( i = 0, 1, 2, ..., NDIM). Поскольку в Фортране индексы массивов не могут начинаться с 0, а только с 1, то при реальной засылке по правилу (I) значений h_{k1,k2,...,kNDIM} в массив ARRAY_ значения k_i должны меняться от 1 до N_i. Это означает, что при такой реальной засылке правило (I) может быть непосредственно применено если при каждом использовании из k_i вычитать 1 (см. также "Пример использования").

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

Unit tfta6r_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, FTA6R_p;

function tfta6r: String; 

implementation

function tfta6r: String;
var
NDIM,N1,N2,K2,K1,I,_i :Integer;
ARRAY_ :Array [0..31] of Real;
MESH :Array [0..1] of Integer;
label
_1;
begin
Result := '';
NDIM := 2;
N1 := 4;
N2 := 4;
MESH[0] := N1;
MESH[1] := N2;
for K2:=1 to N1 do
 begin
  for K1:=1 to N2 do
   begin
    I := Trunc( 1+2*((K1-1)+N1*(K2-1)) );
    ARRAY_[I-1] := 4.0*K1+K2-4.0;
    ARRAY_[I+0] := 0.0;
_1:
   end;
 end;
Result := Result + #$0D#$0A;
for _i:=0 to 31 do
 begin
  Result := Result + Format('%20.16f ',[ARRAY_[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
FTA6R(ARRAY_,NDIM,MESH,1);
Result := Result + #$0D#$0A;
for _i:=0 to 31 do
 begin
  Result := Result + Format('%20.16f ',[ARRAY_[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('tfta6r',Result);  { вывод результатов в файл tfta6r.res }
exit;
end;

end.


Результаты:

    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)