数字信号处理上机实验2.docx
《数字信号处理上机实验2.docx》由会员分享,可在线阅读,更多相关《数字信号处理上机实验2.docx(22页珍藏版)》请在冰豆网上搜索。
![数字信号处理上机实验2.docx](https://file1.bdocx.com/fileroot1/2023-4/20/b5bfb2f6-c416-477b-9b20-2b43fd9ca1a9/b5bfb2f6-c416-477b-9b20-2b43fd9ca1a91.gif)
数字信号处理上机实验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