数字信号处理实验 matlab版 快速傅里叶变换(FFT)

(完美格式版，本人自己完成，所有语句正确，不排除极个别错误，特别适用于山大，勿用冰点等工具下载，否则下载之后的word格式会让很多部分格式错误，谢谢)

XXXX学号姓名处XXXX

一、实验目的

1、加深对双线性变换法设计IIR数字滤波器基本方法的了解。

2、掌握用双线性变换法设计数字低通、高通、带通、带阻滤波器的方法。

3、了解MATLAB有关双线性变换法的子函数。

二、实验内容

1、双线性变换法的基本知识

2、用双线性变换法设计IIR数字低通滤波器

3、用双线性变换法设计IIR数字高通滤波器

4、用双线性变换法设计IIR数字带通滤波器

三、实验环境

      MATLAB7.0  

4、实验原理

1、实验涉及的MATLAB子函数

　　（1）fft

功能：一维快速傅里叶变换(FFT)。

调用格式：

；利用FFT算法计算矢量x的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。当x的长度为2的幂次方时，则fft函数采用基2的FFT算法，否则采用稍慢的混合基算法。

；采用n点FFT。当x的长度小于n时，fft函数在x的尾部补零，以构成n点数据；当x的长度大于n时，fft函数会截断序列x。当x为矩阵时，fft函数按类似的方式处理列长度。

（2）ifft

　　功能：一维快速傅里叶逆变换(IFFT)。

　　调用格式：

；用于计算矢量x的IFFT。当x为矩阵时，计算所得的y为矩阵x中每一列的IFFT。

；采用n点IFFT。当length(x)n时，将x截断，使length(x)＝n。

（3）fftshift

功能：对fft的输出进行重新排列，将零频分量移到频谱的中心。

调用格式：

；对fft的输出进行重新排列，将零频分量移到频谱的中心。当x为向量时，直接将x中的左右两半交换而产生y。

当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所示。

　图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];    %建立x2(n)序列

>> N=length(xn1);

>> 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计算的结果一致。

3、用FFT计算有限长序列的频谱

（1）基本概念

一个序号从n1到n2的时域有限长序列x(n)，它的频谱X(ejw)定义为它的离散傅里叶变换，且在奈奎斯特(Nyquist)频率范围内有界并连续。序列的长度为N，则N＝n2－n1＋1。计算x(n)的离散傅里叶变换(DFT)得到的是X(ejw)的N个样本点X(ejwk)。其中数字频率为

式中：dw为数字频率的分辨率；k取对应－(N－1)/2到(N－1)/2区间的整数。

在实际使用中，往往要求计算出信号以模拟频率为横坐标的频谱，此时对应的模拟频率为

式中：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所示。

图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)];

>> 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));        %横轴化成模拟频率作幅度谱

>> title('幅度频谱');xlabel('rad/s');

>> subplot(1,2,2);plot(k*D,angle(X));        %横轴化成模拟频率作相位谱

>> title('相位频谱');xlabel('rad/s');

此时程序执行的结果如图14-3所示。由图可以看出，图形的分辨率提高，曲线几乎是连续的频谱了。

图14-3 将例14-2有限长序列末尾补0到N＝1000时的频谱

(3)实偶序列如何补0

例14-4 已知一个矩形窗函数序列为

采样周期Ts＝0.5 s，要求用FFT求其频谱。

解 由于该序列是一个实的偶序列，因而补0时需要仔细分析。假定按N＝32补0，则主值区域在n＝0～31，FFT的输入应为

即原来n=[－5：－1]的前五个点移到n＝[27：31]中去了。

下面考虑分别用N=32，，512，观察不同N值代入对频谱的影响。

程序如下，

>> Ts=0.5;C=[32,,512]; %输入不同的N值

>> for r=0:2;

>> N=C(r+1);

>> xn=[ones(1,6),zeros(1,N-11),ones(1,5)];             %建立x(n) 

>> D=2*pi/(N*Ts);

>> k=floor(-(N-1)/2:(N-1)/2);

>> X=fftshift(fft(xn,N));

>> subplot(3,2,2*r+1);plot(k*D,abs(X));         %幅度频谱

>> subplot(3,2,2*r+2);stairs(k*D,angle(X));        %相位频谱

>> end

注意：此处相位频谱使用了stairs，因为该相位频谱变化率比较陡峭。

　　程序执行结果如图14-4所示。

图14-4 将例14-4有限长序列补0到N＝32、、512时的频谱

如果将x(n)的输入写成

xn=[ones(1,11),zeros(1,N-11)];%建立x(n－5)

相当于起点不是取自n＝0而是n＝－5，计算的是x(n－5)的频谱。幅度频谱不受影响，相位频谱引入一个线性相位－5w，如图14-5所示。

图14-5 将有限长位移序列x(n－5)补0到N＝512时的频谱

4、用FFT计算无限长序列的频谱

用FFT进行无限长序列的频谱计算，首先要将无限长序列截断成一个有限长序列。序列长度的取值对频谱有较大的影响，带来的问题是引起频谱的泄漏和波动。

例14-5 已知一个无限长序列为

采样频率Fs＝20 Hz，要求用FFT求其频谱。

解 MATLAB程序如下：

>> Fs=20;C=[8,16,128]; %输入不同的N值

>> for r=0:2;

>> N=C(r+1);

>> n=0:N-1;

>> xn=exp(-0.5*n);%建立x(n) 

