西安交通大学数字信号处理实验报告频率采样型滤波器.docx

西安交通大学数字信号处理实验报告频率采样型滤波器.docx

《西安交通大学数字信号处理实验报告频率采样型滤波器.docx》由会员分享，可在线阅读，更多相关《西安交通大学数字信号处理实验报告频率采样型滤波器.docx（24页珍藏版）》请在冰豆网上搜索。

西安交通大学数字信号处理实验报告频率采样型滤波器

西安交通大学数字信号处理实验报告-频率采样型滤波器

数字信号处理

实验报告

姓名：

______贺云天______

实验二频率采样型滤波器

一、实验目的

1.学习使用频率采样型结构实现FIR滤波器，初步熟悉FIR滤波器的线性相位特点；

2.直观体会频率采样型滤波器所具有的“滤波器组”特性，即在并联结构的每条支路上可以分别得到输入信号的各次谐波；

3.学习使用周期冲激串检测所实现滤波器的频域响应。

二、实验内容

频率采样型滤波器是由一个梳状滤波器和若干路谐振器构成的，可用公式表述如下：

（1）

其中r值理论上为1，实际中取非常接近1的值。

为了使系数为实数，可以将谐振器的共轭复根合并，不失一般性，假设N为偶数，于是可以得到如图1所示的结构。

图1N为偶数的实系数频率采样型结构滤波器

其中

，

以下实验中假设频率采样型滤波器阶数N=16。

1.构造滤波器输入信号

，其中

，基波频率

，

，

，

，

，

，

，

，

。

设时域信号s（t）的采样频率

，绘制出采样时刻从0到L-1的采样信号波形，其中采样点数为

，确认时域信号采样正确。

2.对采样信号的第二个周期（

）进行离散傅里叶变换，画出幅频特性和相频特性图，观察并分析其特点。

3.设

，

，

，

，

，

，计算滤波器抽头系数h（n），n=0,1,…,N-1，画出该滤波器的频谱图，观察并分析其幅频特性和相频特性。

4.编程实现图1所示的频率采样型滤波器结构，其中r=0.999，H（k）取第3步中的值。

为了简化编程，梳妆滤波器可以调用CombFilter.m，谐振器可以调用Resonator2.m，使用helpCombFilter和helpResonator2查看如何配置参数。

将第1步生成的采样信号通过该滤波器，画出输出信号第二个周期（n=N,N+1,…,L-1）的时域波形和频谱，并与第2步的频谱进行对比，观察并分析二者的区别。

5.分别画出图1中前4路谐振器的输出信号第二个周期（n=N,N+1,…,L-1）的时域波形，观察并分析输出信号的特点。

6.将输入信号换成周期为N的冲激串，画出输出信号第二个周期（n=N,N+1,…,L-1）的幅频特性，并与第3步的滤波器幅频特性进行对比，观察并分析二者的关系。

7.思考题

（1）在第2步的幅频特性中，各次谐波的幅度与相应的时域信号幅度有什么关系？

（2）实验中为什么要观察第二个周期，如果直接观察第一个周期会怎么样？

（3）如果取r=0.95，观察会出现什么情况。

三、实验要求

1.按照实验内容编写Matlab程序，给出运行结果（图），并逐项进行分析讨论。

2.提交完整的Matlab源程序。

3.回答思考题，撰写实验报告。

四、实验结果与分析

1.构造滤波器输入信号

，其中

，基波频率

，

，

，

，

，

，

，

，

。

设时域信号s（t）的采样频率

，绘制出采样时刻从0到L-1的采样信号波形，其中采样点数为

，确认时域信号采样正确。

输入信号是直流、一次波、二次波、三次波的叠加，采样频率为

，高于输入信号的2倍。

输入信号的时域图像如上。

2.对采样信号的第二个周期（

）进行离散傅里叶变换，画出幅频特性和相频特性图，观察并分析其特点。

