Подпрограмма:  FTA5R (модуль FTA5R_p)

Назначение

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

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

Пусть известны значения  f _j = f ( t _j ) вещественной функции  f (t) на равномерной сетке  t _j = j Δt,  j = 0, 1, 2, ..., N - 1 (где  Δt - шаг сетки, а  N равняется целой степени двух), заданные в виде одномерного массива вещественных чисел ARRAY длины  N.

Подпрограмма FTA5R имеет два режима работы, задаваемых при обращении к ней значением параметра IREG. В первом режиме ( IREG = 1) выполняется прямое быстрое дискретное косинус - преобразование Фурье заданного массива ARRAY, состоящее в получении  N вещественных чисел  F_k (k = 0, 1, 2, ..., N - 1) из чисел  f _j по формуле:

                              N-1
                  Fk   =   ∑     f j cos(π j k / N) .
                              j =1 

В этом режиме подпрограмма FTA5R из исходного массива ARRAY строит вспомогательный вещественный массив  Y длины  N, компоненты которого  y_j ( j = 0, 1, 2, .., N - 1) определяются формулами:

             y 0  =  f 0 
             y j  =  0.5 ( f j + f N-j ) - sin( j π / N ) ( f j - f N-j ) ,     j > 0 . 

В результате применения к полученному массиву  Y подпрограммы FTA3R строится комплексный массив длины N/2 с вещественными R_k и мнимыми  I_k (k = 0, 1, 2, ..., N/2 - 1) частями. Компоненты искомого массива определяются через R_k и  I_k следующим образом:

                          N-1
               F1   =   ∑    f j cos ( j π / N )
                          j= 0 
               F2k  =  Rk ,    F2k+1  =  F2k-1 + Ik ,    k = 0, 1, 2, ..., N/2-1 

Вычисленные значения  F_k располагаются в том же массиве ARRAY.

Во втором режиме ( IREG = - 1) выполняется обратное быстрое дискретное косинус - преобразование Фурье, состоящее в восстановлении из значений  F_k (k = 0, 1, 2, ..., N - 1), заданных в виде одномерного массива вещественных чисел ARRAY длины  N, значений  f _j ( j = 0, 1, 2, ..., N - 1), которые располагаются в том же массиве ARRAY. Для этого строится массив  f _m ( m = 0, 1, 2, ..., N - 1) по следующей формуле:

                          N-1
               f m  =   ∑    Fk cos ( π k m / N ) .
                          k =0 
   Между значениями  f m и  f m существуют соотношения:

                f 0   =   N f0  +  ∑ f j
                                        по нечетным  j

(1)            f m  =  (N/2) f m  +  ∑ f j    ,            m - нечетные
                                            по четным  j

                f m  =  (N/2) f m  +  ∑ f j    ,             m - четные ,       m ≠ 0
                                            по нечетным  j  

Неизвестные суммы в правых частях этих соотношений определяются следующим образом:

                                                                N-1
                C1 ≡  ∑  f m =  (N/2) ( f 0  +     ∑    f j )
                                                                j =0
                      по четным  m

                                                      N-1
                C2 ≡  ∑  f m =  (N/2)     ∑    f j
                                                      j =0
                      по нечетным  m

  Отсюда следует, что

                (N/2) f 0  =  C1 - C2 .
  Поэтому

                  ∑ f j  =   f 0  - 2 ( C1 - C2 )
          по нечетным  j

                  ∑ f j  =   (2/N) C2 - ∑ f j
          по четным  j             по нечетным  j 

Зная эти суммы, восстанавливаемые значения  f _j вычисляются из соотношений (1).

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

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

procedure FTA5R(var ARRAY_ :Array of Real; N :Integer; IREG :Integer); 

   Параметры 

   Версии 

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

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

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

Unit TFTA5R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, FTA5R_p;

function TFTA5R: String; 

implementation

function TFTA5R: String;
var
N,_i :Integer;
const
ARRAY_ :Array [0..7] of Real = ( 1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0 );
begin
Result := '';
N := 8;
FTA5R(ARRAY_,N,1);
Result := Result + #$0D#$0A;
for _i:=0 to 7 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;
FTA5R(ARRAY_,N,-1);
Result := Result + #$0D#$0A;
for _i:=0 to 7 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('TFTA5R',Result);  { вывод результатов в файл TFTA5R.res }
exit;
end;

end.

Результаты: 

   а) в случае прямого преобразования:
      ARRAY_ = (36, - 8.13707, - 4, 3.38009, - 4, 4.27677, - 4, 4.48022) 

   б) в случае обратного преобразования:
      ARRAY_ = (1, 2, 3, 4, 5, 6, 7, 8)