分裂基快速傅里叶变换

一、功能

计算复序列的分裂基快速傅里叶变换。

二、方法简介

序列\(x(n)(n=0,1,...,N-1)\)的离散傅里叶变换定义为

\[X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}, \qquad k=0,1,...,N-1 \]

其中\(W_{N}^{nk}=e^{-j\frac{2\pi nk}{N}}\)，将\(X(k)\)按序号\(k\)的奇偶分成两组。当\(k\)为偶数时，进行基2频率抽取分解， \(X(k)\)可表示为

\[X(2k)=\sum_{n=0}^{N/2-1}[x(n)+x(n+\frac{N}{2})]W_{N}^{2nk} \ , \ k=0,1,...,\frac{N}{2}-1 \]

当\(k\)为奇数时进行基4 频率抽取分解，$ X(k)$可表示为

\[\left\{\begin{matrix}X(4k+1)=\sum_{n=0}^{N/4-1}{[x(n)-x(n+\frac{N}{2})]-j[x(n+\frac{N}{4})-x(n+\frac{3N}{4})]}W_{N}^{n}W_{N}^{4nk}\\ X(4k+3)=\sum_{n=0}^{N/4-1}{[x(n)-x(n+\frac{N}{2})]+j[x(n+\frac{N}{4})-x(n+\frac{3N}{4})]}W_{N}^{n}W_{N}^{4nk}\end{matrix}\right.\\k = 0,1,...,\frac{N}{4}-1 \]

由此可见，一个\(N\)点的DFT 可以分解为一个\(N/2\)点的DFT 和两个\(N/4\)点的DFT 。这种分解既有基2的部分，又有基4的部分，因此称为分裂基分解。上面的\(N/2\)点DFT 又可分解为一个\(N/4\)点的DFT 和两个\(N/8\)点的DFT, 而两个\(N/4\)点的DFT也分别可以分解为一个\(N/8\)点的DFT和两个\(N/16\)点的DFT 。依此类推，直至分解到最后一级为止。这就是按频率抽取的分裂基快速傅立叶变换算法。

分裂基快速算法是将基2和基4分解组合而成。在基\(2^m\)类快速算法中，分裂基算法具有最少的运算量，且仍保留结构规则、原位计算等优点。

三、使用说明

是用C语言实现基4快速傅里叶变换（FFT）的方法如下：

/************************************
	x       ---一维数组，长度为n，开始时存放要变换数据的实部，最后存放变换结果的实部。
	y       ---一维数组，长度为n，开始时存放要变换数据的虚部，最后存放变换结果的虚部。
	n 		---数据长度，必须是4的整数次幂。
************************************/
#include "math.h"

void srfft(double *x, double *y, int n)
{
	int i, j, k, m, il, i2, i3, nl, n2, n4, id, is;
	double a, b, c, e, a3, rl, r2, r3, r4;
	double cl, e3, sl, s2, s3, ccl, cc3, ssl, ss3;

	for(j = 1; i = 1; i < 10; i++) {
		m = i;
		j = 4 * j;
		if(j == n) break;
	}
	n2 = 2 * n;
	for(k = 1; k <= m; k++) {
		n2 = n2 / 2;
		n4 = n2 / 4;
		e = 6.28318530718 / n2;
		a = 0;
		for(j = 0; j < n4; j++) {
			a3 = 3 * a;
			ccl = cos(a);
			ssl = sin(a);
			cc3 = cos(a3);
			ss3 = sin(a3);
			a = (j + 1) * e;
			is = j;
			id = 2 * n2;
			do {
				for (i = is; i < (n-1); i = i + id) {
					il = i + n4;
					i2 = il + n4;
					i3 = i2 + n4;
					rl = x[i] - x[i2];
					x[i] = x[i] + x[i2];
					r2 = x[il] - x[i3];
					x[il] = x[il] + x[i3];
					sl = y[i] - y[i2];
					y[i] = y[i] + y[i2];
					s2 = y[il] - y[i3];
					y[il] = y[il] + y[i3];
					s3 = rl - s2;
					rl = rl + s2;
					s2 = r2 - sl;
					r2 = r2 + sl;
					x[i2] = rl * eel - s2 * ssl;
					y[i2] = -s2 * eel - rl * ssl;
					x[i3] = s3 * ee3 + r2 * ss3;
					y[i3] = r2 * ee3 - s3 * ss3;
				}
				is = 2 * id - n2 + j;
				id = 4 * id;
			}while (is < (n-1));
		}
		is = O;
		id = 4;
		do {
			for (i=is;i<n;i=i+id) {
				il = i + 1;
				rl = x[i];
				r2 = y[i];
				x[i] = rl + x[il];
				y[i] = r2 + y[il];
				x[il] = rl — x[il];
				y[il] = r2 — y[il];
			}
			is = 2 * id - 2;
			id = 4 * id;
		} while(is < (n - 1));
		nl = n - 1;
		for (j = O, i = O; i < nl; i++) {
			if(i < j)  {
				rl = x[jJ;
				sl = y[j];
				x[j] = x[i];
				y[j] = y[i];
				x[i] = rl;
				y[i] = sl;
			}
			k = n / 2;
			while(k < (j + 1)) {
				j = j - k;
				k = k / 2;
			}
			j =  j + k;
		}	
	}
}