数字信号处理——快速傅里叶变换

队长-Leader

已于 2023-04-12 10:47:14 修改

阅读量1.1w

点赞数 12

分类专栏：算法 MATLAB 文章标签：信号处理 matlab 算法

于 2021-05-13 23:09:23 首次发布

本文链接：https://blog.csdn.net/k331922164/article/details/116770519

版权

MATLAB 同时被 2 个专栏收录

17 篇文章 16 订阅

订阅专栏

算法

15 篇文章 17 订阅

订阅专栏

本文详细介绍了快速傅里叶变换（FFT）的基本原理，包括旋转因子、蝶形图和倒序等关键步骤，并通过具体的Matlab代码展示了FFT的实现过程。同时，对比了自定义实现的FFT与Matlab内置fft函数的计算结果，证明了算法的正确性。此外，还探讨了如何利用复序列FFT计算实序列FFT和逆FFT。

摘要由CSDN通过智能技术生成

一、前言。

快速傅里叶变换不是一种新的变换，而是离散傅里叶变换的快速算法，而且这个快速算法有很多种，都统称为快速傅里叶变换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=%d\r\n',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=%d\r\n',j,k);
		while(k<N)
			% fprintf('k=%d,%d;p=%d\r\n',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=%d\r\n',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=%d\r\n',j,k);
			while(k<N)
				% fprintf('k=%d,%d;p=%d\r\n',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=%d\r\n',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=%d\r\n',j,k);
			while(k<N)
				% fprintf('k=%d,%d;p=%d\r\n',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=%d\r\n',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=%d\r\n',j,k);
			while(k<N)
				% fprintf('k=%d,%d;p=%d\r\n',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