快速傅里叶变化（FFT）是离散傅里叶变化（DFT）的快速算法。一个长度为N的有限长序列的DFT为
式中的 为 。离散傅里叶反变换为
FFT算法的基本思想是利用 的对称性和周期性，即 =- 和 （r为任意整数），通过将长序列分解为短序列，再由短序列的DFT组合来实现整个长序列的DFT，从而减少DFT的复数乘法和复数加法次数，提高DFT的运算速度和效率。将N点的长序列分解为短序列可以有多种方式，常用的有基2时间抽取FFT算法，基2频率抽取FFT算法，基4抽取FFT算法和分裂基FFT算法。基2抽取算法简单有效，基4抽取算法和分裂基抽取算法与基2抽取算法原理类似，但基4算法和分裂基算法效率更高，运算量更少。本文主要介绍基2抽取FFT算法。
1、基2时间抽取FFT算法：
基2时间抽取FFT算法是将N点序列x[n]按序号奇偶分解得到两个 点子序列 和 ,通过计算这两个子序列的DFT 和 实现 的DFT 。算法要点为N点序列x[n]对时域下标按奇数和偶数分解，即
（1）
对频域X[k]按前一半和后一半分解，即 。频域的前一半可由下式求得
由DFT的周期性和 的对称性有
-=
故频域的后一半可由下式求得
将前一半和后一半组合起来，即得原来的长序列x[n]的DFT 。
可对 和 继续按式(1)分解下去，直到每个子序列的长度为1为止。任何一个N= 点的DFT，可由通过M次分解，得到长度为1的子序列。对于N不是2的整数次幂的情况，可在序列后面添0补齐到2的整数次幂。每级分解运算可用一个蝶形运算流程图来表示，如图
下图是对一个长度为8的序列进行基2时间抽取FFT运算的流程图
直接计算N点序列的DFT需要N2次复数乘法，需要N(N-1)次复数加法。而用基2时间抽取FFT算法，则乘法次数为 ，加法次数 ，运算次数大大减少。
2、基2频率抽取FFT算法
基2频率抽取FFT算法则是对频域下标按奇数和偶数分解，对时域下标按前一半和后一半分解。 的DFT 可表示为
将 分解成 和 两个长度为 的子序列，如下：
（2）
则 的DFT 可按频率下标分奇数和偶数表示为
可继续将 和 按式(2)分解下去，直到每个子序列长度为1。任何一个N= 点的DFT，可由通过M次分解，得到长度为1的子序列。对于N不是2的整数次幂的情况，可在序列后面添0补齐到2的整数次幂。每级分解运算也可用一个蝶形运算流程图来表示，如图
下图是对一个长度为8的序列进行基2频率抽取FFT运算的流程图：
3、编程实现
Matlab有一个内置的函数fft()可以实现快速傅里叶变换，由于它是built-in函数，所以无法查看源码来判断它是应用的那种分解方式来实现的快速傅里叶变换。下面将用Matlab语言来编写实现类似fft()的快速傅里叶变换函数myfft()。本函数采用的FFT算法为基2时间抽取FFT算法，函数定义源码如下：
%----------------------------------
%myfft(A,M)有两个形参，A为要进行快速傅里叶变换的序列，M为序列长度对2的对数，
%即如果序列A[n]的长度为8，则M=3
function [A] = myfft(A,M)
N=2^M; %序列的长度为N
LH=N/2; %序列长度的一半
N1=N-2;
J=LH
for I=1:1:N1 %将输入序列A[n]按运算流程图第一列的顺序排好
if I<J
temp=A(I+1);
A(I+1)=A(J+1);
A(J+1)=temp;
end
K=LH;
while J>=K
J=J-K;
K=K/2;
end
J=J+K;
end
for L=1:1:M %L表示序列分解层次，此次分解有M-L个子序列
B=2^(L-1);
for S=0:B-1 %执行第S个蝶形运算，每个子序列共有B个蝶形运算
p=S*2^(M-L); %旋转因子的上标，即此时旋转因子为
for k=S:2^L:N-1 %执行各个子序列中的第S个蝶形运算，有M-L个子序列
temp=A(k+1)+A(k+B+1)*exp(-i*2*pi*p/N);
A(k+B+1)=A(k+1)-A(k+B+1)*exp(-i*2*pi*p/N);
A(k+1)=temp; %蝶形运算
end
end
end
end
python语法简洁清晰，具有强大和丰富的三方库，能把其他语言(如C/C++,Matlab，java等)制作的各种模块轻松联结在一起，开发效率高，很适合快速开发实验原型以验证算法正确性和有效性。下面是基2时间抽取FFT算法的一个python版本实现。
#-------------------------------------
import numpy as np
def FFT(x):
x =np.asarray(x, dtype=float)
N =x.shape[0]
ifN%2>0: #输入序列的长度必须为2的幂
raiseValueError("size of x must be a power of 2")
elif N<= 2: # 子序列长度小于等于2时，直接计算
x =np.asarray(x, dtype=float)
N = x.shape[0]
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
eturn np.dot(M, x)
else:
X_even= FFT(x[::2]) #偶序列继续分解
X_odd= FFT(x[1::2]) #奇序列继续分解
factor= np.exp(-2j * np.pi * np.arange(N/2) / N) #旋转因子
returnnp.concatenate([X_even + factor* X_odd, #返回序列频域前半部分
X_even - factor*X_odd]) #返回序列频域后半部分
4、函数测试和实验
快速傅里叶变换可以得到信号的频谱。通过对一个给定的信号进行快速傅里叶变换，观察它的频谱，并和Matlab内置函数fft()变换得到的结果进行对比，来对myfft()的功能进行测试。测试可用如下脚本来实现：
%----------------------------
%测试脚本，分别用myfft()和fft()对一个信号进行快速傅里叶变换，该信号由一个
%50Hz和120Hz的正弦信号叠加，频谱会在这50Hz和120Hz出现峰值。
Fs = 1000; % 采样频率
T = 1/Fs; % 采样周期
L = 5000; % 信号长度
t = (0:L-1)*T; % 时间序号数组
% x(t)=0.7sin(2 *50t)+sin(2 *150t),为50Hz和150Hz正弦信号的叠加
x = 0.7*sin(2*pi*50*t) + sin(2*pi*150*t);
y = x + 2*randn(size(t)); % 给信号加上噪声
%画出x(t)和y(t)的图像
subplot(2,1,1)
plot(Fs*t(1:50),x(1:50))
title('Sinusoids Signal')
xlabel('time (milliseconds)')
subplot(2,1,2)
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
%对加了噪声的y(t)进行快速傅里叶变换
NFFT = 2^nextpow2(L); %延长信号长度为最小的2的幂
y=[y zeros(1,2^nextpow2(L)-L)]; %延长信号长度，后面用0补齐
Y = myfft(y,log2(NFFT))/L; % 用自己实现的快速傅里叶变化myfft处理信号得到Y（f）
Y2 = fft(y,NFFT)/L; %用Matlab官方实现的fft处理信号得到Y2(f）
f = Fs/2*linspace(0,1,NFFT/2+1);%频域横坐标
%对信号y(t)进行快速傅里叶变换，并画出的频谱
figure
subplot(2,1,1)
plot(f,2*abs(Y(1:NFFT/2+1))) %
title('myfft()处理的y(t)频谱')% y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,1,2)
plot(f,2*abs(Y2(1:NFFT/2+1))) %
title('fft()处理的y(t)频谱')
xlabel('Frequency (Hz)')
ylabel('|Y2(f)|')
%对原信号x(t)进行快速傅里叶变换，并画出频谱
x=[x zeros(1,2^nextpow2(L)-L)]; %延长信号长度，后面用0补齐
X = myfft(x,log2(NFFT))/L; % 用自己实现的快速傅里叶变化myfft处理信号得到X（f）
X2 = fft(x,NFFT)/L; %用Matlab官方实现的fft处理信号得到X2(f）
figure
subplot(2,1,1)
plot(f,2*abs(X(1:NFFT/2+1))) %
title('myfft()处理的x(t)频谱')%
xlabel('Frequency (Hz)')
ylabel('|X(f)|')
subplot(2,1,2)
plot(f,2*abs(X2(1:NFFT/2+1))) %
title('fft()处理的x(t)频谱')
xlabel('Frequency (Hz)')
ylabel('|X2(f)|')
运行脚本，得到x(t)和y(t)的频谱，如图。可见，myfft()和Matlab内置函数fft()的处理结果一致，没有区别。频谱图在50Hz和150Hz出现峰值，幅值分别为0.63和1.0左右。
Python版的测试脚本如下
a=np.array([1,0,1,0,1,0,1,0]) #定义序列a[n]=[1,0,1,0,1,0,1,0]
b=[1,1,1,1,0,0,0,0] #定义序列b[n]=[1,1,1,1,0,0,0,0]
A=FFT(b) #对a[n]进行快速傅里叶变换
B=FFT(a) #对b[n]进行快速傅里叶变换
printA,”\n”,B #输出两个序列的傅里叶变换的结果A[M]和B[m]
python版的傅里叶变换结果和Matlab自带的fft()变换结果对比如下
Python版的运行结果
Matlab版的运行结果
注意python版本的运行结果中，有”e-16”的结果，这是因为计算机浮点数的运算误差导致的，实际上可认为是0。可见，python实现的快速傅里叶变换函数的结果和Matlab的fft()变换结果一致，python实现的快速傅里叶变换结果正确无误。
5、快速傅里叶变换的应用
快速傅里叶变换可以将有限长或者周期离散时间信号从时域表示转换到频域表示。快速傅里叶变换除了可以得到一个信号的频谱外，还可以基于得到的频谱对信号进行一些频域上的处理，并通过快速傅里叶反变换得到经过处理后的信号。比如计算离散时间信号的线性卷积，滤波，信号压缩，降噪，谱估计等。FFT除了用于离散时间信号处理，也可以用于逼近连续时间信号的频谱。
下面脚本展示了利用快速傅里叶变换和反变换进行快速卷积计算
subplot(3,1,1);
n=0:16;
x=0.8.^n; %x[n]=0.8^n的离散信号
stem(n,x);
xlabel('n');ylabel('x[n]');
subplot(3,1,2);
n=0:15;
y=[ones(1,10) zeros(1,6)];
stem(n,y);
xlabel('n');ylabel('y[n]')
subplot(3,1,3);
L=26;n=0:L-1;
X=fft(x,L);
Y=fft(y,L);
Z=X.*Y; %在频域两信号的频谱相乘
z=ifft(Z,L); %快速傅里叶反变换得到两信号的卷积
stem(n,z);
xlabel('n');
ylabel('z[n]')
脚本运行结果如图
在工程领域，FFT可以用于语音处理，工程传感器信号处理，特征提取等。
二维FFT和IFFT可应用到图像处理领域，比如图像降噪，模糊化，清晰化，轮廓化，压缩，还原等。快速傅里叶变换及其反变换是图像处理领域的重要处理手段，是进行其他各项处理的基础，应用场合很广。下面脚本展示了图像的二维快速傅里叶变换及其反变换，及基于快速傅里叶变换得到的频谱图，相位图等一系列频域信息。
C=imread('2.jpg');%载入图片
A=rgb2gray(C);
B=fftshift(fft2(A));% 进行傅立叶变换
E=fft2(A);
D=ifft2(E);
subplot(331)
imshow(C);
title('原始图像');
subplot(332)
imshow(A);
title('原始灰度图像');
subplot(339)
imshow(abs(B),[]);
title('原始频谱图');
subplot(334)
imshow(log(abs(B)),[]);
title('取对数后的频谱图');
subplot(335)
imshow(angle(B),[]);
title('相位图');
subplot(336)
imshow(real(B),[]);
title('实部图');
subplot(337)
imshow(imag(B),[]);
title('虚部图');
subplot(338)
imshow(abs(E),[]);
title('fft2频谱图');
subplot(333)
imshow(D,[]);
title('经快速傅里叶逆变换得到的灰度图');
下面为脚本运行的结果，可见经快速傅里叶反变换后，可以从图像的频域表示得到时域表示，即重新得到一张图像。由于此脚本没有对图像做其他处理，所以快速傅里叶反变换得到的图形与原来的灰度图像一致。
6、总结
快速傅里叶变换大大简化了运算量，并且算法具有迭代和分治的特点，特别适合计算机实现，使得用计算机对离散时间信号进行快速处理成为现实。配合采样定理和信号采样及重构等方法，使得快速傅里叶变换可以进一步推广到连续时间信号处理，大大提高了信号处理的效率和速度。在语音处理领域，傅里叶变换可以用于语音识别，音频分离，降噪等应用。在图像应用领域，快速傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。而灰度分布函数即为灰度图的表达。快速傅里叶变换使得对图像可以用频域处理的方法来处理，对基于计算机的图像处理应用的可行性和实时性具有重大意义，是数字信号处理的基石。
743
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
