Текст подпрограммы и версий
ams1r_p.zip , ams1c_p.zip 		Тексты тестовых примеров
tams1r_p.zip , tams1c_p.zip

Подпрограмма:  AMS1R (модуль AMS1R_p)

Назначение

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

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

1. Задана теплицева матрица А размерности N*N: 

                   | a0      a-1    a-2    ...  a-(N-2)   a-(N-1)  |
                   | a1      a0     a-1    ...  a-(N-3)   a-(N-2)  |
                   | a2      a1      a0    ...  a-(N-4)   a-(N-3)  |
(1)              | . . . . . . . . . . . . . . . . . . . . . . . . . . . .           ,
                   | aN-2  aN-3  aN-4  ...  a0         a-1        |
                   | aN-1  aN-2  aN-3  ...  a1         a0         |

которая определяется 2N - 1 элементами первой строки и первого столбца: 

         a -(N-1), a -(N-2), ... , a -2, a -1, a0, a1, a2, ... , aN-2, aN-1 ;

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

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

Дополним теплицеву матрицу А порядка N до циркулянтной матрицы А0 порядка М, которая определяется М элементами своего первого столбца: 

         a -(N-1), a -(N-2), ... , a -1, a0, a1, ... , aN-1, 0, 0, ... , 0 ,

где М - наименьшая целая степень числа два, удовлетворяющая условию М ≥ 2N - 1, а число нулей K равно разности: К = М - (2N - 1) ≥  0.

Далее определим вектор  x0 длины М, получаемый из заданного вектора  x добавлением К + N - 1 нулей: x0 = (x0, x1,..., xN - 1, 0, 0,..., 0). Тогда исходная задача сводится к задаче умножения циркулянтной матрицы А0 на вектор  x0, которая решается в подпрограмме AMC1R на основе алгоритма Быстрого Преобразования Фурье [2], при этом N компонент вычисленного вектора y0 = А0 x0 = (y00, y10,..., y0M - 1) будут совпадать с искомым вектором  y, а именно:  y0N - 1 + s = ys ,  s = 0, 1, ..., N - 1.

Аналогичный метод умножения теплицевой матрицы на вектор изложен в работе [3].

Реализованный в подпрограмме алгоритм требует памяти ЭВМ не менее 4N и не более 8N (из которых 3N занимает исходная информация). Время вычислений пропорционально М log2 М, т.е. по порядку не превосходит 4N (2 + log2N).

2. Пусть А - ханкелева матрица вида 

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

которая определяется 2N - 1 элементами первой строки и последнего столбца: 

         a0, a1, ... , aN-1, aN, ... , a2N-3, a2N-2 ;

все элементы этой матрицы на любой ее диагонали, параллельной побочной, равны между собой. В данном случае задача вычисления вектора y = Ax = (y0, y1,..., yN - 1) эквивалентна вычислению взаимной корреляции двух последовательностей  as,  xj,  s = 0, 1,..., 2N - 2,  j = 0, 1,..., N - 1 длины 2N - 1 и N соответственно [1], т.е. 

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

в частности, при  xj = aj,  j = 0, 1,..., N - 1 вектор  y = Аx определяет автокорреляцию последовательности  as,  s = 0, 1, ..., 2N - 2.

Легко видеть, что ханкелева матрица А переходит в теплицеву T путем умножения слева на матрицу перестановок U вида 

                               |                   1 |
                               |                1    |
                               |             .        |
                     U =    |          .           |    =   U -1 ,
                               |       .              |
                               |    1                |
                               | 1                   |

т.е. Т = UА. Поэтому для получения вектора  y = Аx сначала вычисляется вектор Tx (по описанному выше алгоритму) и затем делается соответствующая перестановка его компонент, а именно:  y = U(Тx).

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

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

procedure AMS1R(N :Integer; var M :Integer; Q :Real;
                var B :Array of Real; var X :Array of Real); 

   Параметры