>> D=2*pi*Fs/N;

>> k=floor(-(N-1)/2:(N-1)/2);

>> X=fftshift(fft(xn,N));

>> subplot(3,2,2*r+1);plot(k*D,abs(X));

>> axis([-80,80,0,3]);

>> subplot(3,2,2*r+2);stairs(k*D,angle(X));

>> axis([-80,80,-1,1]);

>> end

运行结果如图14-6所示。

图14-6 将无限长序列截断为N＝8，16，128时的频谱

由图14-6可见，N值取得越大，即序列保留得越长，曲线精度越高。

例14-6 用FFT计算下列连续时间信号的频谱，并观察选择不同的Ts和N值对频谱特性的影响。

　　xa(t)＝e－0.01t(sin2t＋sin2.1t＋sin2.2t)  t≥0

解 该题选择了三个非常接近的正弦信号，为了将各频率成分区分出来，在满足奈奎斯特定理的条件下确定采样周期，选择三组数据，分别是Ts＝0.5 s、0.25 s和0.125 s；再确定N值，分别选择N＝256和2048。观察不同Ts和N的组合对频谱的影响。

程序如下：

>> T0=[0.5,0.25,0.125,0.125]; %输入不同的Ts值

>> N0=[256,256,256,2048];%输入不同的N值

>> for r=1:4;

>> Ts=T0(r);N=N0(r);%赋Ts和N值 

>> n=0:N-1;

>> D=2*pi/(Ts*N);%计算模拟频率分辨率 

>> xa=exp(-0.01*n*Ts).*(sin(2*n*Ts)+sin(2.1*n*Ts)+sin(2.2*n*Ts));

>> k=floor(-(N-1)/2:(N-1)/2);

>> Xa=Ts*fftshift(fft(xa,N));

>> [r,Xa(1)]%输出Xa(1)的数值,供误差计算用 

>> subplot(2,2,r);plot(k*D,abs(Xa),'k');

>> axis([1,3,1.1*min(abs(Xa)),1.1*max(abs(Xa))]);

>> end

运行结果如图14-7所示。

图14-7 用FFT计算三个很靠近的谐波分量的频谱图

由图14-7可以得出以下结论：

N同样取256(如前三个图形)，当Ts越大时，时域信号的长度L＝NTs保留得越长，则分辨率越高，频谱特性误差越小；反之，则分辨率越低，频谱特性误差越大，甚至丢失某些信号分量。

　　Ts相同(如后两个图形)，当N越大时，在［0，2p］范围内等间隔抽样点数越多，且时域信号的长度L＝NTs保留得越长，则分辨率越高，频谱特性误差越小；反之，当N越小时，在［0，2p］范围内等间隔抽样点数越少，则有可能漏掉某些重要的信号分量，称为栅栏效应。

5、实验过程

1 已知有限长序列x(n)＝［1，0.5，0，0.5，1，1，0.5，0］，要求：

(1)用FFT算法求该时域序列的DFT、IDFT的图形；

(2)假定采样频率Fs＝20 Hz，序列长度N分别取8、32和，使用FFT来计算其幅度频谱和相位频谱。

解 MATLAB程序如下：

>> xn=[1,0.5,0,0.5,1,1,0.5,0];

>> N=length(xn);

>> n=0:N-1;k=0:N-1;

>> Xk=fft(xn,N);

>> subplot(2,1,1);stem(k,abs(Xk));

>> title('Xk=DFT(x(n))');

>> xn1=ifft(Xk,N);

>> subplot(2,1,2);stem(n,xn1);

>> title('x(n)=IDFT(Xk)');

　　运行结果如图1所示。

图1

>> Fs=20;C=[8,32,];

>> for r=0:2;

>> N=C(r+1);

>> xn=[1,0.5,0,0.5,1,1,0.5,0];

>> D=2*pi*Fs/N;

>> k=floor(-(N-1)/2:(N-1)/2);

>> X=fftshift(fft(xn,N));

>> subplot(3,2,2*r+1);plot(k*D,abs(X));

>> title('幅度频谱');xlabel('rad/s');

>> subplot(3,2,2*r+2);stairs(k*D,angle(X));

>> title('相位频谱');xlabel('rad/s');

>> end

　　运行结果如图2所示。

图2

2 已知一个无限长序列x(n)＝0.5n(n≥0)，采样周期Ts＝0.2 s，要求序列长度N分别取8、32和，用FFT求其频谱。

解 MATLAB程序如下：

>> Ts=0.2;C=[8,32,];

>> for r=0:2;

>> N=C(r+1);

>> n=0:N-1;

>> xn=exp(0.5*n);

>> D=2*pi/(N*Ts);

>> k=floor(-(N-1)/2:(N-1)/2);

>> X=fftshift(fft(xn,N));

>> subplot(3,2,2*r+1);plot(k*D,abs(X));

>> title('幅度频谱');xlabel('rad/s');

>> subplot(3,2,2*r+2);stairs(k*D,angle(X));

>> title('相位频谱');xlabel('rad/s');

>> end

　　运行结果如图3所示。

图3

六、实验感想

    通过此次实验中练习使用matlab语言进行快速傅里叶变换，更为熟悉的掌握了matlab的功能，在实验过程中也遇到很多小问题，并通过仔细检查和查阅相关书籍解决此类问题，让我深刻认识到，细节的重要性。在使用help过程中，深切体会到良好的英语基础和充实的课堂知识的重要性。