一、功能
計算實序列的快速傅里葉變換。
二、方法簡介
實序列\(x(n)\)的離散傅立葉變換為
\[X(k)=\sum_{n=0}^{N-1}x(n)W_{N}^{nk} \ , \ k=0,1,...,N-1 \]
上式可用復序列FFT算法進行計算。但考慮到\(x(n)\)是實數，為進一步提高計算效率，需要對按時間抽取的基2復序列FFT算法進行一定的修改。
如果序列\(x(n)\)是實數，那么其傅立葉變換\(X(k)\)一般是復數，但其實部是偶對稱，虛部是奇對稱，即\(X(k)\)具有如下共輒對稱性： \(X(0)\)和\(X(N/2)\)都是實數，且有
\[X(k)=X^{*}(N-k) \ , \ 1 \leqslant k \leqslant \frac{N}{2} - 1 \]
在計算離散傅立葉變換時，利用這種共輒對稱性，我們就可以不必計算與存儲\(X(k)(N/2 + 1 \leqslant k \leqslant N — 1)\)以及\(X(0)\)和\(X(N/2)\)的虛部，而僅需計算\(X(0)\)到\(X(N/2)\)即可。此處我們選擇的是計算\(X(0)\)到\(X(N/4)\)和\(X(N/2)\)到\(X(3N/4)\), 這樣做可以恰好利用復序列FFT 算法的前\((N/4)+1\)個復數蝶形。這就是按時間抽取的基2實序列FFT算法，它比復序列FFT算法大約可減少一半的運算量和存儲量。
三、使用方法
是用C語言實現實序列快速傅里葉變換的方法如下：
/************************************
x ----長度為n。開始時存放要變換的實數據，最后存放變換結果的前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(k)的共軛對稱性，很容易寫
出后半部分的值。
n ----數據長度，必須是2的整數次冪，即n=2^m。
************************************/
#include "math.h"
void rfft(double *x, int n)
{
int i, j, k, m, il, i2, i3, i4, nl, n2, n4;
double a, e, cc, ss, xt, tl, t2;
for(j = 1, i = 1; i < 16; i++) {
m = i;
j = 2 * j;
if(j == n) break;
}
n1 = n - 1;
for(j = 0, i = 0; i < n1; i++) {
if(i < j) {
xt = x[j];
x[j] = x[i];
x[i] = xt;
}
k = n / 2;
while(k < (j + 1)) {
j = j - k;
k = k / 2;
}
j = j + k;
}
for(i = 0; i < n; i += 2) {
xt = x[i];
x[i] = xt + x[i + 1];
x[i + 1] = xt - x[i + 1];
}
n2 = 1;
for(k = 2; k <= m; k++) {
n4 = n2;
n2 = 2 * n4;
n1 = 2 * n2;
e = 6.28318530718 / nl;
for(i = 0; i < n; i += n1) {
xt = x[i];
x[i] = xt + x[i + n2];
x[i + n2] = xt - x[i + n2];
x[i + n2 + n4] = -x[i + n2 + n4];
a = e;
for(j = 1; j <= (n4-1); j++) {
i1 = i + j;
i2 = i - j + n2;
i3 = i + j + n2;
i4 = i - j + n1;
cc = cos(a);
ss = sin(a);
a = a + e;
t1 = cc * x[i3] + ss * x[i4];
t2 = ss * x[i3] - cc * x[i4];
x[i4] = x[i2] - t2;
x[i3] = -x[i2] - t2;
x[i2] = x[i1] - t1;
x[i1] = x[i1] + t1;
}
}
}
}