matlab实现快速傅里叶变换

傅里叶变换是能将信号从时域转到频域的一种变换形式,是我们分析信号的常用方法,毕竟相比于时域的纷乱复杂,频域显得清爽很多。有关傅里叶变换的理解,可以看我之前转载的一篇文章:【转载】傅里叶分析之掐死教程(完整版)
本篇文章主要介绍快速傅里叶变换(FFT)的一种,基2FFT,讲诉计算原理和提供相应的matlab代码。

DFT和FFT

DFT是离散傅里叶变换,是傅里叶变换在离散系统的表现形式。
离散信号x(n)的傅里叶变换可以表示为:

其中k=0,1,2…N-1,所以DFT的运算量为N的平方,当N取值很大时,计算量就会十分巨大。所以需要找一些算法去优化它的计算,减少计算量,这便是FFT,快速傅里叶变换。

一般来说,FFT算法主要分为两种类型,时域抽取和频域抽取。两者的区别主要是提取数列的方法不一致,在计算上表现为蝶形因子的位置不同,前者蝶形因子出现在输入端,后者出现在输出端。

以N=8为例,时域抽取:
时域抽取
频域抽取
频域抽取

基2FFT计算原理

本文所写的FFT以时域抽取为例,频域抽取不作讲解,下面是计算过程。

时间抽取FFT是将N点输入序列x(n),按照奇偶分为奇数列和偶数列。
奇数列:x(1),x(3),x(5)….x(N-1)
偶数列:x(0),x(2),x(4)….x(N-2)

由于旋转因子具备对称性和周期性(周期为N/2),所以可以得到:


如果对X1(k)和X2(K)继续用这样的办法分解下去,最终一个N点的DFT可以用一组2点的DFT去计算。所以在基2的FFT中,总共有log2(N)级运算,每级有N/2个2点FFT的蝶形运算。

这里我想说得再清晰一些,可是感觉自己表诉得不够准确,害怕误导大家,所以计算过程希望读者还是自己算一遍下来,去感受一下,相信你们会有感悟,如有问题,欢迎大家和在评论区交流。

单个蝶形运算图示

N=8时,时间抽取FFT的运算流图如下:

直接计算DFT和基2FFT算法运算量比较

直接计算DFT由公式:

它的运算量为N^2。

根据FFT的运算规则可见,当N=2^M时,从运算流图可以看出,应该会用M级蝶形,每一级都有N/2个蝶形运算构成。因此,每一级都需要N/2复数乘和N次复数加(每个蝶形需要两次复数加法),所以:

M级运算总的复数乘次数为: MN/2 = log2(N)(N/2)
复数的加次数为: M
N = Nlog2(N)

当 N>>1 时,N^2>>og2(N)*(N/2),所以基2FFT算法比直接计算DFT得到运算量大大减少。

matlab程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
A=input('请输入序列A:');
N=input('请输入N:');
L=length(A(:)); %输入序列长度
b=0;
if(N<L)
P=input('由于N小于序列长度,采样信号不完整,计算无意义,建议重新输入N\n 1:重新输入N\n 0:坚持计算\n');
if(P==1)
N=input('请输入N:');
while(N<L)
N=input('N仍小于序列长度,请重新输入:\n');
end
end
end

while 1
M=log2(N); %M是层数
if(M==fix(M)) % 判断是否适合做快速傅里叶变换
if( log2(L) < M)
while(log2(L+b)~=M)
b=b+1;
end
for n=L+1:1:L+b
A(n)=0;
end
L=length(A);
end
X0=myfft(A,M);
X1=fft(A,N);
X2=myfft(A,M)-fft(A,N);
%%%% 画图
subplot(4,1,1);
stem(A);
title('A序列');

subplot(4,1,2);
stem(X0);
title('A序列的myfft');

subplot(4,1,3);
stem(X1);
title('A序列的fft');

subplot(4,1,4);
stem(X2);
title('myfft-fft');
%%%%%
flag=1;
break;
else
y=input('N点不适合做快速傅里叶变换,现在您有两个选择\n 0:重新输入N\n 1:使用直接傅里叶变换\n');
if(y==0)
N=input('请重新输入N\n');
if(N<L)
P=input('由于N小于序列长度,采样信号不完整,计算无意义,建议重新输入N\n 1:重新输入N\n 0:坚持计算\n');
if(P==1)
N=input('请输入N:');
while(N<L)
N=input('N仍小于序列长度,请重新输入:\n');
end
end
end
else
if(y==1)
X1=fft(A,N);
subplot(2,1,1);
stem(A);
title('A序列');
subplot(2,1,2);
stem(X1);
title('A序列的DFT');
break;
else
fprintf('拜托!只有0和1两个选择,— —#\n');
end
end
end
end

FFT函数文件:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
%%  myfft函数,此程序是用来调用的

function [A] = myfft(A,M)
N=2^M; % M 表示层数
LH=N/2;
J=LH;
N1=N-2;
%倒序程序
for I=1:1:N1
if I<J
T=A(I+1);
A(I+1)=A(J+1);
A(J+1)=T;
end
K=LH;
if J>=K
J=J-K;
K=K/2;
end
J=J+K;
end

for L=1:1:M
B=2^(L-1);
for J=0:B-1
p=J*2^(M-L);
for k=J:2^L:N-1
T=A(k+1)+A(k+B+1)*exp(-1i*2*pi*p/N);
A(k+B+1)=A(k+1)-A(k+B+1)*exp(-1i*2*pi*p/N);
A(k+1)=T;
end
end
end

B=(0);
for i=1:2^M
B(i)=A(i);
end
A=B;

这里说下程序的特点,由于快速傅里叶变换本身因为计算的原因,要求输入的N应该为2^M(M是级数)即类似2,4,8,16…这些数字,所以当输入N=3,6,12时,我们并不能进行快速傅里叶变换,不过在程序中我加入了修改提示,还顺便卖萌装年轻,算是比较人性化的一部分。

此外,还有输入序列x(n)的长度和N之间的关系处理:
当N > 数列长度时,说明采样点过多,需要我们为序列后面补零,这点在程序中我也添加了。
当N < 数列长度,说明采样不足,实际上这样的采样对输入的信号绝对是会有损耗的,即便计算出来也没有太大的意义,除非有专门的算法做处理,类似一些压缩技术,不然计算的结果只是数学上的数字,对实际并没有多大的用途。

当输入一个周期的正弦序列A,信号长度为32,N=32时,效果如下图。

关于程序或者文章的问题,欢迎和我在评论区和我讨论,有打赏就更好了。如文章有错误,请斧正。

参考文献

数字信号处理(高西全,丁美玉).西安电子科技大学出版社
手把手教你理解(FFT).百度文库
还有好几篇不记得了在哪找的PPT(—— ——#)

本文完。

请我喝杯咖啡~
------ 本文结束------
0%