Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://num-anal.srcc.msu.ru/lib_na/cat/r_htm_c/rsc2c_c.htm
Дата изменения: Mon Mar 2 17:08:12 2015 Дата индексирования: Sun Apr 10 01:46:00 2016 Кодировка: Windows-1251 |
Текст подпрограммы и версий rsc2c_c.zip |
Тексты тестовых примеров trsc2c_c.zip |
Вычисление двумерных периодических свеpток двух комплексных матриц с применением алгоритма быстрого преобразования Фурье.
Заданы две комплексные матрицы A, B одинаковой размерности N1 * N2 с элементами
An1, n2 , n1 = 0, 1, ..., N1 - 1 , n2 = 0, 1, ..., N2 - 1 , Br1, r2 , r1 = 0, 1, ..., N1 - 1 , r2 = 0, 1, ..., N2 - 1 .
Требуется вычислить периодическую (круговую) свеpтку C этих матриц, т.е.
N1-1 N2-1 (1) ∑ ∑ As1 - r1, s2 - r2 Br1, r2 = Cs1, s2 r1=0 r2=0 s1 = 0, 1, ..., N1 - 1 , s2 = 0, 1, ..., N2 - 1 ,
где элементы матрицы A продолжены периодически для отрицательных значений индексов с периодом N1 в каждом столбце и с периодом N2 в каждой стpоке:
A- n1, n2 = AN1 - n1, n2 , n1 = 1, 2, ..., N1 - 1 , n2 = 0, 1, ..., N2 - 1 , An1, - n2 = An1, N2 - n2 , n1 = 0, 1, ..., N1 - 1 , n2 = 1, 2, ..., N2 - 1 .
Пусть F - оператор двумерного Дискретного Преобразования Фурье, задаваемый формулой (например, для матрицы A):
Am1, m2 = F(An1, n2) =
N1-1 N2-1
= ∑ ∑ An1, n2 * exp [ - 2π i ( n1*m1/N1 + n2*m2/N2)] ,
n1=0 n2=0
m1 = 0, 1, ..., N1 - 1 , m2 = 0, 1, ..., N2 - 1 , i = √-1 .
Аналогично Bm1, m2 = F (Br1, r2), Cm1, m2 = F (Cs1, s2).
По теоpеме о дискретной периодической свеpтке [1]
(2) Am1, m2 * Bm1, m2 = Cm1, m2 , m1 = 0, 1, ..., N1 - 1 , m2 = 0, 1, ..., N2 - 1 .
Искомая матрица С находится путем обратного Дискретного Преобразования Фурье F - 1 от произведения ДПФ матриц A и B:
Cs1, s2 = F - 1 ( Cm1, m2 ) = N1-1 N2-1 = 1/(N1*N2) ∑ ∑ Am1,m2* Bm1,m2* exp[2π i (s1*m1/N1 + s2*m2/N2)], m1=0 m2=0 s1 = 0, 1, ..., N1 - 1 , s2 = 0, 1, ..., N2 - 1 .
Подпрограмма предусматривает работу в двух режимах (в зависимости от значения параметра L).
1. | L=1. | Производится вычисление свертки (1) с вычислением ДПФ каждой из двух матриц A, B. |
2. | L=2. | Производится вычисление свертки (1) в предположении, что ДПФ матрицы A уже вычислено. Этот режим позволяет избежать повторного получения комплексной матрицы A при многократном вычислении свертки (1) с той же матрицей A. |
Для эффективного вычисления ДПФ и ОДПФ применяется двумерное Быстрое Преобразование Фурье [2]. При этом время вычислений пропорционально N1 N2 (log2 N1 + log2 N2), а объем памяти пропорционален N1 * N2.
Аналогичный алгоритм для вычисления одномерных периодических сверток двух последовательностей описан в подпрограмме amc1r_c.
1. | Л.Рабинер, Б.Гоулд. Теория и применение цифровой обработки сигналов. M., "Мир", 1978. |
2. | В.А.Морозов, Н.Н.Кирсанова, А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов, сб. "Численный анализ на ФОРТРАНЕ", вып.15, Изд - во МГУ, M., 1976. |
int rsc2c_c (real *ar, real *ai, real *br, real *bi, integer *n1, integer *n2, integer *l)
Параметры
ar,ai - | двумерные вещественные массивы размера n1 * n2, содержащие соответственно вещественные и мнимые части элементов заданной матрицы A; |
br, bi - | двумерные вещественные массивы размера n1 * n2, содержащие соответственно вещественные и мнимые части элементов заданной матрицы B; на выходе подпрограммы содержат соответственно вещественные и мнимые части элементов вычисленной матрицы C; |
n1 - | заданное число стpок матриц A и B, n1 ≥ 4 - целая степень числа два (тип: целый); |
n2 - | заданное число столбцов матриц A и B, n2 ≥ 4 - целая степень числа два (тип: целый); |
l - | параметр, определяющий режим использования подпрограммы (тип: целый); |
l=1 - | вычисление свертки C матриц A, B в общем случае, |
l=2 - | вычисление свертки C матриц A, B в предположении, что вещественные и мнимые части элементов ДПФ матрицы A уже содержатся соответственно в массивах ar, ai. |
Версии: нет
Вызываемые подпрограммы
ftftc_c - | вычисление двумерного дискретного или обратного дискретного преобразования Фурье комплексной матрицы размера n1 * n2, n1 = 2 j1, n2 = 2 j2 методом быстрого преобразования Фурье. |
Замечания по использованию
1. |
B результате работы подпрограммы rsc2c_c при l = 1 в массивах ar, ai содержатся соответственно вещественные и мнимые части элементов ДПФ матрицы A. | |
2. |
При первом обращении к подпрограмме (l = 1) необходимо задать значения параметров ar, ai, br, bi, n1, n2. При повторном вычислении свертки (1) с той же матрицей A (l = 2) можно менять только параметры br, bi. | |
3. |
Kак следует из формулы (2), сомножители A, B в свеpтке (1) перестановочны, поэтому при необходимости повторного вычисления свертки с той же матрицей B достаточно при l = 1 матрицу B задать в массивах ar, ai, а матрицу A - в массивах br, bi и далее воспользоваться режимом при l = 2. |
int main(void) { /* Local variables */ extern int rsc2c_c(float *, float *, float *, float *, int *, int *, int *); static int i__, j, l, n1, n2; static float ai[32] /* was [4][8] */, bi[32] /* was [4][8] */, ar[32] /* was [4][8] */, br[32] /* was [4][8] */; int i__1, i__2; #define ai_ref(a_1,a_2) ai[(a_2)*4 + a_1 - 5] #define bi_ref(a_1,a_2) bi[(a_2)*4 + a_1 - 5] #define ar_ref(a_1,a_2) ar[(a_2)*4 + a_1 - 5] #define br_ref(a_1,a_2) br[(a_2)*4 + a_1 - 5] n1 = 4; n2 = 8; i__1 = n1; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = n2; for (j = 1; j <= i__2; ++j) { ar_ref(i__, j) = (2.f - i__) * (j - 4.f); ai_ref(i__, j) = 0.f; br_ref(i__, j) = 6.f - i__ - j; /* l1: */ bi_ref(i__, j) = (float) (i__ - j); } } l = 1; rsc2c_c(ar, ai, br, bi, &n1, &n2, &l); for (i__ = 1; i__ <= 4; ++i__) { printf("\n %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f \n", br_ref(i__, 1), br_ref(i__, 2), br_ref(i__, 3), br_ref(i__, 4), br_ref(i__, 5), br_ref(i__, 6), br_ref(i__, 7), br_ref(i__, 8)); } for (i__ = 1; i__ <= 4; ++i__) { printf("\n %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f \n", bi_ref(i__, 1), bi_ref(i__, 2), bi_ref(i__, 3), bi_ref(i__, 4), bi_ref(i__, 5), bi_ref(i__, 6), bi_ref(i__, 7), bi_ref(i__, 8)); } return 0; } /* main */ Результаты: | - 16. 24. 48. 56. 48. 24. - 16. - 72. | br = | - 8. 32. 56. 64. 56. 32. - 8. - 64. | | - 16. 24. 48. 56. 48. 24. - 16. - 72. | | - 40. 0. 24. 32. 24. 0. - 40. - 96. | | - 16. 24. 48. 56. 48. 24. - 16. - 72. | bi = | - 24. 16. 40. 48. 40. 16. - 24. - 80. | | - 16. 24. 48. 56. 48. 24. - 16. - 72. | | 8. 48. 72. 80. 72. 48. 8. - 48. |