西安交通大学数字信号处理实验报告频率采样型滤波器.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后相频特性'）;
升级会员 