概述
一、前言。
快速傅里叶变换不是一种新的变换，而是离散傅里叶变换的快速算法，而且这个快速算法有很多种，都统称为快速傅里叶变换FFT。
如果直接用公式计算DFT，其时间复杂度为O(n*n)，这是难以应用在工程当中的。
而基2-FFT的时间复杂度为O(n*log(n))，随着序列的点数增加，其运算效率大大提高。
其中最常用的是基2时域抽取的FFT，下面来详细说明。
二、旋转因子。
function [] = calW() %只需要计算N/2点的Wn
clear;close all;clc;
N=8;
for j=1:N
Wn(j) = exp(-i * (2*pi/N))^(j-1); %这里的i为虚数单位
end
for j=1:N
disp(Wn(j))
end
这里可以看出旋转因子是对称的，Wn(1)=1而Wn(5)为-1，这里的位置从1开始（Matlab语法）。
所以只需要得到N/2个点的旋转因子，也就是一半的点数，而另一半乘以-1即可得到。
旋转因子的性质，如下图所示。
其中可约性：下标N和次方nk，同乘以一个数m，旋转因子的值不变（下面会用到）。
三、蝶形图。
蝶形运算如下图所示，其中只需要x2(k)乘以Wn^k，再乘以-1可以得到另一个旋转因子。
8点时间抽取的基2-fft，蝶形图如下图所示。
N=8，那么就有m=log2(N)=3级（竖着看有三级蝶形运算），而每级蝶形运算又包含了N/2=4个蝶形运算。
其中Wn^0=1，为了方便理解才标出来。
第二级中，WN^0、WN^2，这里N=8，实际为W4^0和W4^1，利用可约性，下标和次方同乘以2了，为了方便代码实现。
四、倒序。
我们看到时域的位置并不是从0排到7的。这里可以先将0-7转为2进制，进制的位宽为log2(N)=3，再将2进制整体镜像过来，即可完成倒序。
function [] = reverseIndex() %倒序
clear;close all;clc;
len = 8;
yk = zeros(1,len);
xk=yk;
for i=1:len
indexBin = fliplr(dec2bin(i-1,log2(len)));
xk(i)=i-1;
yk(i) = bin2dec(indexBin);
end
xk
yk
五、FFT算法的Matlab实现。
function [] = fft_base_2()
% 基2时间抽取fft
clear;close all;clc;
% Xm = [1,1,1,1,0,0,0,0]; %8点时域
Xm = [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0]; %16点时域
subplot(211)
stem(Xm);
N = length(Xm); %FFT的点数
Wn = calW(N); %旋转因子,先计算好,后面查表即可
Xm = reverseIndex(Xm); %倒序
Xk = Xm; %频域
m=log2(N);
for i=1:m %log2(N)级运算,8点有3级运算
distNum = 2^(i-1); %级内蝶形运算的距离
% fprintf('distNum=%drn',distNum);
for j=0:distNum-1 %级内Wn的次方,此时wn的n为原始的2,4,8
%p为Wn的次方,也是数组的元素位置
p=j*(2^(m-i))+1; %将wn的n和次方都乘以一个数达到N,此时wn的n均为8,得到的次方p,+1是因为matlab的数组从1开始
k=j+1; %+1是因为matlab的数组从1开始
% fprintf('j=%d,k=%drn',j,k);
while(k<N)
% fprintf('k=%d,%d;p=%drn',k,k+distNum,p);
[Xk(k),Xk(k+distNum)] = butterfly(Wn(p),Xk(k),Xk(k+distNum)); %蝶形运算
k=k+2^i; %每次累加级内蝶形运算的距离
end
end
end
Xk
subplot(212)
stem(abs(Xk));
hold on
plot(abs(Xk));
function Wn = calW(N) %只需要计算N/2点的Wn
Wn(1) = 1; %Wn的0次方为1
for j=2:N/2
Wn(j) = exp(-i * (2*pi/N))^(j-1); %这里的i为虚数单位
end
function [y1k,y2k] = butterfly(wnk,x1k,x2k) %计算蝶形运算
tmp = wnk*x2k;
y1k = x1k + tmp;
y2k = x1k - tmp;
function yk = reverseIndex(xk) %倒序
len = length(xk);
yk = zeros(1,len);
for i=1:len
indexBin = fliplr(dec2bin(i-1,log2(len)));
yk(i) = xk(bin2dec(indexBin) + 1);
end
Xm = [1,1,1,1,0,0,0,0]; 时，得到的幅频特性为：
Xm = [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0]; 时，得到的幅频特性为：
六、与Matlab的fft函数比较。
tic
Xm = [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0]; %16点时域
Xk_fft=fft(Xm)
toc时间已过 0.001903 秒。
Xk_fft =
1 至 7 列
4.0000 + 0.0000i 3.0137 - 2.0137i 1.0000 - 2.4142i -0.2483 - 1.2483i 0.0000 + 0.0000i 0.8341 + 0.1659i 1.0000 - 0.4142i
8 至 14 列
0.4005 - 0.5995i 0.0000 + 0.0000i 0.4005 + 0.5995i 1.0000 + 0.4142i 0.8341 - 0.1659i 0.0000 + 0.0000i -0.2483 + 1.2483i
15 至 16 列
1.0000 + 2.4142i 3.0137 + 2.0137i
时间已过 0.001216 秒。
Xk =
1 至 7 列
4.0000 + 0.0000i 3.0137 - 2.0137i 1.0000 - 2.4142i -0.2483 - 1.2483i 0.0000 + 0.0000i 0.8341 + 0.1659i 1.0000 - 0.4142i
8 至 14 列
0.4005 - 0.5995i 0.0000 + 0.0000i 0.4005 + 0.5995i 1.0000 + 0.4142i 0.8341 - 0.1659i 0.0000 + 0.0000i -0.2483 + 1.2483i
15 至 16 列
1.0000 + 2.4142i 3.0137 + 2.0137i可以看出，两者计算结果一致。
七、利用N点复序列FFT，计算2个N点实序列FFT。
function [] = fft_test1()
clear;close all;clc;
Xm1 = [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0];
Xm2 = [1,2,3,2,1,0,0,0,0,0,0,0,0,0,0,0];
Ym = Xm1 + Xm2*i;
Yk1 = fft_base_2(Ym);
N = length(Yk1);
Yk2 = zeros(1,N);
for m=0:N-1
Yk2(m+1) = conj(Yk1(mod(N-m,N)+1));
end
Yk2;
Xk1 = (Yk1 + Yk2)/2;
Xk1_fft=fft(Xm1); %与matlab的fft函数相比
Xk2 = (Yk1 - Yk2)/(2*i);
Xk2_fft=fft(Xm2);
subplot(611)
stem(Xm1);
subplot(612)
stem(Xm2);
subplot(613)
stem(abs(Xk1));
hold on
plot(abs(Xk1));
subplot(614)
stem(abs(Xk1_fft));
hold on
plot(abs(Xk1_fft));
subplot(615)
stem(abs(Xk2));
hold on
plot(abs(Xk2));
subplot(616)
stem(abs(Xk2_fft));
hold on
plot(abs(Xk2_fft));
function Xk = fft_base_2(Xm) % 基2时间抽取fft
N = length(Xm); %FFT的点数
Wn = calW(N); %旋转因子,先计算好,后面查表即可
Xm = reverseIndex(Xm); %倒序
Xk = Xm; %频域
m=log2(N);
for i=1:m %log2(N)级运算,8点有3级运算
distNum = 2^(i-1); %级内蝶形运算的距离
% fprintf('distNum=%drn',distNum);
for j=0:distNum-1 %级内Wn的次方,此时wn的n为原始的2,4,8
%p为Wn的次方,也是数组的元素位置
p=j*(2^(m-i))+1; %将wn的n和次方都乘以一个数达到N,此时wn的n均为8,得到的次方p,+1是因为matlab的数组从1开始
k=j+1; %+1是因为matlab的数组从1开始
% fprintf('j=%d,k=%drn',j,k);
while(k<N)
% fprintf('k=%d,%d;p=%drn',k,k+distNum,p);
[Xk(k),Xk(k+distNum)] = butterfly(Wn(p),Xk(k),Xk(k+distNum)); %蝶形运算
k=k+2^i; %每次累加级内蝶形运算的距离
end
end
end
function Wn = calW(N) %只需要计算N/2点的Wn
Wn(1) = 1; %Wn的0次方为1
for j=2:N/2
Wn(j) = exp(-i * (2*pi/N))^(j-1); %这里的i为虚数单位
end
function [y1k,y2k] = butterfly(wnk,x1k,x2k) %计算蝶形运算
tmp = wnk*x2k;
y1k = x1k + tmp;
y2k = x1k - tmp;
function yk = reverseIndex(xk) %倒序
len = length(xk);
yk = zeros(1,len);
for i=1:len
indexBin = fliplr(dec2bin(i-1,log2(len)));
yk(i) = xk(bin2dec(indexBin) + 1);
end
八、利用N点复序列FFT，计算2N点实序列FFT。
function [] = fft_test2()
clear;close all;clc;
XmIn = [1,-1,1,-1,2,-1,1,-1];
N2 = length(XmIn);
Xm1 = zeros(1,N2/2);
Xm2 = Xm1;
for m=0:N2/2-1 %先以基2抽取
Xm1(m+1) = XmIn(2*m + 1);
Xm2(m+1) = XmIn(2*m + 1 + 1);
end
Ym = Xm1 + Xm2*i;
Yk1 = fft_base_2(Ym);
N = length(Yk1);
Yk2 = zeros(1,N);
for m=0:N-1
Yk2(m+1) = conj(Yk1(mod(N-m,N)+1));
end
Xk1 = (Yk1 + Yk2)/2;
Xk2 = (Yk1 - Yk2)/(2*i);
Xk_my=zeros(1,2*N);
Wn = calW(2*N);
for m=0:N-1 %以基2合成
[Xk_my(m+1),Xk_my(m+N+1)] = butterfly(Wn(m+1),Xk1(m+1),Xk2(m+1));
end
Xk_my
Xk_fft=fft(XmIn) %与matlab的fft函数相比
subplot(311)
stem(XmIn);
subplot(312)
stem(abs(Xk_my));
hold on
plot(abs(Xk_my));
subplot(313)
stem(abs(Xk_fft));
hold on
plot(abs(Xk_fft));
function Xk = fft_base_2(Xm) % 基2时间抽取fft
N = length(Xm); %FFT的点数
Wn = calW(N); %旋转因子,先计算好,后面查表即可
Xm = reverseIndex(Xm); %倒序
Xk = Xm; %频域
m=log2(N);
for i=1:m %log2(N)级运算,8点有3级运算
distNum = 2^(i-1); %级内蝶形运算的距离
% fprintf('distNum=%drn',distNum);
for j=0:distNum-1 %级内Wn的次方,此时wn的n为原始的2,4,8
%p为Wn的次方,也是数组的元素位置
p=j*(2^(m-i))+1; %将wn的n和次方都乘以一个数达到N,此时wn的n均为8,得到的次方p,+1是因为matlab的数组从1开始
k=j+1; %+1是因为matlab的数组从1开始
% fprintf('j=%d,k=%drn',j,k);
while(k<N)
% fprintf('k=%d,%d;p=%drn',k,k+distNum,p);
[Xk(k),Xk(k+distNum)] = butterfly(Wn(p),Xk(k),Xk(k+distNum)); %蝶形运算
k=k+2^i; %每次累加级内蝶形运算的距离
end
end
end
function Wn = calW(N) %只需要计算N/2点的Wn
Wn(1) = 1; %Wn的0次方为1
for j=2:N/2
Wn(j) = exp(-i * (2*pi/N))^(j-1); %这里的i为虚数单位
end
function [y1k,y2k] = butterfly(wnk,x1k,x2k) %计算蝶形运算
tmp = wnk*x2k;
y1k = x1k + tmp;
y2k = x1k - tmp;
function yk = reverseIndex(xk) %倒序
len = length(xk);
yk = zeros(1,len);
for i=1:len
indexBin = fliplr(dec2bin(i-1,log2(len)));
yk(i) = xk(bin2dec(indexBin) + 1);
end
九、利用N点复序列FFT，计算N点复序列IFFT。
function [] = fft_test3()
clear;close all;clc;
Xm = [1,1,1,1,0,0,0,0];
N = length(Xm);
Yk = conj(fft_base_2(conj(Xm)))/N
Yk_ifft=ifft(Xm) %与matlab的ifft函数相比
subplot(311)
stem(Xm);
subplot(312)
stem(abs(Yk));
hold on
plot(abs(Yk));
subplot(313)
stem(abs(Yk_ifft));
hold on
plot(abs(Yk_ifft));
function Xk = fft_base_2(Xm) % 基2时间抽取fft
N = length(Xm); %FFT的点数
Wn = calW(N); %旋转因子,先计算好,后面查表即可
Xm = reverseIndex(Xm); %倒序
Xk = Xm; %频域
m=log2(N);
for i=1:m %log2(N)级运算,8点有3级运算
distNum = 2^(i-1); %级内蝶形运算的距离
% fprintf('distNum=%drn',distNum);
for j=0:distNum-1 %级内Wn的次方,此时wn的n为原始的2,4,8
%p为Wn的次方,也是数组的元素位置
p=j*(2^(m-i))+1; %将wn的n和次方都乘以一个数达到N,此时wn的n均为8,得到的次方p,+1是因为matlab的数组从1开始
k=j+1; %+1是因为matlab的数组从1开始
% fprintf('j=%d,k=%drn',j,k);
while(k<N)
% fprintf('k=%d,%d;p=%drn',k,k+distNum,p);
[Xk(k),Xk(k+distNum)] = butterfly(Wn(p),Xk(k),Xk(k+distNum)); %蝶形运算
k=k+2^i; %每次累加级内蝶形运算的距离
end
end
end
function Wn = calW(N) %只需要计算N/2点的Wn
Wn(1) = 1; %Wn的0次方为1
for j=2:N/2
Wn(j) = exp(-i * (2*pi/N))^(j-1); %这里的i为虚数单位
end
function [y1k,y2k] = butterfly(wnk,x1k,x2k) %计算蝶形运算
tmp = wnk*x2k;
y1k = x1k + tmp;
y2k = x1k - tmp;
function yk = reverseIndex(xk) %倒序
len = length(xk);
yk = zeros(1,len);
for i=1:len
indexBin = fliplr(dec2bin(i-1,log2(len)));
yk(i) = xk(bin2dec(indexBin) + 1);
end
最后
以上就是笑点低星星为你收集整理的数字信号处理——快速傅里叶变换的全部内容，希望文章能够帮你解决数字信号处理——快速傅里叶变换所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。
本图文内容来源于网友提供，作为学习参考使用，或来自网络收集整理，版权属于原作者所有。
发表评论 取消回复