实验一FFT与谱分析.docx

实验一FFT与谱分析.docx

《实验一FFT与谱分析.docx》由会员分享，可在线阅读，更多相关《实验一FFT与谱分析.docx（24页珍藏版）》请在冰豆网上搜索。

实验一FFT与谱分析

《数字信号处理实验报告》

姓名：

朱道萍

班级：

通信1401

学号：

201403090129

实验一FFT与谱分析

1.实验目的

（1）加深对DFT算法原理和物理意义的理解。

（2）熟悉FFT算法，增强对FFT结果的分析能力。

（3）掌握用FFT对连续时间信号和离散时间信号进行谱分析的方法，了解误差及其产生的原因。

2.实验原理

（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。

（2）复习FFT算法原理与编程思想。

（3）查询Matlab中fft函数的用法，学会通过改变fft的参数计算不同DFT。

（4）编制信号产生子程序，产生一下信号：

（5）编写实验程序：

参考图1-1所示的程序流程图，编写实验程序。

3.实验内容

对2中所给出的11个信号逐一进行谱分析。

各信号的FFT变换区间N、连续时间信号x6（t）的采样频率fs等参数的参考取值如下：

a）x1[n],x2[n],x3[n],x4[n],x5[n],x7[n],x8[n]:

N=8,16

b）x6（t）:

fs=64Hz，N=16,32,64

c）x9[n],x10[n],x11[n]:

N=16

x1[n]:

N=8

n=0:

N-1;

x1=[11110000];

m=0:

2*N-1;

x2=[1111000000000000];

f1=fft（x1,N）;

f2=fft（x2,2*N）;

subplot（2,2,1）

stem（n,x1）

subplot（2,2,2）

stem（n,abs（f1））

subplot（2,2,3）

stem（m,x2）

subplot（2,2,4）

stem（m,abs（f2））

x2[n]:

N=8

n=0:

N-1;

x1=[12344321];

m=0:

2*N-1;

x2=[1234432100000000];

f1=fft（x1,N）;

f2=fft（x2,2*N）;

subplot（2,2,1）

stem（n,x1）

subplot（2,2,2）

stem（n,abs（f1））

subplot（2,2,3）

stem（m,x2）

subplot（2,2,4）

stem（m,abs（f2））

x3[n]:

N=8

n=0:

N-1;

x1=[43211234];

m=0:

2*N-1;

x2=[4321123400000000];

f1=fft（x1,N）;

f2=fft（x2,2*N）

subplot（2,2,1）

stem（n,x1）

subplot（2,2,2）

stem（n,abs（f1））

subplot（2,2,3）

stem（m,x2）

subplot（2,2,4）

stem（m,abs（f2））

x4[n]:

N=8;

n=0:

N-1;

x1=cos（0.25*pi*n）;

m=0:

2*N-1;

x2=cos（0.25*pi*m）;

f1=fft（x1,N）;

f2=fft（x2,2*N）;

subplot（2,2,1）

stem（n,x1）

subplot（2,2,2）

stem（n,abs（f1））

subplot（2,2,3）

stem（m,x2）

subplot（2,2,4）

stem（m,abs（f2））

x5[n]:

N=8;

n=0:

N-1;

x1=sin（0.125*pi*n）;

m=0:

2*N-1;

x2=sin（0.125*pi*m）;

f1=fft（x1,N）;

f2=fft（x2,2*N）;

subplot（2,2,1）

stem（n,x1）

subplot（2,2,2）

stem（n,abs（f1））

subplot（2,2,3）

stem（m,x2）

subplot（2,2,4）

stem（m,abs（f2））

x6[n]:

N=16;

fs=64;

n=0:

N-1;

x1=cos（8*pi*n/fs）+cos（16*pi*n/fs）+cos（20*pi*n/fs）;

m=0:

2*N-1;

x2=cos（8*pi*m/fs）+cos（16*pi*m/fs）+cos（20*pi*m/fs）;

l=0:

4*N-1;

x3=cos（8*pi*l/fs）+cos（16*pi*l/fs）+cos（20*pi*l/fs）;

f1=fft（x1,N）;

f2=fft（x2,2*N）;

f3=fft（x3,4*N）;

subplot（3,2,1）

stem（n,x1）

subplot（3,2,2）

stem（n,abs（f1））

subplot（3,2,3）

stem（m,x2）

subplot（3,2,4）

stem（m,abs（f2））

subplot（3,2,5）

stem（l,x3）

subplot（3,2,6）

stem（l,abs（f3））

x7[n],x8[n]:

N=8;

n=0:

N-1;

x51=sin（0.125*pi*n）;

x41=cos（0.25*pi*n）;

x71=x41+x51;

m=0:

2*N-1;

x52=sin（0.125*pi*m）;