N - размерность исходной матрицы А, N ≥ 2 (тип: целый);
M - минимальная целая степень числа два, удовлетворяющая условию М ≥ 2N - 1 (тип: целый);
Q - управляющий параметр (тип: вещественный):
Q= 0. - если исходная матрица А теплицева,
Q= 1. - если исходная матрица А ханкелева;
B - вещественный одномерный массив длины М, в первых 2N - 1 элементах которого содержатся заданные определяющие элементы теплицевой матрицы А: a- (N - 1),  a- (N - 2),..., a- 1, a0, a1, ..., aN - 1 или ханкелевой матрицы А: a0, a1, ..., aN - 1, aN, ..., a2N - 3, a2N - 2;
X - вещественный одномерный массив длины М, в первых N элементах которого содержатся заданные компоненты вектора x = (x0, x1,..., xN - 1) на входе и вычисленные компоненты вектора y = (y0, y1,..., yN - 1), на выходе. 
   Версии
AMS1C - вычисление произведения комплексной теплицевой или ханкелевой матрицы А на комплексный вектор  x. Для подпрограммы АМS1С параметры B, X имеют тип Complex. 
   Вызываемые подпрограммы
 AMC1R -
 AMC1C   		подпрограммы умножения вещественной циркулянтной матрицы на вещественный вектор и компексной циркулянтной матрицы на комплексный вектор на основе алгоритма Быстрого Преобразования Фурье (используются в АМS1R и АМS1С соответственно). 
   Замечания по использованию
  1. 

После выхода из подпрограмм АМS1R и АМS1С значения определяющих элементов исходной матрицы А в массиве B не сохраняются.

  2. 

Если матрица А циркулянтная (частный случай теплицевой),то целесообразно вместо подпрограмм АМS1R, АМS1С воспользоваться соответствующими подпрограммами АМС1R и АMC1C.

  3. 

Если теплицева матрица А имеет вид 

          | a0         a1        a2         ...   aN-1  |
          | a-1        a0        a1         ...   aN-2  |
          | a-2        a-1       a0         ...   aN-3  |
          | . . . . . . . . . . . . . . . . . . . . . . . . .
          | a-(N-1)   a-(N-2)  a-(N-3)  ...  a0     |
и задается 2N - 1 элементами первого столбца и первой строки:  a- (N - 1), a- 1, a0, a1, a2,..., aN - 1; то до обращения к подпрограммам АМS1R, АМS1С эти элементы следует упорядочить в массиве B по убыванию индексов:  aN - 1, ..., a1, a0, a- 1, a- 2, ..., a- (N - 1) .

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

   1.

Unit tams1r_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, AMS1R_p;

function tams1r: String; 

implementation

function tams1r: String;
var
N,M,I :Integer;
Q :Real;
const
B :Array [0..15] of Real = ( 1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,
0.0,0.0,0.0,0.0,0.0 );
X :Array [0..15] of Real = ( 1.0,-1.0,-2.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0 );
label
_1;
begin
Result := '';
N := 6;
M := 16;
Q := 1.0;
AMS1R(N,M,Q,B,X);
for I:=1 to N do
 begin
_1:
  Result := Result + Format('                    %20.16f ',
[X[I-1]]) + #$0D#$0A;
 end;
UtRes('tams1r',Result);  { вывод результатов в файл tams1r.res }
exit;
end;

end.


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

   2.
Unit tams1c_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AMS1C_p;

function tams1c: String;

implementation

function tams1c: String;
var
N,M,I :Integer;
Q :Real;
const
B :Array [0..7] of Complex = ( ( re:4.0; im:-1.0 ),( re:3.0; im:2.0 ),( 
re:2.0; im:0.0 ),( re:1.0; im:1.0 ),( re:-1.0; 
im:0.0 ),( re:0.0; im:-1.0 ),( re:2.0; im:1.0 ),( 
re:0.0; im: 0.0 ) );
X :Array [0..7] of Complex = ( ( re:1.0; im:-1.0 ),( re:-1.0; im:0.0 ),( 
re:-2.0; im:1.0 ),( re:3.0; im:-2.0 ),( re:0.0; 
im:0.0 ),( re:0.0; im:0.0 ),( re:0.0; im:0.0 ),( 
re:0.0; im: 0.0 ) );
label
_1;
begin
Result := '';  { результат функции }
N := 4;
M := 8;
Q := 0.0;
AMS1C(N,M,Q,B,X);
for I:=1 to N do
 begin
_1:
  Result := Result + Format('          %20.16f %20.16f ',
 [X[I-1].re,X[I-1].im]) + #$0D#$0A;
 end;
UtRes('tams1c',Result);  { вывод результатов в файл tams1c.res }
exit;
end;

end.

Результат:     y  =  ( (2., -12.), (7., 2.), (3., -6.), (10., 0.) )