南邮数字信号处理实验报告(带问题答案小结).doc

南京邮电大学实 验 报 告实验名称 熟悉MATLAB环境 快速傅里叶变换(FFT)及其应用 IIR数字滤波器的设计 FIR数字滤波器的设计课程名称 数字信号处理A 班级学号_ 12006311_ 姓 名_张文欣_开课时间 2014/2015学年， 第 二 学期 实验一 熟悉MATLAB环境一、 实验目的（1）熟悉MATLAB的主要操作命令。（2）学会简单的矩阵输入和数据读写。（3）掌握简单的绘图命令。（4）用MATLAB编程并学会创建函数。（5）观察离散系统的频率响应。二、 实验内容(1) 数组的加、减、乘、除和乘方运算。输入A=1 2 3 4，B=3,4,5,6，求C=A+B， D=A-B，E=A.*B,F=A./B,G=A.B 。并用stem语句画出A、B、C、D、E、F、G。n = 0:1:3;A=1 2 3 4;subplot(4,2,1)stem(n,A)xlabel(n)ylabel(A)B=3,4,5,6;subplot(4,2,2)stem(n,B)xlabel(n)ylabel(B)C=A+B;subplot(4,2,3)stem(n,C)xlabel(n)ylabel(C)D=A-BSubplot(4,2,4)stem(n,D)xlabel(n)ylabel(D)E=A.*Bsubplot(4,2,5)stem(n,E)xlabel(n)ylabel(E)F=A./Bsubplot(4,2,6)stem(n,F)xlabel(n)ylabel(F)G=A.B subplot(4,2, 7)stem(n,G)xlabel(n)ylabel(G) (2) 用MATLAB实现下列序列：a) n=0:1:15;x1=0.8.n;stem(n,x1)xlabel(n)ylabel(x(n)title(2(a) b) n=0:1:15;i=sqrt(-1);a = 0.2+3*i;x2=exp(a*n);figuresubplot(1,2,1)stem(n,real(x2)xlabel(n)ylabel(x(n)实部)subplot(1,2,2)stem(n,imag(x2)xlabel(n)ylabel(x(n)虚部)c) n=0:1:15;x3=3*cos(0.125*pi*n+0.2*pi) + 2*sin(0.25*pi*n+0.1*pi);stem(n,x3)xlabel(n)ylabel(x(n)title(2(c)(4) 绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注：a) t=0:0.001:10;x=sin(2*pi*t);plot(t,x,r-)xlabel(t),ylabel(x(t),title(sin(2pit)b) t=0:0.001:4;x=cos(100*pi*t).*sin(pi*t); plot(t,x,b-)xlabel(t),ylabel(x(t),title(cos(100pit)*sin(pit)(6) 给定一因果系统，求出并绘制H(z)的幅频响应和相频响应。b=1,sqrt(2),1;a=1,-0.67,0.9;h,w=freqz(b,a);am=20*log10(abs(h);% am=20*log10(abs(h)为幅频响应取dBsubplot(2,1,1);plot(w,abs(h);xlabel(w);ylabel(幅频响应);title(系统响应)ph=angle(h);subplot(2,1,2); plot(w,ph);xlabel(w);ylabel(相频响应);(7) 计算序列8 -2 -1 2 3和序列2 3 -1 -3的离散卷积，并作图表示卷积结果。N=50;a=1 -2;b=1 0.1 -0.06;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y);xlabel(n);ylabel(幅度 );(8) 求以下差分方程所描述系统的单位脉冲响应h(n), N=50;a=1 -2;b=1 0.1 -0.06;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y);xlabel(n);ylabel(幅度 );每一小题均给出实验过程与结果（含实验程序、运行的数据结果和图形）；实验二 快速傅里叶变换(FFT)及其应用一、 实验目的(1) 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。(2) 应用FFT对典型信号进行频谱分析。(3) 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。(4) 应用FFT实现序列的线性卷积和相关。二、 实验内容实验中用到的信号序列a)高斯序列b)衰减正弦序列c)三角波序列d)反三角波序列(1) 观察高斯序列的时域和幅频特性，固定信号中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。n=0:15;p=8;q=2;xa=exp(-(n-p).2/q);subplot(3,2,1);stem(n,xa,.);title(q=2高斯序列时域特性);subplot(3,2,2);stem(fft(xa),.);title(q=2高斯序列幅频特性); n=0:15; p=8;q=4;xa=exp(-(n-p).2/q);subplot(3,2,3);stem(n,xa,.);title(q=4高斯序列时域特性);subplot(3,2,4);stem(fft(xa),.);title(q=4高斯序列幅频特性); n=0:15;p=8;q=8;xa=exp(-(n-p).2/q);subplot(3,2,5);stem(n,xa,.);title(q=8高斯序列时域特性);subplot(3,2,6);stem(fft(xa),.);title(q=8高斯序列幅频特性);由上图比较可知，当保持P不变时，随q的增大时域波形展宽且变得平滑，而频域波形且频域分量变少，不容易发生混叠，可见高斯序列中的q表示时域波形的陡峭程度。n=0:15;p=8;q=8;xa=exp(-(n-p).2/q);subplot(3,2,1);stem(n,xa,.);title(p=8高斯序列时域特性);subplot(3,2,2);stem(fft(xa),.);title(p=8高斯序列幅频特性); n=0:15;p=13;q=8;xa=exp(-(n-p).2/q);subplot(3,2,3);stem(n,xa,.);title(p=13高斯序列时域特性);subplot(3,2,4);stem(fft(xa),.);title(p=13高斯序列幅频特性); n=0:15;p=14;q=8;xa=exp(-(n-p).2/q);subplot(3,2,5);stem(n,xa,.);title(p=14高斯序列时域特性);subplot(3,2,6);stem(fft(xa),.);title(p=14高斯序列幅频特性);由上图比较可知当保持q不变时，随p的增大时域波形整体向右移动，由上图可见时域的窗口固定时在p=13,p=14有几点被泄露了，可见p值表示时域波形峰值的位置，p值越大泄露越多，而频域波形随p的增大频域分量会增多，容易产生混叠。(3) 观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列和的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。在和末尾补零，用N=32点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两种情况的FFT频谱还有相同之处吗？这些变化说明了什么？n=0:3;k=1:8;Xc(n+1) = n;Xc(n+5) =4-n;Xd(n+1) = 4-n;Xd(n+5) =n; %三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel(n);ylabel(时域特性);text(1,3,三角波);subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel(k);ylabel(幅频特性);text(4,10,三角波);%反三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel(n);ylabel(时域特性);text(3,3,反三角波);subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel(k);ylabel(幅频特性);text(4,10,反三角波);legend(N=8)由图可知，N=8时，正反三角的频域图形是相同的，因为做DFT时要先周期延拓，做完后取主值部分，而正反三角形周期延拓后是相同的，只差一个相位，因此得到的频域图形也是相同的 。n=0:3;k=1:32;Xc(n+1) = n;Xc(n+5) =4-n;Xd(n+1) = 4-n;Xd(n+5) =n; Xc(9:32)=0;Xd(9:32)=0;k=1:32;figure; %三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel(n);ylabel(时域特性);text(1,3,三角波);subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel(k);ylabel(幅频特性);text(4,10,三角波); %反三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel(n);ylabel(时域特性);text(3,3,反三角波);subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel(k);ylabel(幅频特性);text(4,10,反三角波);legend(N=32) N=32时，两者的频谱是不同的，因为此时再做周期延拓就不相同了。在后面补零对于正三角形在n=8时是连续的，而反三角波在n=8时有一个突变，时域中出现了陡峭的地方，在频域中频谱分量会增多。通过N=8 和N=32比较得，通过在原序列的末端补零，增加了采样的点数，是谱线增多，弱化了栅栏效应，但增多后的谱线形状是与时域信号的形状有关的。(5) 用FFT分别实现（p8，q2）和（a0.1，f0.0625）的16点循环卷积和线性卷积。n1=0:1:15;p1=8;q1=2;x=exp(-(n1-p1).2/q1);n2=0:1:15;a=0.1;f2=0.0625;y=(exp(-a*n2).*sin(2*pi*f2*n2);N=length(x);n=0:N-1;n3=0:30;X=fft(x);Y=fft(y);x32=x zeros(1,16);y32=y zeros(1,16);X32=fft(x32);Y32=fft(y32);z16=ifft(X.*Y);z32=ifft(X32.*Y32);subplot(2,2,1);plot(n,z16,-*);xlabel(n);ylabel(z(n);title(循环卷积结果);subplot(2,2,2);plot(n3,z32(1:2*N-1),-o);xlabel(n);ylabel(z(n);title(线性卷积结果);rm16=real(ifft(conj(X).*Y);rm32_0=real(ifft(conj(X32).*Y32);rm32=rm32_0(N+2:2*N) rm32_0(1:N);m=n;subplot(2,2,3);plot(m,rm16,-);xlabel(m);ylabel(rm);title(循环相关结果);m=-(N-1):N-1;subplot(2,2,4);plot(m,rm32,-o);xlabel(m);ylabel(rm);title(线性相关结果);每一小题均给出实验过程与结果（含实验程序、运行的数据结果和图形）；实验三 IIR数字滤波器的设计一、 实验目的(1) 掌握双线性变换法及脉冲响应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。 (2) 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。 (3) 熟悉巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器的频率特性。二、 实验内容（1）P162 例4.4设采样周期T=250(采样频率=4kHz)，分别用脉冲响应不变法和双线性变换法设计一个三阶巴特沃思低通滤波器，其3dB边界频率为=1kHz。%脉冲响应不变法fc=1000;fs=4000;B,A=butter(3, 2*pi*fc,s);num1,den1=impinvar(B,A,fs);h1,w=freqz(num1,den1);f = w/pi*fs/2;%双线性变换法fc=1000;fs=4000; B,A=butter(3, 2*fs*tan(pi*fc/fs),s);num2,den2=bilinear(B,A,fs);h2,w=freqz(num2,den2);f = w/pi*fs/2;plot(f,abs(h1),-.,f,abs(h2),-);grid onlegend(脉冲响应不变法,双线性变换法)xlabel(频率/Hz)ylabel(幅度)(2)=0.2kHz,=1dB,=0.3kHz,At=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。clear;fc = 200;Rp = 1;fr = 300;Rs = 25;T = 0.001;fs = 1/T;%脉冲响应不变法Wp=2*pi*fc;Ws=2*pi*fr;N,Wn = buttord(Wp, Ws, Rp, Rs,s)B,A=butter(N, Wn,s);num1,den1=impinvar(B,A,fs);h1,w=freqz(num1,den1);%双线性变换法Wp=2*fs*tan(pi*fc/fs);Ws=2*fs*tan(pi*fr/fs);N,Wn = buttord(Wp, Ws, Rp, Rs,s)B,A=butter(N, Wn,s);num2,den2=bilinear(B,A,fs);h2,w=freqz(num2,den2);f = w/pi*fs/2;plot(f,abs(h1),-.,f,abs(h2),-);grid onlegend(脉冲响应不变法,双线性变换法)xlabel(频率/Hz)ylabel(幅度)N = 9Wn = 1.3693e+003N = 6Wn = 1.7043e+003由图可以看出，利用脉冲响应不变法得到的滤波器在高频部分（阻带）衰减特性不好；而双线性变换法将模拟频率的无穷大映射到数字频率的pi,这样使得其对应的滤波器在高频部分衰减很快。(3)利用双线性变换法设计满足下列指标的巴特沃思数字低通滤波器，并作图验证设计结果：=1.2kHz,=2kHz,At,=8kHz。clear;fc =1200;Rp = 0.5;fr = 2000;Rs = 40;fs = 8000;Wp=2*fs*tan(pi*fc/fs);Ws=2*fs*tan(pi*fr/fs);N,Wn = buttord(Wp, Ws, Rp, Rs,s)B,A=butter(N, Wn,s);num2,den2=bilinear(B,A,fs);h2,w=freqz(num2,den2);f = w/pi*fs/2;plot(f,abs(h2),-);grid onlegend(双线性变换法)xlabel(频率/Hz)ylabel(幅度)N = 9Wn = 9.5918e+003每一小题均给出实验过程与结果（含实验程序、运行的数据结果和图形）；实验四 FIR数字滤波器的设计一、 实验目的(1) 掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉相应的计算机编程； (2) 熟悉线性相位FIR滤波器的幅频特性和相频特性； (3) 了解各种不同窗函数对滤波器性能的影响。二、 实验内容(1) N=45，计算并画出矩形窗、汉明窗、布莱克曼窗的归一化的幅度谱。N = 45;n=0:1:N-1;win=boxcar(N); %矩形窗 H,w=freqz(win,1); H = H/H(1); stem(n,win) figure plot(w,abs(H); N = 45;n=0:1:N-1;win=hamming(N); %汉明窗 H,w=freqz(win,1); H = H/H(1); stem(n,win) figure plot(w,abs(H); N = 45;n=0:1:N-1;win=blackman(N); %布莱克曼窗 H,w=freqz(win,1); H = H/H(1); stem(n,win) figure plot(w,abs(H); (2) 用矩形窗设计一个21阶的线性相位低通FIR数字滤波器，截止频率Wc0.25，求出滤波器系数，并绘出滤波器的幅频特性。修改程序，分别得到阶次为N41，61的滤波器，并显示其各自的幅频曲线。N=21矩形窗Wc=0.25*pi; N=21; M=(N-1)/2; %位移量 for n=0:(N-1) if (n= fix(M) %中间的点单独算 hd(n+1)=Wc/pi; else hd(n+1)=sin(Wc*(n-M) /(pi*(n-M);end;end; win=boxcar(N);%不同窗函数 h=hd.*win; H,w=freqz(h,1); n=0:1:N-1; subplot(3,1,1); stem(n,h) xlabel(n),ylabel(h(n)subplot(3,1,2); plot(w,abs(H);xlabel(omega),ylabel(幅度谱) subplot(3,1,3); plot(w,angle(H);xlabel(omega),ylabel(相位谱) legend(N=21 矩形窗)汉明窗Wc=0.25*pi; N=21; M=(N-1)/2; %位移量 for n=0:(N-1) if (n= fix(M) %中间的点单独算 hd(n+1)=Wc/pi; else hd(n+1)=sin(Wc*(n-M) /(pi*(n-M);end;end; win=hamming(N); %汉明窗 h=hd.*win; H,w=freqz(h,1); n=0:1:N-1; subplot(3,1,1); stem(n,h) xlabel(n),ylabel(h(n)subplot(3,1,2); plot(w,abs(H);xlabel(omega),ylabel(幅度谱) subplot(3,1,3); plot(w,angle(H);xlabel(omega),ylabel(相位谱) legend(N=21 汉明窗)N=41时矩形窗Wc=0.25*pi; N=41; M=(N-1)/2; %位移量 for n=0:(N-1) if (n= fix(M) %中间的点单独算 hd(n+1)=Wc/pi; else hd(n+1)=sin(Wc*(n-M) /(pi*(n-M);end;end; win=boxcar(N);%不同窗函数 h=hd.*win; H,w=freqz(h,1); n=0:1:N-1; subplot(3,1,1); stem(n,h) xlabel(n),ylabel(h(n)subplot(3,1,2); plot(w,abs(H);xlabel(omega),ylabel(幅度谱) subplot(3,1,3); plot(w,angle(H);xlabel(omega),ylabel(相位谱) legend(N=41 矩形窗)汉明窗Wc=0.25*pi; N=41; M=(N-1)/2; %位移量 for n=0:(N-1) if (n= fix(M) %中间的点单独算 hd(n+1)=Wc/pi; else hd(n+1)=sin(Wc*(n-M) /(pi*(n-M);end;end; win=hamming(N); %汉明窗 h=hd.*win; H,w=freqz(h,1); n=0:1:N-1; subplot(3,1,1); stem(n,h) xlabel(n),ylabel(h(n)subplot(3,1,2); plot(w,abs(H);xlabel(omega),ylabel(幅度谱) subplot(3,1,3); plot(w,angle(H);xlabel(omega),ylabel(相位谱) legend(N=41 汉明窗)N=61时矩形窗Wc=0.25*pi; N=61; M=(N-1)/2; %位移量 for n=0:(N-1) if (n= fix(M) %中间的点单独算 hd(n+1)=Wc/pi; else hd(n+1)=sin(Wc*(n-M) /(pi*(n-M);end;end; win=boxcar(N);%不同窗函数 h=hd.*win; H,w=freqz(h,1); n=0:1:N-1; subplot(3,1,1); stem(n,h) xlabel(n),