x42=cos（0.25*pi*m）;

x72=x42+x52;

f41=fft（x41,N）;

f51=fft（x51,N）;

f71=fft（x71,N）;

f42=fft（x42,2*N）;

f52=fft（x52,2*N）;

f72=fft（x72,2*N）;

figure

（1）

subplot（2,2,1）

stem（n,real（f71））

subplot（2,2,2）

stem（n,f41）

subplot（2,2,3）

stem（n,imag（f71））

subplot（2,2,4）

stem（n,imag（f51））

figure

（2）

subplot（2,2,1）

stem（m,real（f72））

subplot（2,2,2）

stem（m,f42）

subplot（2,2,3）

stem（m,imag（f72））

subplot（2,2,4）

stem（m,imag（f52））

x9[n],x10[n]:

N=8;

n=0:

N-1;

x51=sin（0.125*pi*n）;

x41=cos（0.25*pi*n）;

x81=x41+j*x51;

m=0:

2*N-1;

x52=sin（0.125*pi*m）;

x42=cos（0.25*pi*m）;

x82=x42+j*x52;

f41=fft（x41,N）;

f51=fft（x51,N）;

f81=fft（x81,N）;

f83=f81;

fori=2:

N

f83（i）=conj（f81（N+2-i））;

end

f42=fft（x42,2*N）;

f52=fft（x52,2*N）;

f82=fft（x82,2*N）;

f84=f82;

fori=2:

2*N

f84（i）=conj（f82（2*N+2-i））;

end

figure

（1）

subplot（2,2,1）

stem（n,（f81+f83）/2）

subplot（2,2,2）

stem（n,f41）

subplot（2,2,3）

stem（n,abs（（f81-f83）*1i/2））

subplot（2,2,4）

stem（n,abs（f51））

figure

（2）

subplot（2,2,1）

stem（m,（f82+f84）/2）

subplot（2,2,2）

stem（m,f42）

subplot（2,2,3）

stem（m,abs（（f82-f84）*1i/2））

subplot（2,2,4）

stem（m,abs（f52））

x11[n]:

N=16;

n=0:

N-1;

x1=[1234000000000000];

x2=[1234123412341234];

x3=[1000200030004000];

f1=fft（x1,N）;

f2=fft（x2,N）;

f3=fft（x3,N）;

subplot（3,2,1）

stem（n,x1）

subplot（3,2,2）

stem（n,abs（f1））

subplot（3,2,3）

stem（n,x2）

subplot（3,2,4）

stem（n,abs（f2））

subplot（3,2,5）

stem（n,x3）

subplot（3,2,6）

stem（n,abs（f3））

4.思考题

（1）在N=8时，x2[n]和x3[n]的幅度谱会相同吗？

为什么？

N=16时，x2[n]和x3[n]的幅度谱会

相同吗？

为什么？

答：

不相同，由DFT的定义可知,x[n]的作用是对相同的K值的幅度值，所以虽然不影响频

谱分，布但是影响频率幅值.由于x2[n]和x3[n]的数值分布范围不同,所以导致在相同K

值时幅值不同，因此x2[n]和x3[n]的幅频特性不相同。

N=16时也不相同。

（2）比较x4[n]的8点和16点的DFT波形，这说明什么？

再比较x5[n]的8点和16点的DFT

波形，这说明了什么？

分析x4[n]和x5[n]的这些频谱异同产生的原因。

答：

x4[n]的8点和16点波形相同，x5[n]的8点和16点波形相同，说明仅当函数为偶函数

时，8点和16点DFT波形相同。

（3）比较x6（t）的64Hz采样时16，32，64点DFT波形，分析它们的异同产生的原因。

答：

波形相同，因为是对同一个数字频谱的采样。

（4）比较x6（t）的16Hz采样和64Hz采样时64点DFT波形，这说明什么？

答：

64Hz是16Hz扩张四倍后的波形，说明频率改变频谱的压缩或扩张。

（5）找出X7[k]与X4[k]=DFT[x4[n]]、X5[k]=DFT[x5[n]]的关系。

答：

X7[k]的实部和X4[k]相等，X7[k]的虚部和X5[k]的虚部相等。

（6）找出X8[k]与X4[k]=DFT[x4[n]]、X5[k]=DFT[x5[n]]的关系。

答：

X4[k]=（X8[k]+X8*[n-k]）/2,X5[k]=j/2（X8[k]-X8*[n-k]）。

（7）比较x9[n]、x10[n]和x11[n]的DFT结果，分析信号末尾补零、周期性延拓和时域插值在频

谱上的变化。

答：

信号末尾补零频谱不发生变化，周期性延拓频谱变化是周期性的，时域插值频谱压缩。

