|
Текст подпрограммы и версий rsa1r_c.zip , rsa1d_c.zip |
Тексты тестовых примеров trsa1r_c.zip , trsa1d_c.zip |
Вычисление прямой и обратной дискретной периодической (круговой) свертки двух одномерных массивов вещественных чисел.
Пусть известны N значений sj = s ( tj ) функции сигнала s ( t ) на равномерной сетке tj = j Δt, j = 0, 1, 2, ..., N - 1 (где Δt - шаг сетки, а N равняется целой степени двух). Предполагается, что сигнал sj является периодическим с периодом N (отсюда, в частности, следует, что значения sj известны и для отрицательных значений j, т.е. s- j = sN - j, j = 1, 2, ..., N - 1). Значения sj задаются в вещественном одномерном массиве ARRAY длины N: ARRAY ( J ) = sj ( J = 1, 2, ..., N; j = 0, 1, 2, ..., N - 1).
Пусть известны также M значений rk = r ( tk ) импульсной функции отклика r ( t ) на той же равномерной сетке tk = k Δt, k = - (M - 1)/2, ..., 0, ..., (M - 1)/2 (где Δt - тот же шаг сетки, что и для функции сигнала, а M - нечетное натуральное число такое, что M < N). Значения rk задаются в M первых элементах одномерного вещественного массива RESP длины N следующим образом:
RESP(1) = r0
RESP( I ) = rk , I = 2, 3, ..., (M+1)/2 ; k = 1, 2, ..., (M-1)/2
RESP( I ) = rk , I = (M+3)/2, ..., M ; k = (M-1)/2, ..., - 1
Подпрограмма rsa1r_c имеет два режима работы, задаваемых при обращении к ней значением параметра IREG. В первом режиме ( IREG = 1) выполняется прямая дискретная периодическая свертка массивов ARRAY и RESP (L = 1, 2, ..., N):
N
(r*s)L = ∑ ARRAY(L - J) × RESP( J )
j =1
где * - знак свертки, а
× - знак умножения.
Предварительно в массив RESP добавляются нули следующим образом:
RESP( I ) = RESP( I ) , I = 1, 2, ..., (M+1)/2 RESP( I ) = 0.0 , I = (M+3)/2, ..., N - (M-1)/2 RESP(N+1-I ) = RESP(M+1-I ) , I = 1, 2, ..., (M-1)/2
Для вычисления указанной свертки используется быстрое дискретное преобразование Фурье в соответствии с теоремой о дискретной периодической (круговой) свертке. Полученные значения (r * s)L помещаются в первых N компонентах массива ANS.
Во втором режиме ( IREG = - 1) выполняется восстановление значений sj функции сигнала по известным значениям rk импульсной функции ответа, задаваемым в массиве RESP, и значениям соответствующей свертки, задаваемым в массиве ARRAY. Полученные значения sj помещаются в первых N элементах массива ANS.
Р.Отнес, Л.Эноксон. Прикладной анализ временных рядов. Изд - во "Мир", 1982.
int rsa1r_c (real *array, integer *n, real *resp, integer *m,
complex *ans, complex *fft, integer *ireg, integer *iflag)
Параметры
| array - | вещественный одномерный массив длины n, содержащий либо значения sj ( ireg = 1), либо значения свертки ( ireg = - 1); |
| n - | длина массива array, равная целой степени двух (тип: целый); |
| resp - | вещественный одномерный массив длины n, в первых m компонентах которого задаются значения rk; |
| m - | количество заданных значений rk (тип: целый); |
| ans - | вещественный одномерный массив длины 2n, в первых n элементах которого в результате работы подпрограммы размещаются либо компоненты вектора свертки (ireg = 1), либо восстановленные значения sj (ireg = - 1); |
| fft - | комплексный одномерный массив длины n, используемый в подпрограмме в качестве рабочего; |
| ireg - | задает режим работы подпрограммы (тип: целый); при этом: |
| ireg= 1 - | когда вычисляется свертка; |
| ireg= -1 - | когда восстанавливается функция сигнала sj |
| iflag - | целая переменная, служащая для сообщения о характере окончания работы подпрограммы во втором режиме ( ireg = - 1); при этом: |
| iflag=0 - | когда восстановлены все значения sj; |
| iflag=i - | когда i - ая компонента преобразования Фурье импульсной функции отклика rk равна нулю; это означает, что значение si не может быть восстановлено из - за потери информации в векторе свертки; значения sj в этом случае подпрограммой не восстанавливаются |
Версии
| rsa1d_c - | вычисление прямой и обратной дискретной периодической (круговой) свертки двух одномерных массивов вещественных чисел в режиме удвоенной точности; при этом параметры array, resp и ans должны иметь тип double, параметр fft - тип doublecomplex (это означает, что данная версия существует только для тех компиляторов, входной язык которых допускает комплексные данные удвоенной точности) |
Вызываемые подпрограммы
|
fta2r_c - fta2d_c | одновременное выполнение прямого быстрого дискретного преобразования Фурье двух одномерных массивов вещественных чисел, каждый из которых имеет длину n, равной целой степени двух, в режимах одинарной и удвоенной точности; используются в подпрограммах rsa1r_c и rsa1d_c соответственно |
|
fta3r_c - fta3d_c | выполнение прямого и обратного быстрых преобразований Фурье одномерного массива вещественных чисел длины 2n, где n равняется целой степени двух, в режимах одинарной и удвоенной точности; используются в подпрограммах rsa1r_c и rsa1d_c соответственно |
Замечания по использованию
| В подпрограммах rsa1r_c и rsa1d_c проверка того, что значение n должно быть целой степенью двух, а m должно быть нечетным числом, меньшим n, не производится |
int main(void)
{
/* Initialized data */
static float array[16] = { 0.f,0.f,.5f,0.f,0.f,-.5f,0.f,0.f,.5f,.5f,.5f,
.5f,.5f,0.f,0.f,0.f };
/* Local variables */
static float resp[16];
extern int rsa1r_c(float *, int *, float *, int *, float *, complex *,
int *, int *);
static int m, n, iflag, i;
static complex fft[16];
static float ans[32];
n = 16;
m = 3;
resp[0] = .7f;
resp[1] = .4f;
resp[2] = .4f;
rsa1r_c(array, &n, resp, &m, ans, fft, &c__1, &iflag);
printf("\n %5i \n", iflag);
for (i = 0; i <= 28; i += 4) {
printf("\n %16.7e %16.7e %16.7e %16.7e \n",
ans[i], ans[i+1], ans[i+2], ans[i+3]);
}
return 0;
} /* main */
Результаты:
ans = (- 0.596046e - 7, 0.2, 0.35, 0.2, - 0.2, - 0.35, - 0.2,
0.2, 0.55, 0.75, 0.75, 0.75, 0.55, 0.2, - 0.745058e - 8,
0.521541e - 07)