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

..数字信号处理实验报告姓名：______贺云天______学号：____2010035015____班级：_____电信硕23_____实验日期：2014年10月22日提交日期：2014年10月30日1实验二频率采样型滤波器一、实验目的1.学习使用频率采样型结构实现FIR滤波器，初步熟悉FIR滤波器的线性相位特点；2.直观体会频率采样型滤波器所具有的“滤波器组”特性，即在并联结构的每条支路上可以分别得到输入信号的各次谐波；3.学习使用周期冲激串检测所实现滤波器的频域响应。二、实验内容频率采样型滤波器是由一个梳状滤波器和若干路谐振器构成的，可用公式表述如下：𝐻(𝑧)=1−𝑟𝑁𝑧−𝑁𝑁∑𝐻(𝑘)1−𝑟𝑊𝑁−𝑘𝑧−1𝑁−1𝑘=0(1)其中r值理论上为1，实际中取非常接近1的值。为了使系数为实数，可以将谐振器的共轭复根合并，不失一般性，假设N为偶数，于是可以得到如图1所示的结构。图1N为偶数的实系数频率采样型结构滤波器2其中𝛼0𝑘=2𝑅𝑒[𝐻(𝑘)]，𝛼1𝑘=−2𝑟𝑅𝑒[𝐻(𝑘)𝑊𝑁𝑘]以下实验中假设频率采样型滤波器阶数N=16。1.构造滤波器输入信号𝑠(𝑡)=∑𝑠𝑘(𝑡)3𝑘=0，其中𝑠𝑘(𝑡)=𝐴𝑘cos⁡(2𝜋𝑘𝑓0𝑡+𝜙𝑘)，基波频率𝑓0=50𝐻𝑧，𝐴0=0.5，𝐴1=1，𝐴2=0.5，𝐴3=2，𝜙0=0，𝜙1=𝜋2，𝜙0=𝜋，𝜙3=−𝜋2。设时域信号s(t)的采样频率𝑓𝑠=𝑁𝑓0，绘制出采样时刻从0到L-1的采样信号波形，其中采样点数为𝐿=2𝑁，确认时域信号采样正确。2.对采样信号的第二个周期（𝑛=𝑁,𝑁+1,⋯,𝐿−1）进行离散傅里叶变换，画出幅频特性和相频特性图，观察并分析其特点。3.设𝐻(0)=1，𝐻(1)=e−𝑗𝜋(𝑁−1)𝑁，𝐻(2)=e−𝑗2𝜋(𝑁−1)𝑁，𝐻(3)=𝐻(4)=⋯=𝐻(13)=0，𝐻(14)=−e−𝑗14𝜋(𝑁−1)𝑁，𝐻(15)=−e−𝑗15𝜋(𝑁−1)𝑁，计算滤波器抽头系数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.构造滤波器输入信号𝑠(𝑡)=∑𝑠𝑘(𝑡)3𝑘=0，其中𝑠𝑘(𝑡)=𝐴𝑘cos⁡(2𝜋𝑘𝑓0𝑡+𝜙𝑘)，基波频率𝑓0=50𝐻𝑧，𝐴0=0.5，𝐴1=1，𝐴2=0.5，𝐴3=2，𝜙0=0，𝜙1=𝜋2，𝜙0=𝜋，𝜙3=−𝜋2。设时域信号s(t)的采样频率𝑓𝑠=𝑁𝑓0，绘制出采样时刻从0到L-1的采样信号波形，其中采样点数为𝐿=2𝑁，确认时域信号采样正确。3输入信号是直流、一次波、二次波、三次波的叠加，采样频率为16𝑓0，高于输入信号的2倍。输入信号的时域图像如上。2.对采样信号的第二个周期（𝑛=𝑁,𝑁+1,⋯,𝐿−1）进行离散傅里叶变换，画出幅频特性和相频特性图，观察并分析其特点。4由幅频特性图可知，在第0、1、2、3、13、14、15个点上有幅值，数值对称、相位反对称，其余点的幅度为0。第0点表示基波分量；第1点和第15点为直流分量，其幅值和是基波分量的2倍；第2点和第14点为二次分量，其幅值和与基波分量相等；第3点和第13点为三次分量，幅值和是基波分量的4倍。同时，从FFT结果可发现，幅值为0的分量相位并不相同，这是由计算机计算精度造成的。3.设𝐻(0)=1，𝐻(1)=e−𝑗𝜋(𝑁−1)𝑁，𝐻(2)=e−𝑗2𝜋(𝑁−1)𝑁，𝐻(3)=𝐻(4)=⋯=𝐻(13)=0，𝐻(14)=−e−𝑗14𝜋(𝑁−1)𝑁，𝐻(15)=−e−𝑗15𝜋(𝑁−1)𝑁，计算滤波器抽头系数h(n)，n=0,1,…,N-1，画出该滤波器的频谱图，观察并分析其幅频特性和相频特性。5滤波器抽头系数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步的频谱进行对比，观察并分析二者的区别。6第四题输出图像：信号通过滤波器后第二个周期的时域图像、FFT后的频谱信号直接作FFT后的频谱7幅值：与第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步的滤波器幅频特性进行对比，观察并分析二者的关系。8输入信号是周期为N的冲激串时，输出信号的时域图像和频谱滤波器幅频特性冲激函数在频域上为1，通过滤波器后被滤掉k=3~13区域，所以其幅值和滤波器本身的幅值相等，保留了直流、基波、二次分量。9五、思考题1.在第2步的幅频特性中，各次谐波的幅度与相应的时域信号幅度有什么关系？频率为0时的幅度等于时域波形中直流的幅度，直流、二次、三次谐波的谱线的幅度为时域波形相应谐波对应幅度的一半。因为正弦信号的频谱中有jwe和jwe两条谱线，幅度都为信号的一半，cos()=𝑗−𝑗2。2.实验中为什么要观察第二个周期，如果直接观察第一个周期会怎么样？第二个周期信号与系统的单位脉冲响应卷积后，输出信号开始变得周期，不会发生失真；直接观察第一个周期，信号没有完全进入系统，卷积结果不完整，输出会发生失真。3.如果取r=0.95，观察会出现什么情况。所有谐振器的极点从单位圆向内收缩一点，同时梳状滤波器的零点也移到r圆上，滤波器的幅频特性的幅度减小，通带内的仍然保持线性相位特性。104.如何理解第3步与第6步在工程使用中的区别？两者均为低通滤波器，但是第三步是进行频率采样后由内插函数得到传递函数H(z)，从而得到系统频响H(jwe)；而题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);11stem(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);hH_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:Ne=exp(-j*w(i)*(n-1));H2(i)=H2(i)+h(n)*e;endendw=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;%初始化数据12f0=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));