（8）如果周期信号的周期预先不知道，如何用FFT进行谱分析？

答：

周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较

结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差

要求就继续将截取长度加倍,重复比较,直到结果满足要求。

5.分析与结论

FFT变换即快速傅里叶变换的性质同DFT即离散傅里叶变换相同。

离散傅里叶变换有两个物理意义，一是，是对该序列的傅里叶变换w的抽样或者说对Z变换单位圆内的抽样。

二是，将该序列进行周期延拓后的傅里叶级数变换的主值序列。

用FFT作谱分析时，若给出的是连续信号，需要根据其最高频率确定采样速率，根据频率分辨率选取采样点数N，对其进行软件采样，产生对应序列。

若给出的是周期序列，则应该截取周期的整数倍进行谱分析，这样可以尽量减少实验造成的误差。

实验二IIR数字滤波器设计

1.实验目的

（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法。

（2）掌握数字滤波器的计算机仿真方法。

（3）通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。

2.实验原理

（1）IIR数字滤波器的系统函数可以写成

设计经典的IIR数字滤波器，就是要确定该系统函数中的2组系数：

Num:

b0,b1,...,bM

Den:

a1,a2,...,aN

也就是说，这两组系数确定了，IIR数字滤波器的理论设计就完成了。

接下去可以进行

结构变换和系数量化。

（2）IIR数字滤波器的设计可以借助模拟滤波器的设计方法。

先设计一个模拟滤波器，然后转

换成数字滤波器。

（3）模拟滤波器的原型有很多种。

常用的低通模拟滤波器原型有四种：

巴特沃斯（Butterworth）

滤波器、切比雪夫I型（ChebyshevⅠ）滤波器、切比雪夫II型（ChebyshevⅡ）滤波器、

椭圆（Elliptic）滤波器。

这些滤波器的增益（Gain）曲线如图2-1所示。

（4）滤波器设计指标可分为模拟指标和数字指标，绝对指标和相对指标。

一般来说，具体指

标有通带的边界频率、阻带的边界频率，通带最大衰减，阻带最小衰减四种。

低通滤波器

指标具体含义如图2-2所示。

（5）根据指标设计模拟滤波器，得到传输函数Hc（s）。

用双线性变换法，

其中，为了防止频响曲线发生畸变，设计指标的模拟频率与数字频率间的关系需做预畸变：

3.实验内容

（1）复习有关巴特沃斯、切比雪夫I型、切比雪夫II型、椭圆四种模拟滤波器设计和用双线

性变换法设计IIR数字滤波器的内容。

（2）人体心电图信号在测量过程中往往受到工业高频（主要是50Hz市电及其高次谐波）干扰，

所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。

x[n]为一实际心电图信号采样序列样本:

{x[n]}={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,

12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0-2-2,0,0,-2,-2,-2,-2,0}，

画出该信号的时域波形和频谱。

从频谱图中可以看到该信号存在高频干扰，噪声主要

集中在大于0.36πrad的频段。

（3）用双线性变换法设计低通IIR数字滤波器，滤除信号x[n]中的干扰成分。

根据待处理信号

的频谱，确定滤波器的设计指标。

例如：

通带[0,0.2π]内，最大衰减小于1dB,

阻带[0.3π,π]内，最小衰减大于15dB。

（4）利用Matlab函数设计IIR滤波器

a.先采用巴特沃兹原型设计。

根据设计指标，调用Matlab信号处理工具箱函数buttord（）

和butter（），得到以巴特沃兹模拟低通为原型的数字滤波器H（z）。

以0.02π为采样间隔，