由幅频特性图可知，在第0、1、2、3、13、14、15个点上有幅值，数值对称、相位反对称，其余点的幅度为0。

第0点表示基波分量；第1点和第15点为直流分量，其幅值和是基波分量的2倍；第2点和第14点为二次分量，其幅值和与基波分量相等；第3点和第13点为三次分量，幅值和是基波分量的4倍。

同时，从FFT结果可发现，幅值为0的分量相位并不相同，这是由计算机计算精度造成的。

3.设

，

，

，

，

，

，计算滤波器抽头系数h（n），n=0,1,…,N-1，画出该滤波器的频谱图，观察并分析其幅频特性和相频特性。

滤波器抽头系数h（n）可由H（k）作IFFT运算得出，第一行两幅图表示出H（k）的幅频特性和相频特性。

将H（k）内插后观察，发现H2为低通滤波器。

4.编程实现图1所示的频率采样型滤波器结构，其中r=0.999，H（k）取第3步中的值。

为了简化编程，梳妆滤波器可以调用CombFilter.m，谐振器可以调用Resonator2.m，使用helpCombFilter和helpResonator2查看如何配置参数。

将第1步生成的采样信号通过该滤波器，画出输出信号第二个周期（n=N,N+1,…,L-1）的时域波形和频谱，并与第2步的频谱进行对比，观察并分析二者的区别。

第四题输出图像：

信号通过滤波器后第二个周期的时域图像、FFT后的频谱

信号直接作FFT后的频谱

幅值：

与第2题对比，第0、1、2、14、15点不变，第3点和第13点的数值为0，这是因为三次分量在通过滤波器时被滤掉了，而滤波器在第0、1、2、14、15点的幅值为1，所以通过滤波器后的幅值没有发生变化。

相位：

直流分量的初始相位与滤波器H（0）相乘，得到新的相位；一次分量的相位由基波初始相位0与H

（1）相位相乘而得；其他分量的相位无法确定。

5.分别画出图1中前4路谐振器的输出信号第二个周期（n=N,N+1,…,L-1）的时域波形，观察并分析输出信号的特点。

H（0）输出波形为直流分量，所以波形的幅度不变；因为H（3）=0，所以H（3）输出波形也为0，幅度不变；其他两路分别为基波和二次分量的波形。

6.将输入信号换成周期为N的冲激串，画出输出信号第二个周期（n=N,N+1,…,L-1）的幅频特性，并与第3步的滤波器幅频特性进行对比，观察并分析二者的关系。

输入信号是周期为N的冲激串时，输出信号的时域图像和频谱

滤波器幅频特性

冲激函数在频域上为1，通过滤波器后被滤掉k=3~13区域，所以其幅值和滤波器本身的幅值相等，保留了直流、基波、二次分量。

五、思考题

1.在第2步的幅频特性中，各次谐波的幅度与相应的时域信号幅度有什么关系？

频率为0时的幅度等于时域波形中直流的幅度，直流、二次、三次谐波的谱线的幅度为时域波形相应谐波对应幅度的一半。

因为正弦信号的频谱中有

和

两条谱线，幅度都为信号的一半，

。

2.实验中为什么要观察第二个周期，如果直接观察第一个周期会怎么样？

第二个周期信号与系统的单位脉冲响应卷积后，输出信号开始变得周期，不会发生失真；直接观察第一个周期，信号没有完全进入系统，卷积结果不完整，输出会发生失真。

3.如果取r=0.95，观察会出现什么情况。

所有谐振器的极点从单位圆向内收缩一点，同时梳状滤波器的零点也移到r圆上，滤波器的幅频特性的幅度减小，通带内的仍然保持线性相位特性。

4.如何理解第3步与第6步在工程使用中的区别？

两者均为低通滤波器，但是第三步是进行频率采样后由内插函数得到传递函数H（z），从而得到系统频响H（

