信号处理实验报告.doc

数字信号处理课程设计报告 通信0903北方的狼课程设计所选题目：1 3 5 7 题。课程设计要求：使用MATLAB（或其它开发工具）编程实现上述内容，滤波器设计题目应尽量避免使用现成的工具箱函数。为便于分析与观察，设计中所有频谱显示中的频率参数均对应折叠频率归一化。课程设计内容：第一题： 1. 给定模拟信号：，式中，。对进行采样，可得采样序列 1） 选择采样频率=1 kHz，观测时间ms，观测所得序列及其幅频特性2） 改变采样频率=300Hz，观测此时的变化3） 令采样频率=200Hz，观测此时的变化要求分析说明原理，绘出相应的序列及其它们对应的幅频特性曲线，指出的变化，说明为什么？题目分析及回答： 分别以所要求采样频率对进行等间隔采样，得到，（ 为采样周期），可知，采样信号的频谱函数是在原模拟信号频谱函数的周期延拓；若以频率为自变量，则以采样频率为延拓周期；对频带限于的模拟信号，根据采样定理，只有当时，采样后才不会发生频谱混叠失真。但是在Matlab上是无法计算连续函数，只有在当足够大的时候，我们才将频谱混叠忽略不计，从而可对采样序列进行傅里叶变换；最后应用subplot()命令实现画图，并注意到要归一化。 由以下实验得到的图形结果可以看到：，当采样频率越大的时候，采样信号的频谱越陡峭，而其失真情况也越来越小，根据采样定理“采样频率越接近信号频率，其失真情况就越小”可以得到上述结论。 源代码及图形：（1） 采样频率fs=1000HZ：源代码：figure(1)a=50*sqrt(2)*pi;A=444.128;Tp=0.05;w=a;fs=1000;T=1/fs;n=0:(Tp*fs);t=T*n;x=A*exp(-a*t).*sin(w*t);X=fft(x);magy=abs(X);subplot(121);stem(n,x);xlabel(n);ylabel(x(n);title(x(n) 1000HZ);subplot(122);plot(2/Tp*t(1:Tp*fs/2),magy(1:Tp*fs/2);xlabel(omega);ylabel(|X(ejomega)|);title(幅频特性);波形如下：（2）采样频率fs=300HZ:源代码：figure(1)a=50*sqrt(2)*pi;A=444.128;Tp=0.05;w=a;fs=300;T=1/fs;n=0:(Tp*fs);t=T*n;x=A*exp(-a*t).*sin(w*t);X=fft(x);magy=abs(X);subplot(121);stem(n,x);xlabel(n);ylabel(x(n);title(x(n)300HZ);subplot(122);plot(2/Tp*t(1:Tp*fs/2),magy(1:Tp*fs/2);xlabel(omega);ylabel(|X(ejomega)|);title(幅频特性);波形如下：（3）采样频率fs为200HZ：源代码：figure(1)a=50*sqrt(2)*pi;A=444.128;Tp=0.05;w=a;fs=200;T=1/fs;n=0:(Tp*fs);t=T*n;x=A*exp(-a*t).*sin(w*t);X=fft(x);magy=abs(X);subplot(121);stem(n,x);xlabel(n);ylabel(x(n);title(x(n)200HZ);subplot(122);plot(2/Tp*t(1:Tp*fs/2),magy(1:Tp*fs/2);xlabel(omega);ylabel(|X(ejomega)|);title(幅频特性);（4）综合比较的代码：figure(1)a=50*sqrt(2)*pi;A=444.128;Tp=0.05;w=a;fs=1000;T=1/fs;n=0:(Tp*fs);t=T*n;x=A*exp(-a*t).*sin(w*t);X=fft(x);magy=abs(X);subplot(321);stem(n,x);xlabel(n);ylabel(x(n);title(1000HZ);subplot(322);plot(2/Tp*t(1:Tp*fs/2),magy(1:Tp*fs/2);xlabel(omega);ylabel(|X(ejomega)|);title(幅频特性);a=50*sqrt(2)*pi;A=444.128;Tp=0.05;w=a;fs=300;T=1/fs;n=0:(Tp*fs);t=T*n;x=A*exp(-a*t).*sin(w*t);X=fft(x);magy=abs(X);subplot(323);stem(n,x);xlabel(n);ylabel(x(n);title(300HZ);subplot(324);plot(2/Tp*t(1:Tp*fs/2),magy(1:Tp*fs/2);xlabel(omega);ylabel(|X(ejomega)|);title(幅频特性);a=50*sqrt(2)*pi;A=444.128;Tp=0.05;w=a;fs=200;T=1/fs;n=0:(Tp*fs);t=T*n;x=A*exp(-a*t).*sin(w*t);X=fft(x);magy=abs(X);subplot(325);stem(n,x);xlabel(n);ylabel(x(n);title(200HZ);subplot(326);plot(2/Tp*t(1:Tp*fs/2),magy(1:Tp*fs/2);xlabel(omega);ylabel(|X(ejomega)|);title(幅频特性);综合比较的波形：第三题：一个连续信号含两个频率分量，经采样得x(n)=sin2*0.125n+cos2*(0.125+f)n n=0,1,N-1已知N=16，f分别为1/16和1/64，观察其幅频特性；当N=128时，f不变，其结果有何不同，为什么？分析说明原因，并打印出相应的幅频特性曲线题目分析及回答：首先要在MATLAB中实现信号x(n)，然后分别按照进行等间隔采样；再调用函数FFT对它进行离散傅里叶变换；图画调用subplot（）来完成，并且要注意频率的归一化。由以下得到的波形可以看出：数字时域中每格1/N，频域中不失真时应该在0.125和0.125+f的地方有幅值，如当N=16, f =1/16时，两频率分量应该是0.125和0.125+1/16，对应到谱线上是2*(1/N)、3*(1/N)即第2和3条，分辨率足够把两谱线分开，但当f =1/64时由于两谱线靠的太近，1/N的分辨率不能分开他们，出现了栅栏效应。由于时域窗太小，所以有泄漏现象。当N=128时，分辨率等于 fs/N，f1/N时，可以分辨，不会发生栅栏效应，图中两谱线都能分清，且时域窗增大了，泄漏现象有所改善，可见提高分辨率可以增加信号的截取长度。源代码及波形：（1）n=16：figure(1) N=16;n=0:(N-1);df1=1/16;df2=1/64;x1=sin(2*pi*0.125*n)+cos(2*pi*(0.125*n+df1).*n);x2=sin(2*pi*0.125*n)+cos(2*pi*(0.125*n+df2).*n);H1=fft(x1);H2=fft(x2); subplot(221);stem(n,x1);xlabel(n);ylabel(x1(n);title(df=1/16); subplot(222);plot(2/N*n(1:N/2),abs(H1(1:N/2);xlabel(omega/pi);ylabel(|X(ejomega)|);title(df=1/16的频谱分析); subplot(223);stem(n,x2);xlabel(n);ylabel(x2(n);title(df=1/64);subplot(224);plot(2/N*n(1:N/2),abs(H2(1:N/2);xlabel(omega/pi);ylabel(|X(ejomega)|);title(df=1/64频谱分析);（2）n=128：figure(1)N=128;n=0:(N-1);df1=1/16;df2=1/64;x1=sin(2*pi*0.125*n)+cos(2*pi*(0.125*n+df1).*n);x2=sin(2*pi*0.125*n)+cos(2*pi*(0.125*n+df2).*n);H1=fft(x1);H2=fft(x2);subplot(221);stem(n,x1);xlabel(n);ylabel(x1(n);title(df=1/16);subplot(222);plot(2/N*n(1:N/2),abs(H1(1:N/2);xlabel(omega/pi);ylabel(|X(ejomega)|);title(df=1/16的频谱分析);subplot(223);stem(n,x2);xlabel(n);ylabel(x2(n);title(df=1/64);subplot(224);plot(2/N*n(1:N/2),abs(H2(1:N/2);xlabel(omega/pi);ylabel(|X(ejomega)|);title(df=1/64频谱分析);波形如下：第五题：一个序列为，使用FFT分析其频谱：1） 使用不同宽度的矩形窗截短该序列为M点长度，取M分别为： a) M=20 b) M=40 c)M=160 ；观察不同长度的窗对谱分析结果的影响；2）使用汉宁窗、哈明窗重做1）3）对三种窗的结果进行理论分析及比较。并绘出相应的幅频特性曲线题目分析及回答：首先在MATLAB编辑窗口中实现函数，然后分别用窗函数（矩形窗boxcar，汉宁窗hann,哈明窗hamming）进行窗口化，然后对窗口化所得到的序列进行快速的傅里叶（fft）变化，便可得到各个窗口化的频谱图形。源代码及图形：（1） 矩形窗（从左到右M依次为20,40，160）源代码：figure(1)M1=20n=0:M1-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(boxcar(M1);Y=fft(y);magy=abs(Y);subplot(1,3,1);stem(2/M1*n(1:M1/2),magy(1:M1/2);title(幅频特性)clear; M2=40n=0:M2-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(boxcar(M2);Y=fft(y);magy=abs(Y);subplot(1,3,2);stem(2/M2*n(1:M2/2),magy(1:M2/2);title(幅频特性)clear; M3=160n=0:M3-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(boxcar(M3);Y=fft(y);magy=abs(Y);subplot(1,3,3);stem(2/M3*n(1:M3/2),magy(1:M3/2);title(幅频特性)clear; （2） 哈明窗：（从左到右M依次为20,40，160） 源代码：figure(1)M1=20n=0:M1-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hamming(M1);Y=fft(y);magy=abs(Y);subplot(1,3,1);stem(2/M1*n(1:M1/2),magy(1:M1/2);title(幅频特性)clear;M2=40n=0:M2-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hamming(M2);Y=fft(y);magy=abs(Y);subplot(1,3,2);stem(2/M2*n(1:M2/2),magy(1:M2/2);title(幅频特性)clear;M3=160n=0:M3-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hamming(M3);Y=fft(y);magy=abs(Y);subplot(1,3,3);stem(2/M3*n(1:M3/2),magy(1:M3/2);title(幅频特性)clear;（3）汉宁窗：（从左到右M依次为20,40，160）源代码：figure(1)M1=20n=0:M1-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hann(M1);Y=fft(y);magy=abs(Y);subplot(1,3,1);stem(2/M1*n(1:M1/2),magy(1:M1/2);title(幅频特性)clear;M2=40n=0:M2-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hann(M2);Y=fft(y);magy=abs(Y);subplot(1,3,2);stem(2/M2*n(1:M2/2),magy(1:M2/2);title(幅频特性)clear;M3=160n=0:M3-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hann(M3);Y=fft(y);magy=abs(Y);subplot(1,3,3);stem(2/M3*n(1:M3/2),magy(1:M3/2);title(幅频特性)clear;（3） 综合比较：源代码：figure(1)M1=20n=0:M1-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(boxcar(M1);Y=fft(y);magy=abs(Y);subplot(3,3,1);stem(2/M1*n(1:M1/2),magy(1:M1/2);title(矩形窗 M=20)clear; M2=40n=0:M2-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(boxcar(M2);Y=fft(y);magy=abs(Y);subplot(3,3,2);stem(2/M2*n(1:M2/2),magy(1:M2/2);title(矩形窗 M=40)clear; M3=160n=0:M3-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(boxcar(M3);Y=fft(y);magy=abs(Y);subplot(3,3,3);stem(2/M3*n(1:M3/2),magy(1:M3/2);title(矩形窗 M=160)clear; M1=20n=0:M1-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hamming(M1);Y=fft(y);magy=abs(Y);subplot(3,3,4);stem(2/M1*n(1:M1/2),magy(1:M1/2);title(哈明窗 M=20)clear;M2=40n=0:M2-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hamming(M2);Y=fft(y);magy=abs(Y);subplot(3,3,5);stem(2/M2*n(1:M2/2),magy(1:M2/2);title(哈明窗 M=40)clear;M3=160n=0:M3-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hamming(M3);Y=fft(y);magy=abs(Y);subplot(3,3,6);stem(2/M3*n(1:M3/2),magy(1:M3/2);title(哈明窗 M=160)clear; M1=20n=0:M1-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hann(M1);Y=fft(y);magy=abs(Y);subplot(3,3,7);stem(2/M1*n(1:M1/2),magy(1:M1/2);title(汉宁窗 M=20)clear;M2=40n=0:M2-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hann(M2);Y=fft(y);magy=abs(Y);subplot(3,3,8);stem(2/M2*n(1:M2/2),magy(1:M2/2);title(汉宁窗 M=40)clear;M3=160n=0:M3-1;x=0.5*cos(11*pi*n/20)+cos(9*pi*n/20);y=x.*(hann(M3);Y=fft(y);magy=abs(Y);subplot(3,3,9);stem(2/M3*n(1:M3/2),magy(1:M3/2);title(汉宁窗 M=160)clear;第七题：分别利用矩形窗、汉宁窗、哈明窗设计一个N=11的线性相位FIR 低通和高通数字滤波器，截止频率，要求：求出各滤波器的单位脉冲响应h（n）；绘出各滤波器的幅频及相频响应曲线；观察各滤波器的通带波纹和阻带波纹；比较不同窗函数对滤波特性的影响。分析及回答：本题中的矩形窗，汉宁窗，哈明窗可以直接在MATLAB中的boxcar(N),hanning(N),hanming(N)来设计，再通过各滤波器的幅频及相频响应曲线，观察各滤波器的通带波纹和阻带波纹以及不同窗函数对滤波特性的影响。由下面的波形可以观察到不同的窗函数生成的过渡带宽度和阻带最小衰减是不同的，我们还可以根据这一点来选择窗函数。源代码及图形：（1） 矩形窗：源代码：figure(1)N=10;wn=1/3;b1lp=fir1(N,wn,boxcar(N+1);b1hp=fir1(N,wn,high,boxcar(N+1);M=1024;w=0:(M-1);t=0:10;h1lp=freqz(b1lp,1,w/M*2*pi);h1hp=freqz(b1hp,1,w/M*2*pi);subplot(321);stem(t,b1lp);xlabel(n)ylabel(h(n)title(低通：h(n)subplot(322);stem(t,b1hp);xlabel(n)ylabel(h(n)title(高通：h(n)subplot(323);plot(w(1:(M/2)/M*2,abs(h1lp(1:(M/2);title(低通：幅频特性)subplot(324);plot(w(1:(M/2)/M*2,abs(h1hp(1:(M/2);title(高通：幅频特性)subplot(325);plot(w/M*2,angle(h1lp);title(低通：相频特性)subplot(326);plot(w/M*2,angle(h1hp);title(高通：相频特性)clear;波形如下：（1） 哈明窗：源代码：figure(1)N=10;wn=1/3;b3lp=fir1(N,wn,hamming(N+1);b3hp=fir1(N,wn,high,hamming(N+1);M=1024;w=0:(M-1);t=0:10;h3lp=freqz(b3lp,1,w/M*2*pi);h3hp=freqz(b3hp,1,w/M*2*pi);subplot(321);stem(t,b3lp);xlabel(n)ylabel(h(n)title(低通：h(n)subplot(322);stem(t,b3hp);xlabel(n)ylabel(h(n)title(高通：h(n)subplot(323);plot(w(1:(M/2)/M*2,abs(h3lp(1:(M/2);title(低通：幅频特性)subplot(324);plot(w(1:(M/2)/M*2,abs(h3hp(1:(M/2);title(高通：幅频特性)subplot(325);plot(w/M*2,angle(h3lp);title(低通：相频特性)subplot(326);plot(w/M*2,angle(h3hp);title(高通：相频特性)clear;波形如下：（2） 汉宁窗：源代码：figure(1)N=10;wn=1/3;b2lp=fir1(N,wn,hanning(N+1);b2hp=fir1(N,wn,high,hanning(N+1);M=1024;w=0:(M-1);t=0:10;h2lp=freqz(b2lp,1,w/M*2*pi);h2hp=freqz(b2hp,1,w/M*2*pi);subplot(321);stem(t,b2lp);xlabel(n)ylabel(h(n)title(低通：h(n)subplot(322);stem(t,b2hp);xlabel(n)ylabel(h(n)title(高通：h(n)subplot(323);plot(w(1:(M/2)/M*2,abs(h2lp(1:(M/2);title(低通：幅频特性)subplot(324);plot(w(1:(M/2)/M*2,abs(h2hp(1:(M/2);title(高通：幅频特性)subplot(325);plot(w/M*2,angle(h2lp);title(低通：相频特性)subplot(326);plot(w/M*2,angle(h2hp);title(高通：相频特性)clear;波形如下：（3） 综合比较：源代码：figure(1)N=10;wn=1/3;b1lp=fir1(N,wn,boxcar(N+1);b1hp=fir1(N,wn,high,boxcar(N+1);b2lp=fir1(N,wn,hanning(N+1);b2hp=fir1(N,wn,high,hanning(N+1);b3lp=fir1(N,wn,hamming(N+1);b3hp=fir1(N,wn,high,hamming(N+1);t=0:10;figure(1)subplot(321);stem(t,b1lp);xlabel(n)ylabel(h(n)title(矩形窗低通)subplot(322);stem(t,b1hp);xlabel(n)ylabel(h(n)title(矩形窗高通)subplot(323);stem(t,b2lp);xlabel(n)ylabel(h(n)title(汉宁窗低通)subplot(324);stem(t,b2hp);xlabel(n)ylabel(h(n)title(汉宁窗高通)subplot(325);stem(t,b3lp);xlabel(n)ylabel(h(n)title(哈明窗低通)subplot(326);stem(t,b3hp);xlabel(n)ylabel(h(n)title(哈明窗高通)M=1024;w=0:(M-1);h1lp=freqz(b1lp,1,w/M*2*pi);h1hp=freqz(b1hp,1,w/M*2*pi);h2lp=freqz(b2lp,1,w/M*2*pi);h2hp=freqz(b2hp,1,w/M*2*pi);h3lp=freqz(b3lp,1,w/M*2*pi);h3hp=freqz(b3hp,1,w/M*2*pi);