画出数字滤波器在频率区间[0,π）上的频率响应特性曲线。

b.编写滤波器仿真程序，调用filter（）函数计算H（z）对心电图信号采样序列x[n]的响应y[n]。

画出滤波后心电图信号的时域波形和频谱。

c.调用Matlab的cheb1ord（），cheby1（），cheb2ord（），cheby2（），ellipord（），ellip（）等函数，

分别得到切比雪夫I型、切比雪夫II型、椭圆滤波器为模拟原型的数字滤波器。

重复

步骤（4）a，（4）b。

（5）利用Matlab的滤波器设计分析工具FDATool（FilterDesign&AnalysisTool）设计IIR滤波器。

a.在Matlab的“CommadWindow”里输入fdatool，打开滤波器设计分析工具。

b.在设计窗口中输入设计指标，如图2-3所示：

其中“Matchexactly”这项，如果选择“stopband”，则得到的滤波器阻带指标正好匹配，通带留有裕量；如果选择“passband”，则得到的滤波器通带指标正好匹配，阻带留有裕量。

c.点击DesignFilter，就可以得到IIR滤波器。

设计窗口中会显示设计得到的滤波器的频响

曲线。

d.在“CurrentFilterInformation”一栏里，可以看到滤波器的结构信息，IIR滤波器可以有很

多种结构，选择不同的结构，得到的滤波器系数形式也各不相同。

滤波器的结构可以通过

菜单“Edit”转换。

如同一个6阶直接II型转置型IIR滤波器，图2-4是单一节点结构（Single

Section），设计结果为2组系数，分子分母各7个系数；而图2-5是二阶节点级联结构

（Second-OrderSections，SOS），设计结果为3个二阶节点，每个节点分子分母各3个系数。

e.点击菜单“File”-“Export…”，可以把设计得到的滤波器系数导出到Matlab的工作空间（Workspace）或者存为文件。

f.使用设计得到的滤波器对心电图信号做滤波。

可参考步骤（4）b。

g.在“DesignMethod”的下拉菜单中选用不同的模拟滤波器原型，分别设计相应的IIR数字滤波器。

分别画出这些数字滤波器在频率区间［0,π）上的频率响应特性曲线。

①滤波前：

Wp=0.2*pi;

Ws=0.3*pi;

Rp=1;

Rs=15;

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

n=0:

length（x）-1;

[N,Wn]=buttord（Wp/pi,Ws/pi,Rp,Rs）;

[b,a]=butter（N,Wn）;

y=filter（b,a,x）;

figure

（1）;

freqz（b,a,50）;

figure

（2）

subplot（2,2,1）

stem（n,x,'.'）;

axis（[056-10050]）;

holdon;

xlabel（'\itn'）;

ylabel（'\itx[n]'）;

title（'Cardiogramsequence'）;

subplot（2,2,2）

stem（n,y,'.'）;

axis（[056-10050]）;

xlabel（'\itn'）;

ylabel（'\ity[n]'）;

title（'Cardiogramsequencefiltered'）;

m=0:

pi/1024:

pi/1024*1023;

subplot（2,2,3）

plot（m,abs（fft（x,1024）））;

axis（[0511-50500]）;

xlabel（'\it\omega'）;

ylabel（'\itX（e^j^\omega）'）;

title（'Spectrumofcardiogram'）;

subplot（2,2,4）

plot（m,abs（fft（y,1024）））;

axis（[0511-50500]）;

xlabel（'\it\omega'）;

ylabel（'\itY（e^j^\omega）'）;

title（'Spectrumofcardiogramfiltered'）;

②巴特沃兹滤波器

Wp=0.2*pi;

Ws=0.3*pi;

Rp=1;

Rs=15;

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

n=0:

55;

[N,Wn]=cheb1ord（Wp/pi,Ws/pi,Rp,Rs）;

[b,a]=cheby1（N,Rp,Wn）;

y=filter（b,a,x）;

figure

（1）;

freqz（b,a,50）;

figure

（2）

subplot（2,2,1）

stem（n,x,'.'）;

axis（[056-10050]）;

holdon;

xlabel（'n'）;

ylabel（'x[n]'）;

title（'Cardiogramsequence'）;

subplot（2,2,2）

stem（n,y,'.'）;

axis（[056-10050]）;

xlabel（'n'）;

ylabel（'y[n]'）;

title（'Cardiogramsequencefiltered'）;

m=0:

1023;

subplot（2,2,3）

plot（m,abs（fft（x,1024）））;

axis（[0511-50500]）;

xlabel（'w'）;

ylabel（'X（e^j^w）'）;

title（'Spectrumofcardiogram'）;

subplot（2,2,4）

plot（m,abs（fft（y,1024）））;

axis（[0511-50500]）;

xlabel（'w'）;

ylabel（'Y（e^j^w）'）;

title（'Spectrumofcardiogramfiltered'）

③切比雪夫Ⅰ滤波器

④切比雪夫Ⅱ滤波器

Wp=0.2*pi;

Ws=0.3*pi;

Rp=1;

Rs=15;

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

n=0:

55;

[N,Wn]=cheb2ord（Wp/pi,Ws/pi,Rp,Rs）;

[b,a]=cheby2（N,Rs,Wn）;

y=filter（b,a,x）;

figure

（1）;

freqz（b,a,50）;

figure

（2）

subplot（2,2,1）

stem（n,x,'.'）;

axis（[056-10050]）;

holdon;

xlabel（'n'）;

ylabel（'x[n]'）;

title（'Cardiogramsequence'）;

subplot（2,2,2）

stem（n,y,'.'）;

axis（[056-10050]）;

xlabel（'n'）;

ylabel（'y[n]'）;

title（'Cardiogramsequencefiltered'）;

m=0:

1023;

subplot（2,2,3）

plot（m,abs（fft（x,1024）））;

axis（[0511-50500]）;

xlabel（'w'）;

ylabel（'X（e^j^w）'