数字信号处理上机实验2.docx
《数字信号处理上机实验2.docx》由会员分享，可在线阅读，更多相关《数字信号处理上机实验2.docx（22页珍藏版）》请在冰点文库上搜索。
数字信号处理上机实验2
数字信号处理实验报告
班级：
信息25
姓名：
赵梦然
学号：
2120502123
频率采样型滤波器
一．实验目的
1.学习使用频率采样型结构实现FIR滤波器，初步熟悉FIR滤
波器的线性相位特点。
2.直观体会频率采样型滤波器所具有的“滤波器组”特性，即
在并联结构的每条支路上可以分别得到输入信号的各次谐波。
3.学习如何使用周期冲激串检测所实现滤波器的频域响应。
2．实验内容
频率采样型滤波器是由一个梳状滤波器和若干路谐振器构成的，可用公式表
述如下：
其中r值理论上为1，实际中取非常接近1的值。
为了使系数为实数，可以将谐振器的共轭复根合并，不失一般性，假设N为
偶数，于是可以得到如图1所示的结构。
其中，
，
。
以下实验中假设频率采样型滤波器阶数
。
1.构造滤波器输入信号
，其中，
。
基波频率
，
，
，
，
，
，
，
，
。
设时域信号
的采样频率
，绘制出采样时刻从0到
的采样信号波形，其中采样点数为
，确认时域信号采样正确。
2.对采样信号的第二个周期
，进行离散傅里叶变换，画出幅频特性和相频特性图，观察并分析其特点。
3.设
，
，
，
，
，
，计算滤波器抽头系数
画出该滤波器的频谱图，观察并分析其幅频特性和相频特性。
4.编程实现图1所示的频率采样型滤波器结构，其中
取第3步中的值。
为了简化编程，梳状滤波器可以调用CombFilter.m，谐振器可以调用Resonator2.m，使用helpCombFilter和helpResonatoe2查看如何配置参数。
将第1步生成的采样信号通过该滤波器，画出输出信号第二个周期
的时域波形和频谱，并与第2步的频谱图进行对比，观察并分析二者的区别。
5.分别画出图1中前4路谐振器的输出信号第二个周期
的时域波形，观察并分析输出信号的特点。
6.将输入信号换成周期为
的冲激串，画出输出信号第二个周期
的幅频特性，并与第3步的滤波器幅频特性进行对比，观察并分析二者的关系。
7.思考并回答下列问题
（1）在第2步的幅频特性中，各次谐波的幅度与相应的时域信号幅度有什么关系？
（2）实验中为什么要观察第二个周期，如果直接观察第一个周期会怎么样？
（3）如果取
，观察会出现什么情况。
3．实验报告要求
1.按照实验内容编写Matlab程序，给出运行结果（图），并逐项进行分析讨
论。
2.提交完整的Matlab源程序。
3.回答思考题，撰写实验报告。
四．实验过程及结果分析
第一问：
MATLAB代码：
%%?
?
?
f0=50;
N=16;
fs=N*f0;
L=2*N;
T=1/fs;
t=0:
T:
（L-1）*T;
s=0.5*cos（2*pi*0*f0*t+0）+1*cos（2*pi*1*f0*t+pi/2）+0.5*cos（2*pi*2*f0*t+pi）+2*cos（2*pi*3*f0*t-pi/2）;
%%?
?
i=0:
1:
L-1;
stem（i,s（i+1）,'.'）;
title（'?
?
?
?
?
?
'）;
运行结果：
第二问：
MATLAB代码：
%%?
?
?
f0=50;
N=16;
fs=N*f0;
L=2*N;%%?
?
?
?
?
?
?
?
T=1/fs;
t=N*T:
T:
（L-1）*T;%%?
?
?
?
?
?
?
?
?
?
?
s=0.5*cos（2*pi*0*f0*t+0）+1*cos（2*pi*1*f0*t+pi/2）+0.5*cos（2*pi*2*f0*t+pi）+2*cos（2*pi*3*f0*t-pi/2）;
a=fft（s）;%%?
FFT
y=abs（a）;%%?
?
?
?
?
?
ang=angle（a）;%%?
?
?
?
?
?
i=0:
1:
N-1;
subplot（2,1,1）;%%subplot（m,n,p）,m?
?
?
?
?
?
m?
n?
?
?
?
?
n?
p?
?
?
?
?
?
?
?
.?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
stem（i,y（i+1）,'.'）;
title（'?
?
?
?
?
'）;
subplot（2,1,2）;
stem（i,ang（i+1）,'.'）;
title（'?
?
?
?
?
'）;
运行结果：
第三问：
MATLAB代码：
%%?
?
?
N=16;
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））];%%?
?
?
H（n）,?
?
?
?
?
?
?
?
h=ifft（H）;%%?
H（n）?
IFFT
h%%?
?
?
?
?
?
h
%%?
?
?
?
y=abs（H）;
ang=angle（H）;
i=0:
1:
N-1;
subplot（2,1,1）;
stem（i,y（i+1）,'.'）;%%stem?
?
?
?
i?
?
?
0?
?
?
?
?
?
?
?
?
?
?
1?
?
?
?
1,?
?
y（i+1）
title（'?
?
?
?
?
?
?
'）;
subplot（2,1,2）;
stem（i,ang（i+1）,'.'）;%stem?
?
?
?
i?
?
?
0?
?
?
?
?
?
?
?
?
?
?
1?
?
?
?
1,?
?
ang（i+1）
title（'?
?
?
?
?
?
?
'）;
运行结果：
输出的h：
h=
Columns1through6
0.0554-0.0000i0.0064-0.0000i-0.0548-0.0000i-0.0774-0.0000i-0.0286+0.0000i0.0841+0.0000i
Columns7through12
0.2143+0.0000i0.3006+0.0000i0.3006+0.0000i0.2143+0.0000i0.0841+0.0000i-0.0286+0.0000i
Columns13through16
-0.0774-0.0000i-0.0548-0.0000i0.0064-0.0000i0.0554-0.0000i
为了分析，将其做成连续的图形：
N=16;
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））];%%?
?
?
H（n）,?
?
?
?
?
?
?
?
h=ifft（H）;%%?
H（n）?
IFFT
h
%%?
?
?
?
?
?
?
?
?
w=0:
pi/200:
2*pi;
h1=zeros（1,length（w））;
ford=1:
length（w）
forn=1:
N
re=cos（w（d）*（n-1））;%%?
?
im=-sin（w（d）*（n-1））;%%?
?
e=complex（re,im）;%%?
?
?
?
.c=complex（a,b）?
?
c=a+bi
h1（d）=h1（d）+h（n）*e;
end
end
w=0:
pi/200:
2*pi;
subplot（2,1,1）;
plot（abs（h1））;%%plot?
?
?
?
title（'?
?
?
?
?
?
?
'）;
subplot（2,1,2）;
plot（angle（h1））;%%plot?
?
?
?
title（'?
?
?
?
?
?
?
'）;
运行结果：
第四问：
MATLAB代码：
%%?
?
?
f0=50;
N=16;
fs=N*f0;
L=2*N;
T=1/fs;
t=0:
T:
（L-1）*T;
s=0.5*cos（2*pi*0*f0*t+0）+1*cos（2*pi*1*f0*t+pi/2）+0.5*cos（2*pi*2*f0*t+pi）+2*cos（2*pi*3*f0*t-pi/2）;
x=s;%%CombFilter（x,N,r）?
?
?
?
?
?
?
x,?
?
?
s?
?
?
x
r=0.999;
y=CombFilter（x,N,r）;%%?
?
CombFilter（x,N,r）?
?
x=y;%%Resonator2（x,N,r,Order,H（i+1））?
?
?
?
?
?
?
x,?
?
?
y?
?
?
x
z=zeros（1,48）;%%?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
48?
?
?
?
?
?
?
?
48?
?
0?
?
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））;%%?
?
Resonator2（x,N,r,Order,H（i+1））?
?
z=z+y;
end
y=z/N;%%length（y）=48
fori=N:
1:
L-1;
y1（i-N+1）=y（i+1）;%%?
?
N?
L-1?
?
?
end
%%?
?
?
?
i=0:
1:
N-1;
subplot（3,1,1）;
stem（i,y1（i+1）,'.'）;
title（'?
?
?
?
'）;
%%?
?
y=fft（y1）;%%?
FFT
y1=abs（y）;%%?
?
?
ang=angle（y）;%%?
?
?
i=0:
1:
N-1;
subplot（3,1,2）;
stem（i,y1（i+1）,'.'）;%i?
?
?
0?
?
?
?
?
?
?
?
?
?
?
1?
?
?
?
1,?
?
y（i+1）
title（'?
?
?
?
'）;
subplot（3,1,3）;
stem（i,ang（i+1）,'.'）;%i?
?
?
0?
?
?
?
?
?
?
?
?
?
?
1?
?
?
?
1,?
?
ang（i+1）
title（'?
?
?
?
'）;
运行结果：
第五问：
MATLAB代码：
f0=50;
N=16;
fs=N*f0;
L=2*N;
T=1/fs;
t=0:
T:
（L-1）*T;
s0=0.5*cos（2*pi*0*f0*t+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）;
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））];
i=N:
1:
L-1;
y0=Resonator2（x,N,r,0,H
（1））;
subplot（2,2,1）,stem（i,y0（i+1）,'.'）;
title（'?
?
?
?
'）;
y1=Resonator2（x,N,r,1,H
（2））;
subplot（2,2,2）,stem（i,y1（i+1）,'.'）;
title（'?
?
?
?
'）;
y2=Resonator2（x,N,r,2,H（3））;
subplot（2,2,3）,stem（i,y2（i+1）,'.'）;
title（'?
?
?
?
'）;
y3=Resonator2（x,N,r,3,H（4））;
subplot（2,2,4）,stem（i,y3（i+1）,'.'）;
title（'?
?
?
?
'）;
运行结果：
第六问：
MATLAB代码：
clear
%%?
?
?
f0=50;
N=16;
fs=N*f0;
L=2*N;
s=[10000000000000001000000000000000];
%%?
?
?
?
?
?
?
?
?
?
?
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））;%%?
?
Resonator2?
?
z=z+y;
end
运行结果：
5．实验思考题：
1.在第2步的幅频特性中，各次谐波的幅度与相应的时域信号幅度有什么关系？
分析：
零频分量的幅度等于时域波形中直流的幅度，直流、二次、三次谐波的谱线的幅度为时域波形相应谐波对应幅度的一半。
因为正弦信号的频谱中有两根谱线：
和
，每一根的幅度为信号的一半（cos（t）=（
+
）/2）
2.实验中为什么要观察第二个周期，如果直接观察第一个周期会怎么样？
分析：
第二个周期信号输入进系统后不会发生失真，与系统的单位脉冲响应卷积后，输出信号开始变得周期，进入稳态。
直接观察第一个周期，信号还没有完全进入系统，卷积出的信号不够完整，输出的就不是一个完整的周期信号，有一段过渡段，故输出会发生失真。
3.如果取
，观察会出现什么情况。
分析：
当r=0.999时,第四题实验结果如下：
当r=0.95时,重做第四题实验结果如下：
4.如何理解第3步与第6步在工程使用中的区别？
分析：
两者均为低通滤波器，但是第三步是进行频率采样后得到H（k）后由内插函数得到传递函数H（z），从而得到系统频响H（ejw）,内插时用到了正弦函数。
而冲激串是一种理想化条件，冲激串的频响即为系统的频响，但是这在工程中时难以实现的，因为难以产生一个时间趋于零的信号。
六．实验总结
通过本次试验，深入理解了频率采样型滤波器的原理。
直观体会频率采样型滤波器所具有的“滤波器组”特性，即在并联结构的每条支路上可以分别得到输入信号的各次谐波，并学会了使用周期冲激串检测所实现滤波器的频域响应。
七．参考程序：
梳状滤波器函数CombFilter.m
functiony=CombFilter（x,N,r）
L=length（x）;
y=zeros（1,L+N）;
ifL>=N
y（1:
N）=x（1:
N）;
forn=N+1:
L
y（n）=x（n）-x（n-N）*r^N;
end
y（L+1:
L+N）=-x（L+1-N:
L）*r^N;
else
y（1:
L）=x（1:
L）;
forn=N+1:
N+L
y（n）=-x（n-N）*r^N;
end
end
谐振器函数Resonator2.m
functiony=Resonator2（x,N,r,Order,H）
ifOrder<0
display（'OrdershouldnotbenegativeinfunctionResonator2!
'）;
elseifOrder>N/2
display（'OrdershouldnotbegreaterthanN/2infunctionResonator2!
'）;
else
L=length（x）;
f=zeros（1,L）;
y=zeros（1,L）;
ifmod（N,2）==0
ifOrder==0
y
（1）=x
（1）*H;
forn=2:
L
y（n）=H*x（n）+y（n-1）*r;
end
elseifOrder==N/2
y
（1）=x
（1）*H;
forn=2:
L
y（n）=H*x（n）-y（n-1）*r;
end
else
w=exp（1i*2*pi*Order/N）;
alfa0=2*real（H）;
alfa1=-2*r*real（H.*w）;
beta1=2*r*cos（2*pi*Order/N）;
beta2=-r^2;
f
（1）=x
（1）;y
（1）=f
（1）*alfa0;
f
（2）=x
（2）+f
（1）*beta1;y
（2）=f
（2）*alfa0+f
（1）*alfa1;
forn=3:
L
f（n）=x（n）+f（n-1）*beta1-f（n-2）*r^2;
y（n）=f（n）*alfa0+f（n-1）*alfa1;
end
end
else
ifOrder==0
y
（1）=x
（1）*H;
forn=2:
L
y（n）=H*x（n）+y（n-1）*w*r;
end
else
w=exp（1i*2*pi*Order/N）;
alfa0=2*real（H）;
alfa1=-2*r*real（H.*w）;
beta1=2*r*cos（2*pi*Order/N）;
beta2=-r^2;
f
（1）=x
（1）;y
（1）=f
（1）*alfa0;
f
（2）=x
（2）+f
（1）*beta1;y
（2）=f
（2）*alfa0+f
（1）*alfa1;
forn=3:
L
f（n）=x（n）+f（n-1）*beta1-f（n-2）*r^2;
y（n）=f（n）*alfa0+f（n-1）*alfa1;
end
end
end
end
升级会员 