共轭对称序列的快速傅里叶反变换

一、功能

计算共轭对称复序列的快速傅里叶反变换，其变换结果是实数。

二、方法简介

序列\(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;
}