）；而题6中使用了冲激串，这是因为在理想条件下，冲激串的频响就是系统的频响，但是在实际生活中很难实现。

六、MATLAB程序代码

1.以下是MATLAB程序代码：

clc;clear;clf;

f0=50;

N=16;

fs=N*f0;%时域信号s（t）的采样频率

T=1/fs;%采样间隔

L=N*2;%采样点数

t=0:

T:

（L-1）*T;%构造输入信号

s0=0.5*cos（0）;

s1=1*cos（2*pi*f0*t+pi/2）;

s2=0.5*cos（2*pi*2*f0*t+pi）;

s3=2*cos（2*pi*3*f0*t-pi/2）;

s=s0+s1+s2+s3;

i=0:

1:

L-1;

stem（i,s（i+1））;

gridon;

2.以下是MATLAB程序代码：

clc;clear;clf;

f0=50;

N=16;

fs=N*f0;%时域信号s（t）的采样频率

T=1/fs;%采样间隔

L=N*2;%采样点数

t=N*T:

T:

（L-1）*T;

s0=0.5*cos（2*pi*0*f0*t+0）;

s1=1*cos（2*pi*f0*t+pi/2）;

s2=0.5*cos（2*pi*2*f0*t+pi）;

s3=2*cos（2*pi*3*f0*t-pi/2）;

s=s0+s1+s2+s3;

a=fft（s）;%FFT运算

a%输出a的各项值

y=abs（a）;

z=angle（a）;

i=0:

1:

N-1;

subplot（2,1,1）;

stem（i,y（i+1））;

subplot（2,1,2）;

stem（i,z（i+1））;

3.以下是MATLAB程序代码；

clear;clc;clf;

N=16;

H=[1,exp（-1i*pi*（N-1）/N）,exp（-j*2*pi*（N-1）/N）,0,0,0,0,0,0,0,0,0,0,0,-exp（-j*14*pi*（N-1）/N）,-exp（-j*15*pi*（N-1）/N）];

h=ifft（H,N）;

h

H_abs=abs（H）;

H_angle=angle（H）;

i=0:

1:

N-1;

subplot（2,2,1）;

stem（i,H_abs（i+1））;

title（'H幅频特性'）;

gridon;

subplot（2,2,2）;

stem（i,H_angle（i+1））;

title（'H相频特性'）;

gridon;

w=0:

pi/200:

2*pi;

H2=zeros（1,length（w））;

fori=1:

length（w）

forn=1:

N

e=exp（-j*w（i）*（n-1））;

H2（i）=H2（i）+h（n）*e;

end

end

w=0:

pi/200:

2*pi;

subplot（2,2,3）;

plot（abs（H2））;

title（'H2幅频特性'）;

gridon;

subplot（2,2,4）;

plot（angle（H2））;

title（'H2相频特性'）;

gridon;

4.以下是MATLAB程序代码；

clc;

%初始化数据

f0=50;

N=16;

fs=N*f0;

L=2*N;

T=1/fs;

t=0:

T:

（L-1）*T;

s0=0.5*cos（0）;

s1=1*cos（2*pi*1*f0*t+pi/2）;

s2=0.5*cos（2*pi*2*f0*t+pi）;

s3=2*cos（2*pi*3*f0*t-pi/2）;

s=s0+s1+s2+s3;

x=s;

r=0.999;

y=CombFilter（x,N,r）;

x=y;

z=zeros（1,48）;

fori=0:

1:

（N/2）

Order=i;

H=[1exp（-j*pi*（N-1）/N）exp（（-j*2*pi*（N-1）/N））00000000000-exp（（-j*14*pi*（N-1）/N））-exp（（-j*15*pi*（N-1）/N））];

y=Resonator2（x,N,r,Order,H（i+1））;

z=z+y;

end

y=z/N;

fori=N:

1:

L-1;

y1（i-N+1）=y（i+1）;%取出N到L-1的元素

end

%作时域图

i=0:

1:

N-1;

