发布时间 : 星期一 文章数字信号处理实验 matlab版 快速傅里叶变换(FFT)更新完毕开始阅读7cc9b6d776a20029bd642d48
实验14 快速傅里叶变换(FFT)
(完美格式版,本人自己完成,所有语句正确,不排除极个别错误,特别适用于山大,勿用冰点等工具下载,否则下载之后的word格式会让很多部分格式错误,谢谢)
XXXX学号姓名处XXXX
一、实验目的
1、加深对双线性变换法设计IIR数字滤波器基本方法的了解。
2、掌握用双线性变换法设计数字低通、高通、带通、带阻滤波器的方法。 3、了解MATLAB有关双线性变换法的子函数。
二、实验内容
1、双线性变换法的基本知识
2、用双线性变换法设计IIR数字低通滤波器 3、用双线性变换法设计IIR数字高通滤波器 4、用双线性变换法设计IIR数字带通滤波器
三、实验环境
MATLAB7.0
四、实验原理
1、实验涉及的MATLAB子函数
(1)fft
功能:一维快速傅里叶变换(FFT)。 调用格式:
y?fft(x);利用FFT算法计算矢量x的离散傅里叶变换,当x为矩阵时,y为矩阵x
每一列的FFT。当x的长度为2的幂次方时,则fft函数采用基2的FFT算法,否则采用稍慢的混合基算法。
y?fft(x,n);采用n点FFT。当x的长度小于n时,fft函数在x的尾部补零,以构成n
点数据;当x的长度大于n时,fft函数会截断序列x。当x为矩阵时,fft函数按类似的方式处理列长度。
(2)ifft
功能:一维快速傅里叶逆变换(IFFT)。 调用格式:
y?ifft(x);用于计算矢量x的IFFT。当x为矩阵时,计算所得的y为矩阵x中每一
列的IFFT。
;采用n点IFFT。当length(x)
将x截断,使length(x)=n。
(3)fftshift
功能:对fft的输出进行重新排列,将零频分量移到频谱的中心。 调用格式:
y?fftshift(x)y?ifft(x,n);对fft的输出进行重新排列,将零频分量移到频谱的中心。当x为向
量时,fftshift(x)直接将x中的左右两半交换而产生y。
当x为矩阵时,fftshift(x)同时将x的左右、上下进行交换而产生y。
2、用MATLAB提供的子函数进行快速傅里叶变换
从理论学习可知,DFT是唯一在时域和频域均为离散序列的变换方法,它适用于有限长序列。尽管这种变换方法是可以用于数值计算的,但如果只是简单的按照定义进行数据处理,当序列长度很大时,则将占用很大的内存空间,运算时间将很长。
快速傅里叶变换是用于DFT运算的高效运算方法的统称,FFT只是其中的一种。FFT主要有时域抽取算法和频域抽取算法,基本思想是将一个长度为N的序列分解成多个短序列,如基2算法、基4算法等,大大缩短了运算的时间。
MATLAB中提供了进行快速傅里叶变换(FFT)的子函数,用fft计算DFT,用ifft计算IDFT。
例14-1 已知一个长度为8点的时域离散信号,n1=0,n2=7,在n0=4前为0,n0以后为1。对其进行FFT变换,作时域信号及DFT、IDFT的图形。
解 MATLAB程序如下:
>> n1=0;n2=7;n0=4;
>> n=n1:n2;
>> N=length(n);
>> xn=[(n-n0)>=0];%建立时域信号 >> subplot(2,2,1); >> Stem(n,xn); >> title('x(n)');
>> k=0:N-1;
>> Xk=fft(xn,N);%用FFT计算信号的DFT >> subplot(2,1,2); >> Stem(k,abs(Xk));
>> title('Xk=DFT((n))');
>> xn1=ifft(Xk,N);%用IFFT计算信号的IDFT
>> subplot(2,2,2);stem(n,xn1); >> title('x(n)=IDFT(Xk)');
运行结果如图14-1所示。
x(n)11x(n)=IDFT(Xk)0.50.504050100Xk=DFT((n))5102001234567
图14-1 例14-1用FFT求有限长序列的傅里叶变换
例14-2 将例13-5已知的两个时域周期序列分别取主值,得到x1=[1,1,1,0,0,0],x2=[0,1,2,3,0,0],求时域循环卷积y(n)并用图形表示。
解 本例将例13-5使用DFT处理的计算,改为用FFT和IFFT进行循环卷积。
程序如下:
>> xn1=[0,1,2,3,0,0]; %建立x1(n)序列
>> xn2=[1,1,1,0,0,0]; >> N=length(xn1);
%建立x2(n)序列
>> n=0:N-1;k=0:N-1;
>> Xk1=fft(xn1,N);%由x1(n)的FFT求X1(k) >> Xk2=fft(xn2,N);%由x2(n)的FFT求X2(k) >> Yk=Xk1.*Xk2;%Y(k)=X1(k)X2(k)
>> yn=ifft(Yk,N);%由Y(k)的IFFT求y(n)
>> yn=abs(yn); >> stem(n,yn);
运行结果如图所示,与例13-5用DFT计算的结果一致。
6420012345
3、用FFT计算有限长序列的频谱
(1)基本概念
一个序号从n1到n2的时域有限长序列x(n),它的频谱X(ejw)定义为它的离散傅里叶变换,且在奈奎斯特(Nyquist)频率范围内有界并连续。序列的长度为N,则N=n2-n1+1。计算x(n)的离散傅里叶变换(DFT)得到的是X(ejw)的N个样本点X(ejwk)。其中数字频率为
ωk?k(2πN)?kdω
式中:dw为数字频率的分辨率;k取对应-(N-1)/2到(N-1)/2区间的整数。
在实际使用中,往往要求计算出信号以模拟频率为横坐标的频谱,此时对应的模拟频率为
Ωk?ωk/Ts?k(2πNTs)?k(2πL)?kD
式中:D为模拟频率的分辨率或频率间隔;Ts为采样信号的周期,Ts=1/Fs;定义信号时域长度L=NTs。
在使用FFT进行DFT的高效运算时,一般不直接用n从n1到n2的x(n),而是取 的主值区间(n=0,1,…,N-1)的数据,经FFT将产生N个数据,定位在k=0,1,…,N-1的数字频率点上,即对应[0,2p]。如果要显示[-p,p]范围的频谱,则可以使用fftshift(X)进行位移。
(2)频谱的显示及分辨率问题
例14-3 已知有限长序列x(n)=[1,2,3,2,1],其采样频率Fs=10 Hz。请使用FFT计算其频谱。
解 MATLAB程序如下:
>> Fs=10;
>> xn=[1,2,3,2,1];N=length(xn);
>> D=2*pi*Fs/N; %计算模拟频率分辨率 >> k=floor(-(N-1)/2:(N-1)/2);
%频率显示范围对应[-p,p]
>> X=fftshift(fft(xn,N)); %作FFT运算且移位p >> subplot(1,2,1);plot(k*D,abs(X),'o:'); %横轴化成模拟频率作幅度谱 >> title('幅度频谱');xlabel('rad/s');
>> subplot(1,2,2);plot(k*D,angle(X),'o:'); >> title('相位频谱');xlabel('rad/s'); 程序运行结果:
absX=
0.3820 2.6180 9.0000 2.6180 0.3820
angleX=
-1.2566 2.5133 0 -2.5133 1.2566
运行结果如图14-2所示。
幅度频谱104250-20-500rad/s50-4-500rad/s50相位频谱 %横轴化成模拟频率作相位谱
图14-2 例14-3有限长序列的频谱
由图14-2可知,当有限长序列的长度N=5时,频谱的频率样本点数也为5,如图上用“。”表示的点位。频率点之间的间距非常大,即分辨率很低。即使使用了plot命令的插值功能,显示出的曲线仍是断断续续的,与真实曲线有较大的误差。
改变分辨率的基本方法是给输入序列补零,即增加频谱的密度。注意,这种方法只是改善了图形的视在分辨率,并不增加频谱的细节信息。
将上述有限长序列x(n)=[1,2,3,2,1]末尾补0到N=1000点,将程序改为:
>> Fs=10;N=1000;
>> xn=[1,2,3,2,1];Nx=length(xn); >> xn=[1,2,3,2,1,zeros(1,N-Nx-1)];