数字信号处理实验二

Z.XYZ.XY数字信号处理实验二-频率采样型滤波器姓名：周翔宇班级：电信硕31学号：2140508028Z.XYZ.XY一实验目的1.1学习使用频率采样型结构实现FIR滤波器，初步熟悉FIR滤波器的线性相位特点。1.2直观体会频率采样型滤波器所具有的“滤波器组”特性，即在并联结构的每条支路上可以分别得到输入信号的各次谐波。1.3学习如何使用周期冲激串检测所实现滤波器的频域响应。二实验内容频率采样型滤波器是由一个梳状滤波器和若干路谐振器构成的，可用公式表述如下：1-01---N-1Nr-1NkkNNzrWkHzzH其中r值理论上为1，实际中取非常接近1的值。为了使系数为实数，可以将谐振器的共轭复根合并，不失一般性，假设N为偶数，于是可以得到如图1所示的结构。Z.XYZ.XY其中，]k[HRe20k，][H(k)WRe2-kN1rk。以下实验中假设频率采样型滤波器阶数16N。三．实验过程及结果分析1.构造滤波器输入信号30kkt(t)ss，其中，)2(cos0kkktkfAts。基波频率Hzf500，5.00A，11A，5.02A，23A，00，21，2，2-3。设时域信号ts的采样频率0Nffs，绘制出采样时刻从0到1-L的采样信号波形，其中采样点数为NL2，确认时域信号采样正确。Matlab源码：A=[0.510.52];Q=[0pi/2pi-pi/2];N=16;L=2*N;f0=50;fs=N*f0;T=1/fs;t=[0:T:(L-1)*T];s=zeros(1,length(t));fork=1:1:4s=s+A(k)*cos(2*pi*(k-1)*f0*t+Q(k));endk=0:L-1stem(k,s(k+1));运行结果：Z.XYZ.XY实验结果分析：信号由直流分量，一次谐波（基波），二次谐波，三次谐波组成；信号基波频率为50HZ，采样频率为fs=f0*N，采样点数为2*N；相当于对直流分量采样了32个幅度相等的点；对一次谐波（基波）采样了两个周期，每个周期16个点；对二次谐波采样了四个周期，每个周期8个点；对三次谐波采样了六个周期，每个周期6个点。2.对采样信号的第二个周期1-,...,1,LNNn，进行离散傅里叶变换，画出幅频特性和相频特性图，观察并分析其特点。Matlab源代码：A=[0.510.52];Q=[0pi/2pi-pi/2];N=16;L=2*N;f0=50;fs=N*f0;Z.XYZ.XYT=1/fs;t=[0:T:(L-1)*T];s=zeros(1,length(t));fork=1:1:4s=s+A(k)*cos(2*pi*(k-1)*f0*t+Q(k));endk=N:L-1;s_analys=s(k+1);y=fft(s_analys);show_module=abs(y);show_angle=angle(y);figure(1);stem(k,show_module(k-N+1));figure(2);stem(k,show_angle(k-N+1));结果：Z.XYZ.XY结果分析：采样信号由直流分量，一次谐波，二次谐波，三次谐波组成；对应频域有四个冲激，依次为直流分量，一次谐波，二次谐波，三次谐波的冲激；这与他们的幅度矩阵[0.510.52]组成，由于一次谐波，二次谐波，三次谐波到频域冲激变为1/2；所以频域四种信号对应得幅度比为1:1:0.5:2，这与幅频特性所显示的是一致的；相频特性由其初始相位[00.5pipi-0.5pi]得，它与图中幅频特性也是一致的。3.10H，NNjH1--exp1，NNjH1-2-exp2，013...43HHH，NNjH1-14-exp-14，N1-15-exp-15NjH，计算滤波器抽头系数1,-N,...,1,0,nnh画出该滤波器的频谱图，观察并分析其幅频特性和相频特性。Z.XYZ.XYMatlab源代码：N=16;H=zeros(1,16);H(1)=1;H(2)=exp(-1i*pi*(N-1)/N);H(3)=exp(-2i*pi*(N-1)/N);H(15)=exp(-14i*pi*(N-1)/N);H(16)=exp(-15i*pi*(N-1)/N);h=ifft(H)show_module_h=abs(H);show_angle_h=angle(H);figure(1);k=0:15;stem(k,show_module_h(k+1));figure(2);stem(k,show_angle_h(k+1));w=0:pi/200:2*pi;h1=zeros(1,length(w));ford=1:length(w)forn=1:Nre=cos(w(d)*(n-1));im=-sin(w(d)*(n-1));e=complex(re,im);h1(d)=h1(d)+h(n)*e;endendfigure(3);w=0:pi/200:2*pi;subplot(2,1,1),plot(abs(h1));xlabel('w');ylabel('·ùÖµ');title('·ùƵÌØÐÔ£¨×÷Á¬Ðøͼ£©');gridon;subplot(2,1,2),plot(angle(h1));xlabel('w');ylabel('·ùÖµ');title('ÏàƵÌØÐÔ£¨×÷Á¬Ðøͼ£©');gridon;运行结果：h=1至7列0.0625+0.0234i0.0625+0.0460i0.0625+0.0116i0.0625-0.0748i0.0625-0.1704i0.0625-0.2194i0.0625-0.1849iZ.XYZ.XY8至14列0.0625-0.0722i0.0625+0.0722i0.0625+0.1849i0.0625+0.2194i0.0625+0.1704i0.0625+0.0748i0.0625-0.0116i15至16列0.0625-0.0460i0.0625-0.0234iZ.XYZ.XY结果分析：由题目给出的H(K)利用ifft函数计算出抽头函数h(n)值；由该滤Z.XYZ.XY波器幅频特性看出此滤波器为低通窄带滤波器，只有处于低通的信号能通过；从图中观察的到高频被滤去了。4.编程实现图1所示的频率采样型滤波器结构，其中999.0r,kH取第3步中的值。为了简化编程，梳状滤波器可以调用CombFilter.m，谐振器可以调用Resonator2.m，使用helpCombFilter和helpResonatoe2查看如何配置参数。将第1步生成的采样信号通过该滤波器，画出输出信号第二个周期1-,...1,LNNn的时域波形和频谱，并与第2步的频谱图进行对比，观察并分析二者的区别。Matlab源代码:clear;A=[0.510.52];Q=[0pi/2pi-pi/2];N=16;L=2*N;f0=50;fs=N*f0;T=1/fs;t=0:T:(L-1)*Ts=zeros(1,length(t));fork=1:1:4s=s+A(k)*cos(2*pi*(k-1)*f0*t+Q(k));endH=zeros(1,16);H(1)=1;H(2)=exp(-1i*pi*(N-1)/N);H(3)=exp(-2i*pi*(N-1)/N);H(15)=exp(-14i*pi*(N-1)/N);H(16)=exp(-15i*pi*(N-1)/N);r=0.95;y=CombFilter(s,N,r);temp=y;Filter=zeros(1,48);fork=0:1:(N/2)temp_y=Resonator2(temp,N,r,k,H(k+1));Z.XYZ.XYFilter=Filter+temp_y;endFilter_temp=Filter/N;fori=N:1:L-1Filter_end(i-N+1)=Filter_temp(i+1);endi=0:N-1;figure(1);stem(i,Filter_end(i+1));title('时域输出信号');y_analys=fft(Filter_end);module_y=abs(y_analys);angle_y=angle(y_analys);figure(2);subplot(211);stem(i,module_y(i+1));title('幅频特性');subplot(212);stem(i,angle_y(i+1));title('相频特性');运行结果：Z.XYZ.XY结果分析：我们由图中可以看出幅频特性有五个冲激这与给出的H（k）是一致的；同时经由滤波器后输出信号发生变化；这说明四种信号有些被滤掉了；有幅频特性看出三次谐波被滤去了；最终输出的信号为直流，一次谐波，二次谐波的叠加。5.分别画出图1中前4路谐振器的输出信号第二个周期1-,...1,LNNn的时域波形，观察并分析输出信号的特点。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);Z.XYZ.XYs2=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);i=0;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;y=z/N;fori=N:1:L-1;y1(i-N+1)=y(i+1);endi=0:1:N-1;subplot(2,2,1),stem(i,y1(i+1));xlabel('n1');ylabel('振幅');title('时域波形');gridon;clearf0=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);Z.XYZ.XYx=y;z=zeros(1,48);i=1;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;y=z/N;fori=N:1:L-1;y1(i-N+1)=y(i+1);endi=0:1:N-1;subplot(2,2,2),stem(i,y1(i+1));xlabel('n2');ylabel('振幅');title('时域波形');gridon;clearf0=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+p