subplot（2,2,1）,stem（i,y1（i+1））;

xlabel（'n'）;

ylabel（'振幅'）;

title（'时域波形'）;

gridon;

%作频谱

y=fft（y1）;

y

y1=abs（y）;

ang=angle（y）;

i=0:

1:

N-1;

subplot（2,2,3）,stem（i,y1（i+1））;

xlabel（'n'）;

ylabel（'幅值'）;

title（'幅频特性'）;

gridon;

subplot（2,2,4）,stem（i,ang（i+1））;

xlabel（'n'）;

ylabel（'幅值'）;

title（'相频特性'）;

gridon;

5.以下是MATLAB程序代码；

clc;

f0=50;

N=16;

fs=N*f0;

T=1/fs;

L=N*2;

t=0:

T:

（L-1）*T

s0=0.5*cos（0）;

s1=1*cos（2*pi*f0*t+pi/2）;

s2=0.5*cos（2*pi*2*f0*t+pi）;

s3=2*cos（2*pi*3*f0*t-pi/2）;

s=s0+s1+s2+s3;%生成信号

%前4路谐振器的输出

r=0.999;

y1=CombFilter（s,N,r）;

H=[1,exp（-1i*pi*（N-1）/N）,exp（-1i*2*pi*（N-1）/N）,0,0,0,0,0,0,0,0,0,0,0,exp（-1i*14*pi*（N-1）/N）,exp（-1i*15*pi*（N-1）/N）];

y_h0=Resonator2（y1,N,r,0,H

（1））;

y_h0=y_h0/N;

y_h0=y_h0（N:

2*N-1）;

y_h1=Resonator2（y1,N,r,1,H

（2））;

y_h1=y_h1/N;

y_h1=y_h1（N:

2*N-1）;

y_h2=Resonator2（y1,N,r,2,H（3））;

y_h2=y_h2/N;

y_h2=y_h2（N:

2*N-1）;

y_h3=Resonator2（y1,N,r,3,H（4））;

y_h3=y_h3/N;

y_h3=y_h3（N:

2*N-1）;

%输出图像

k=0:

N-1;

subplot（2,2,1）;

stem（k,y_h0）;

title（'H（0）路输出信号第二个周期'）;

subplot（2,2,2）;

stem（k,y_h1）;

title（'H

（1）路输出信号第二个周期'）;

subplot（2,2,3）;

stem（k,y_h2）;

title（'H

（2）路输出信号第二个周期'）;

subplot（2,2,4）;

stem（k,y_h3）;

title（'H（3）路输出信号第二个周期'）;

6.以下是MATLAB程序代码；

clc;

N=16;

x=zeros（1,2*N）;

x

（1）=1;

x（N+1）=1;%冲激串

H=zeros（1,N）;

h=zeros（1,N）;

H=[1,exp（-j*pi*（N-1）/N）,exp（（-j*2*pi*（N-1）/N））,0,0,0,0,0,0,0,0,0,0,0,-exp（（-j*14*pi*（N-1）/N））,-exp（（-j*15*pi*（N-1）/N））];

r=0.999;

y=CombFilter（x,N,r）;%冲激串滤波后

y1=zeros;

fork=0:

N/2;

y1=y1+Resonator2（y,N,r,k,H（k+1））;

end

y0=y1/N;%经谐振器后

k=N:

2*N-1;

holdon;

subplot（2,2,1）;

stem（k,y0（k+1））;

title（'输出信号时域'）;

gridon;

Y=zeros（1,N）;

fork=N:

2*N-1;

Y（k-N+1）=y0（k+1）;%输出信号第二个周期

end

F=fft（Y）;

F

k=0:

1:

N-1;

holdon

subplot（2,2,2）

stem（k,abs（F（k+1）））

title（'FFT后幅频特性'）;

gridon;

subplot（2,2,3）

stem（k,angle（F（k+1）））

title（'FFT后相频特性'）;