一、功能
計算共軛對稱復序列的快速傅里葉反變換，其變換結果是實數。
二、方法簡介
序列\(x(n)\)的離散傅里葉變換為
\[X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}, \ k=0,1,...,N-1 \]
序列\(X(k)\)的離散傅里葉反變換為
\[x(n)=\frac{1}{N}\sum_{n=0}^{N-1}X(k)W_{N}^{-nk}, \ n=0,1,...,N-1 \]
共軛對稱復序列的傅里葉反變換，可用復序列快速傅里葉反變換算法進行計算。但考慮到\(X(k)\)是共軛對稱的，其傅里葉反變換\(x(n)\)是實數，因此，為進一步提高計算效率，需要對一般的復序列IFFT算法進行一定的修改。
共軛對稱序列\(X(k)\)具有如下性質：\(X(0)\)和\(X(\frac{N}{2})\)都是實數，且有
\[X(k)=X^{*}(N-k), \ 1\leqslant k\leqslant \frac{N}{2}-1 \]
即\(X(k)\)的實部是偶對稱，虛部是奇對稱。在計算傅里葉反變換時，利用這種共軛對稱性，我們就可以不必計算和存儲\(X(k)(\frac{N}{2}+1\leqslant k\leqslant N-1)\)以及\(X(0)\)和\(X(\frac{N}{2})\)的虛部，這比一般形式的快速傅里葉反變換算法大約可減少一半的運算量和存儲量。具體計算時采用的是分裂基算法。
三、使用說明
C語言實現方式如下
/************************************
x ----長度為n。開始時存放要具有共軛對稱性的復序列X(k)的前n/2+1個值，
其存儲順序為[Re(0),Re(1),...,Re(n/2),Im(n/2-1),...,Im(1)],
其中Re(0)=X(0),Re(n/2)=X(n/2)。最后存放變換結果x(i)(i=0,0,1,...,n-1)，
這里x(i)是實數。
n ----數據長度，必須是2的整數次冪，即n=2^m。
************************************/
#include "math.h"
void irfft(double *x, int n)
{
int i, j, k, m, i1, i2, i3, i4, i5, i6, i7, i8, n2, n4, n8, id, is;
double a, e, a3, t1, t2, t3, t4, t5, cc1, cc3, ss1, ss3;
for (j = 1, i = 1; i < 16; i++) {
m = i;
j = 2 * j;
if(j == n)
break;
}
n2 = 2 * n;
for(k = 1; k < m; k++){
is = 0;
id = n2;
n2 = n2 / 2;
n4 = n2 / 4;
n8 = n4 / 2;
e = 6. 28318530718 / n2;
do {
for(i = is; i < n; i += id){
i1 = i;
i2 = i1 + n4;
i3 = i2 + n4;
i4 = i3 + n4;
t1 = x[i1] - x[i3];
x[i1] = x[i1] + x[i3];
x[i2] = 2 * x[i2];
x[i3] = t1 - 2 * x[i4];
x[i4] = t1 + 2 * x[i4];
if(n4 == 1)
continue;
i1 += n8;
i2 += n8;
i3 += n8;
i4 += n8;
t1 = (x[i2] - x[i1]) / sqrt(2.0);
t2 = (x[i4] + x[i3]) / sqrt(2.0);
x[i1] = x[i1] + x[i2];
x[i2] = x[i4] - x[i3];
x[i3] = 2 * (-t2 - t1);
x[i4] = 2 * (-t2 + t1);
}
is = 2 * id - n2;
id = 4 * id;
}while(is < (n - 1));
a = e;
for(j = 1; j < n8; j++) {
a3 = 3 * a;
cc1 = cos(a);
ss1 = sin(a);
cc3 = cos(a3);
ss3 = sin(a3);
a = (j + 1) * e;
is = 0;
id = 2 * n2;
do {
for(i = is; i <= (n - 1); i = i + id) {
i1 = i - j;
i2 = i1 + n4;
i3 = i2 + n4;
i4 = i3 + n4;
i5 = i + n4 - j;
i6 = i5 + n4;
i7 = i6 + n4;
i8 = i7 + n4;
t1 = x[i1] - x[i6];
x[i1] = x[i1] + x[i6];
t2 = x[i5] - x[i2];
x[i5] = x[i2] + x[i5];
t3 = x[i8] + x[i3];
x[i6] = x[i8] - x[i3];
t4 = x[i4] + x[i7];
x[i2] = x[i4] - [i7];
t5 = t1 - t4;
t1 = t1 + t4;
t4 = t2 - t3;
t2 = t2 + t3;
x[i3] = t5 * cc1 + t4 * ss1;
x[i7] = -t4 * cc1 + t5 * ss1;
x[i4] = t1 * cc3 - t2 * ss3;
x[i8] = t2 * cc3 + t1 * ss3;
}
is = 2 * id - n2;
id = 4 * id;
}
while(is < (n - 1));
}
}
is = 0;
id = 0;
do {
for((i = is; i < n; i = i + id)){
i1 = i + 1;
t1 = x[i];
x[i] = t1 + x[i1];
x[i1] = t1 - x[i1];
}
is = 2 * id - 2;
id = 4 * id;
}while(is < (n - 1));
for(j = 0, i = 0; i < (n - 1); i++){
if(i < j){
t1 = x[j];
x[j] = x[i];
x[i] = t1;
}
k = n / 2;
while(k < (j + 1)){
j = j - k;
k = k / 2;
}
j = j + k;
}
for(i = 0; i < n; i++)
x[i] = x[j] / n;
}