|
Текст подпрограммы и версий amc1r_p.zip , amc1c_p.zip |
Тексты тестовых примеров tamc1r_p.zip , tamc1c_p.zip |
Вычисление произведения вещественной циркулянтной матрицы на вещественный вектор или произведения двух вещественных циркулянтных матриц на основе алгоритма Быстрого Преобразования Фурье. Подпрограмма может использоваться также для вычисления вещественных периодических сверток.
Задана циркулянтная матрица А размерности N*N:
| a0 aN-1 aN-2 ... a2 a1 |
| a1 a0 aN-1 ... a3 a2 |
| a2 a1 a0 ... a4 a3 |
(1) A = | . . . . . . . . . . . . . . . . . . . . . . . ,
| aN-2 aN-3 aN-4 ... a0 aN-1 |
| aN-1 aN-2 aN-3 ... a1 a0 |
которая определяется N элементами первого столбца a0, a1, ..., aN-1, переставляемыми циклически. Требуется вычислить вектор y = Ax = (y0, y1,..., yN - 1), где вектор x = (x0, x1,..., xN - 1) также задан. Эта задача сводится к вычислению периодической (круговой) свертки двух последовательностей as, xs, s = 0, 1, ..., N - 1 одинаковой длины N [1], т.е.
N-1
ys = ∑ as -j x j , s = 0, 1, ... , N-1 ,
j=0
где ряд as продолжен периодически с периодом N для отрицательных значений индекса s: a - s = aN - s , s = 1, 2, ..., N - 1.
Пусть F - оператор Дискретного Преобразования Фурье, задаваемый формулой (например, для ряда as):
N-1
Bm = F(as) = ∑ as exp (-2 π i m s / N) , m = 0, 1,..., N-1 , i = (-1)1/2;
s=0
аналогично Vm = F(xs), Wm = F(ys), s = 0, 1, ..., N - 1. На основании теоремы о дискретной круговой свертке [1]
Wm = Bm* Vm , m = 0, 1, ..., N-1 .
Искомый вектор y находится с помощью Обратного Дискретного Преобразования Фурье F -1 от произведения ДПФ рядов as и xs :
ys = F -1(Wm) =
N-1
= 1/N ∑ Bm* Vm exp (2 π i m s / N) , s = 0, 1,..., N-1,
m=0
Для эффективного вычисления ДПФ и ОДПФ применяется быстрое преобразование Фурье [2], за счет которого число операций в pеализованном алгоритме пропорционально N log2 N и требуется память длины 2N.
Описанный метод можно также применять для быстрого умножения двух циркулянтных матриц: А и матрицы C вида (1), определяемой вектором x = (x0, x1,..., xN - 1), при этом вектор y = Ax = (y0, y1,..., yN - 1), определяет циркулянтную матрицу вида (1) произведения A C.
| 1. | Б.Голд, Ч.Pейдер. Цифровая обработка сигналов. М.,"Советское радио", 1973. |
| 2. | В.А.Морозов, Н.Н.Кирсанова, А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов, Сб. "Численный анализ на ФОPТPАНе", вып.15, Изд - во МГУ, М., 1976. |
procedure AMC1R(var N :Integer; var B :Array of Real;
var X :Array of Real);
Параметры
| N - | размерность исходной матрицы А, равная целой степени числа два, N ≥ 4 (тип: целый); |
| B - | вещественный одномерный массив длины N, содержащий заданные определяющие элементы матрицы А: a0, a1, ..., aN - 1, |
| X - | вещественный одномерный массив длины N, содержащий заданные компоненты вектора x = (x0, x1,..., xN - 1) на входе и вычисленные компоненты вектора y = (y0, y1,..., yN - 1) на выходе. |
Версии
| AMC1C - | вычисление произведения комплексной циркулянтной матрицы А на комплексный вектор X или произведения двух комплексных циркулянтных матриц. Для подпрограммы АMC1C параметры B, X имеют тип Complex. |
Вызываемые подпрограммы
|
FTF1C - FTF2C | подпрограммы вычисления Дискретного или Обратного Дискретного Преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье (отличаются друг от друга способом представления исходного ряда); используются в АМС1R и АMC1C соответственно. |
Замечания по использованию
| 1. |
После выхода из подпрограммы АМС1R и АMC1C значения определяющих элементов исходной матрицы А в массиве B не сохраняются. | |
| 2. |
Если циркулянтная матрица А имеет вид | a0 a1 a2 ... aN-1 |
| aN-1 a0 a1 ... aN-2 |
| aN-2 aN-1 a0 ... aN-3 |
| . . . . . . . . . . . . . . . . . . . . .
| a1 a2 a3 ... a0 |
и задается N элементами первой строки a0, a1, ..., aN - 1, то до обращения к подпрограммам АМС1R и АMC1C эти элементы нужно упорядочить в массиве B следующим образом: a0, aN - 1, aN - 2, ..., a1. |
1.
Unit tamc1r_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, AMC1R_p;
function tamc1r: String;
implementation
function tamc1r: String;
var
N,I :Integer;
const
B :Array [0..7] of Real = ( 1.0,8.0,7.0,6.0,5.0,4.0,3.0,2.0 );
X :Array [0..7] of Real = ( 1.0,-1.0,-2.0,3.0,0.0,2.0,-3.0,0.0 );
label
_1;
begin
Result := '';
N := 8;
AMC1R(N,B,X);
for I:=1 to N do
begin
_1:
Result := Result + Format(' %20.16f ',[X[I-1]]) + #$0D#$0A;
end;
UtRes('tamc1r',Result); { вывод результатов в файл tamc1r.res }
exit;
end;
end.
Результат: y = (-4., 4., -4., -20., 4., 4., 20., -4.)
2.
Unit tamc1c_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AMC1C_p;
function tamc1c: String;
implementation
function tamc1c: String;
var
N,I :Integer;
const
B :Array [0..3] of Complex = ( ( re:1.0; im:1.0 ),( re:4.0; im:-1.0 ),(
re:3.0; im:2.0 ),( re:2.0; im: 0.0 ) );
X :Array [0..3] 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 ) );
label
_1;
begin
Result := ''; { результат функции }
N := 4;
AMC1C(N,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('tamc1c',Result); { вывод результатов в файл tamc1c.res }
exit;
end;
end.
Результат: y = ( (2., -12.), (11., -4.), (4., -5.), (-3., 3.) )