一、功能
計算復序列的分裂基快速傅里葉變換。
二、方法簡介
序列\(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;
}
}
}