DSP实验报告-南邮.doc

上机实验内容:数组的加、减、乘、除和乘方运算。输入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。(A)clear;A=1 2 3 4;stem(A);A: (B)B=3,4,5,6;stem(B);B:(C)A=1 2 3 4;B=3,4,5,6;C=A+BC = 4 6 8 10A=1 2 3 4;B=3,4,5,6;C=A+BC(1)ans = 4 A=1 2 3 4;B=3,4,5,6;C=A+B;C(1:3)ans = 4 6 8 A=1 2 3 4;%stem(A);B=3,4,5,6;%stem(B);C=A+B;%C(1);%C(1:3);C(0)? Subscript indices must either be real positive integers or logicals.A=1 2 3 4;B=3,4,5,6;C=A+B;stem(C);A=1 2 3 4;B=3,4,5,6;C=A+B; n=0:1:3; n; stem(n,C);C = 4 6 8 10(D)A=1 2 3 4;B=3,4,5,6;D=A-B;stem(D);D:D = -2 -2 -2 -2(E)A=1 2 3 4;B=3,4,5,6;C=A+B;E=A.*B;stem(E);E:E=A.*BE = 3 8 15 24(F)A=1 2 3 4;B=3,4,5,6;F=A./B;stem(F); F: F=A./BF = 0.3333 0.5000 0.6000 0.6667(G)A=1 2 3 4;B=3,4,5,6;G=A.B;stem(G); G=A.BG = 1 16 243 4096 (2) 用MATLAB实现下列序列:clear;n=0:1:15;x1=0.8.n;n=0:1:15;n;stem(n,x1); n=0:1:15;a=(0.2+3*j)*n;x2=exp(a);n=0:1:15;n;figure;stem(n,x2); n=0:1:15;x3=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi);figure;stem(n,x3); d) 将(c)中的x(n)扩展为以16为周期的函数 , 绘出四个周期。n=0:1:63;x4=3*cos(0.125*pi*rem(n,16)+0.2*pi)+2*sin(0.25*pi*rem(n,16)+0.1*pi);figure;stem(n,x4);e) 将(c)中的x(n)扩展为以10为周期的函数，绘出四个周期。n=0:1:39;x5=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.25*pi*rem(n,10)+0.1*pi);figure;stem(n,x5); (3) 产生并绘出下列序列的样本: n=0:3; x=1 -1 3 5 x1=circshift(x,0,-2); x2=circshift(x,0,1); x3=2*x1-x2-2*x; stem(n,x3);x=1 -1 3 5;x1=zeros(1,4);x2=zeros(1,4);for k=1:5for n=0:3x3=circshift(x,0,-k)t=n.*x3;x1=x1+t;endx2=x2+x1;endstem(x2);(4) 绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注t=0:0.001:10;x=sin(2*pi*t);plot(t,x)title(x=sin(2*pi*t);xlabel(x);ylabel(t);t=0:0.01:10;x=sin(2*pi*t);plot(t,x);plot(t,x,r_);title(sin);xlabel(t);ylabel(x(t);title(sin);t=0:0.01:4;x1=cos(100*pi*t);x2=sin(pi*t);x=x1.*x2;plot(t,x)title(x=cos(100*pi*t)*sin(pi*t);xlabel(x);ylabel(t);(5) 编写函数 实现，绘出该函数的图形，起点为n1，终点为n2。function stepshift(n0,n1,n2); %单位阶跃序列,n0为时移量n=n1:n0-1; %n1、n2为序列的起止序列号nn=length(n);x=zeros(1,nn); %n0前信号赋值为0stem(n,x,fill) %绘出n1n0-1的波形(0值)hold onk=n0:n2;kk=length(k);x=ones(1,kk); %n0后信号赋值为1stem(k,x,fill) %绘出n1n0-1的波形(1值)hold offaxis(n1,n2,0,1.1)title(单位阶跃序列) 运行程序后，在Commond Windows中输入stepshift(6,-3,24)则得到平移后序列图形如下：(6) 给定一因果系统 求出并绘制H(z)的幅频响应与相频响应。clear all;k=64;b=1 sqrt(2);a=1 -0.67 0.9;w=0:pi/k:pi;h=freqz(b,a,w);subplot(2,1,1);plot(abs(h)grid onsubplot(2,1,2);plot(angle(h)grid on(7) 计算序列和序列的离散卷积，并作图表示卷积结果。A=8 -2 -1 2 3 B=2 3 -1 -1;C=conv(A ,B)plot(C) (8) 求以下差分方程所描述系统的单位脉冲响应 clear all;N=50;a=1 -2;b=1 0.1 -0.06;x1=1 zeros(1,N-1);n=0:1:N-1;h=filter(a,b,x1);stem(n,h)axis(-1 53 -2.5 1.2)三：思考题（1）对于有限长序列，如何用MATLAB计算其DTFT?fs=1000;t=0:1/fs:0.6; f1=100; f2=300; x=sin(2*pi*f1*t)+sin(2*pi*f2*t); subplot(711) plot(x); title(f1(100hz)f2(300hz)的正弦信号，初相0) xlabel(序列（n）) grid on number=512 y=fft(x,number); n=0:length(y)-1; f=fs*n/length(y); subplot(713) plot(f,abs(y); title(f1f2的正弦信号的fft（512点）) xlabel(频率hz) grid on x=x+randn(1,length(x); subplot(715) plot(x); title(原f1f2的正弦信号（含随机噪声）) xlabel(序列（n）) grid on y=fft(x,number); n=0:length(y)-1; f=fs*n/length(y); subplot(717) plot(f,abs(y); title(原f1f2的正弦信号（含随机噪声）的fft（512点）) xlabel(频率hz) grid on(3) 对于由两个子系统级联或并联的系统，如何用MATLAB计算它们的幅频响应与相频响应?答：级联转换为直接型：cas2dir 并联转换为直接型：par2dir 然后用freqz（）就行了。 实验二 快速傅里叶变换(FFT)及其应用 一、实验目的 (1) 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。 (2) 应用FFT对典型信号进行频谱分析。 (3) 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。 (4) 应用FFT实现序列的线性卷积和相关。 二、实验内容 实验中用到的信号序列: a) 高斯序列 b) 衰减正弦序列 c) 三角波序列d) 反三角波序列上机实验内容: 方法二：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高斯序列幅频特性);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高斯序列幅频特性);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);plot(abs(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);plot(abs(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);plot(abs(fft(xa);title( p=14高斯序列幅频特性); (2) n=0:1:15;A=1;f=0.0625;a=0.1;x=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,1);stem(n,x,.);axis(0 15 -2 2);title(f=0.0625衰减正弦序列时域特性);subplot(3,2,2);stem(fft(x),.);axis(0 15 -5 2);title(f=0.0625衰减正弦序列幅频特性); n=0:1:15;A=1;f=0.4375;a=0.1;x=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,3);stem(n,x,.);axis(0 15 -2 2);title(f=0.4375衰减正弦序列时域特性);subplot(3,2,4);stem(fft(x),.);axis(0 15 -5 2);title(f=0.4375衰减正弦序列幅频特性); n=0:1:15;A=1;f=0.5625;a=0.1;x=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,5);stem(n,x,.);axis(0 15 -2 2);title(f=0.5625衰减正弦序列时域特性);subplot(3,2,6);stem(fft(x),.);axis(0 15 -5 2);title(f=0.5625衰减正弦序列幅频特性); n=0:1:15;A=1;f=0.0625;a=0.1;x=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,1);stem(n,x,.);axis(0 15 -2 2);title(f=0.0625衰减正弦序列时域特性);subplot(3,2,2);plot(abs(fft(x);axis(0 20 0 5);title(f=0.0625衰减正弦序列幅频特性); n=0:1:15;A=1;f=0.4375;a=0.1;x=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,3);stem(n,x,.);axis(0 15 -2 2);title(f=0.4375衰减正弦序列时域特性);subplot(3,2,4);plot(abs(fft(x);axis(0 20 0 5);title(f=0.4375衰减正弦序列幅频特性); n=0:1:15;A=1;f=0.5625;a=0.1;x=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,5);stem(n,x,.);axis(0 15 -2 2);title(f=0.5625衰减正弦序列时域特性);subplot(3,2,6);plot(abs(fft(x);axis(0 20 0 5);title(f=0.5625衰减正弦序列幅频特性); (3)clear;for n=1:4;xc(n)=n-1;endfor n=5:8;xc(n)=8-(n-1);endm=0:7;subplot(2,3,1);stem(m,xc,.);title(三角波序列); Xc=fft(xc,8);k=0:1:7;subplot(2,3,2);stem(k,abs(Xc),.);title(三角波序列8点FFT);Xc=fft(xc,32);k=0:1:31;subplot(2,3,3);stem(k,abs(Xc),.);title(三角波系列32点FFT); clear;for n=1:4;xc(n)=4-(n-1);endfor n=5:8;xc(n)=(n-1)-4;endm=0:7;subplot(2,3,4);stem(m,xc,.);title(反三角波序列); Xc=fft(xc,8);k=0:1:7; subplot(2,3,5);stem(k,abs(Xc),.);title(反三角波序列8点FFT);Xc=fft(xc,32);k=0:1:31;subplot(2,3,6);stem(k,abs(Xc),.);title(反三角波序列32点FFT);(4) N=16,f=1/16时，其频谱：N=16;n=0:1:N-1;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/16)*n);subplot(2,1,1);plot(n,x,.);axis(0 15 -2 2);grid;title(连续时间信号);subplot(2,1,2);plot(abs(fft(x);axis(0 20 0 10);grid;title(连续时间信号频谱);N=16,f=1/64时，其频谱：N=16;n=0:1:N-1;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/64)*n);subplot(2,1,1);plot(n,x,.);axis(0 15 -2 2);grid;title(连续时间信号);subplot(2,1,2);plot(abs(fft(x);axis(0 20 0 10);grid;title(连续时间信号频谱);N=128,f=1/16时，其频谱：N=128,f=1/64时，其频谱：(5)循环卷积：n=0:15;p=8;q=2;xa=exp(-(n-p).2/q);subplot(2,3,1);stem(n,xa,.);title(xa波形);Xa=fft(xa,16);subplot(2,3,4);stem(abs(Xa),.);title(Xa(k)=FFTxa(n)的波形 ); n=0:1:15;A=1;f=0.0625;a=0.1;xb=exp(-a*n).*sin(2*pi*f*n);subplot(2,3,2);stem(n,xb,.);title(xb波形);Xb=fft(xb,16);subplot(2,3,5);stem(abs(Xb),.);title(Xb(k)=FFTxb(n)的波形 ); Y=Xa.*Xb;y=ifft(Y,16);subplot(2,3,3);stem(abs(Y),.);title(Y(k)=Xa(k)Xb(k)的波形); subplot(2,3,6);stem(y,.);title(y=IFFTY(k)的波形 );线性卷积：n=0:15;p=8;q=2;xa=exp(-(n-p).2/q);subplot(2,3,1);stem(n,xa,.);title(xa波形);Xa=fft(xa,16);subplot(2,3,4);stem(abs(Xa),.);title(Xa(k)=FFTxa(n)的波形 ); n=0:1:15;f=0.0625;a=0.1;xb=exp(-a*n).*sin(2*pi*f*n);subplot(2,3,2);stem(n,xb,.);title(xb波形);Xb=fft(xb,16);subplot(2,3,5);stem(abs(Xb),.);title(Xb(k)=FFTxb(n)的波形 ); N1=length(xa); N2=length(xb); N=N1+N2-1; xa=xa zeros(1,N2-1); xb=xb zeros(1,N1-1); n=1:N; k=n; Xa=fft(xa); Xb=fft(xb); Y=Xa.*Xb; subplot(2,3,3); stem(abs(Y),.); title(Y(k)=Xa(k).*Xb(k)的波形); y=ifft(Y(k);subplot(2,3,6);stem(y,.);title(y=ifftY(k)的波形 );六clear;xe=rand(1,512);for i=1:4 n(i)=i-1; xc(i)=n(i);endfor i=5:8 n(i)=i-1; xc(i)=8-n(i);end%重叠相加法yn=zeros(1,519);for j=0:7 xj=xe(64*j+1:64*(j+1); xak=fft(xj,71); xck=fft(xc,71); yn1=ifft(xak.*xck); temp=zeros(1,519); temp(64*j+1:64*j+71)=yn1; yn=yn+temp;end;n=0:518;figure(1)subplot(211);stem(n,yn);xlabel(n);ylabel(y(n);title(重叠相加法:xc(n)与xe(n)的线性卷积的时域波形);subplot(212);stem(n,abs(fft(yn);xlabel(k);ylabel(Y(k);axis(0,600,0,300);title(重叠相加法:xc(n)与xe(n)的线性卷积的幅频特性);%重叠保留法k=1:7;xe1=k-k;xe_1=xe1,xe;yn_1=zeros(1,519);for j=0:7 xj_1=xe_1(64*j+1:64*j+71); xak_1=fft(xj_1); xck_1=fft(xc,71); yn1_1=ifft(xak_1.*xck_1); temp=zeros(1,519); temp(64*j+1:64*j+64)=yn1_1(8:71); yn_1=yn_1+temp;end;n=0:518;figure(2)subplot(211);stem(n,yn_1);xlabel(n);ylabel(y(n);title( 重叠保留法:xc(n)与xe(n)的线性卷积的时域波形);subplot(212);stem(n,abs(fft(yn_1);xlabel(k);ylabel(Y(k);axis(0,600,0,300);title(重叠保留法:xc(n)与xe(n)的线性卷积的幅频特性);七clear;N=16;L=16;n=0:N-1;for i=1:N n(i)=i-1; xa(i)=exp(-(n(i)-8).*(n(i)-8)/2); xb(i)=exp(-0.1*n(i)*sin(2*pi*0.0625*n(i);endXa=fft(xa,2*N);Xb=fft(xb,2*N);rm=ifft(conj(Xa).*Xb); rm=rm(N+2:2*N) rm(1:N); m=(-N+1):(N-1);subplot(221)stem(m,rm);title(线性相关rxy);rm=ifft(conj(Xb).*Xa);rm=rm(N+2:2*N) rm(1:N);m=(-N+1):(N-1);subplot(222)stem(m,rm);title(线性相关ryx);Xa16=fft(xa,N);Xb16=fft(xb,N);rm=real(ifft(conj(Xa16).*Xb16);subplot(223)stem(rm);title(循环相关rxy);rm=real(ifft(conj(Xb16).*Xa16);subplot(224)stem(rm);title(循环相关ryx);八clear;N=16;L=16;n=0:N-1;for i=1:N n(i)=i-1; xa(i)=exp(-(n(i)-8).*(n(i)-8)/2); xb(i)=exp(-0.1*n(i)*sin(2*pi*0.0625*n(i);endXa=fft(xa,2*N);Xb=fft(xb,2*N);rm=real(ifft(conj(Xa).*Xa); rm=rm(N+2:2*N) rm(1:N); m=(-N+1):(N-1);subplot(211)stem(m,rm);title(xa自相关);rm=ifft(conj(Xb).*Xb);rm=rm(N+2:2*N) rm(1:N);m=(-N+1):(N-1);subplot(212)stem(m,rm);title(xb自相关);实验三 IIR数字滤波器的设计一、 实验目的1 掌握双线性变换法及脉冲响应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。2 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法设计的特点。3 熟悉巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器的频率特性。二、 实验内容：（1）fc=300;fr=200;Rb=0.8;At=20;T=0.001;wc=2*1000*tan(2*pi*300/(2*1000);wt=2*1000*tan(2*pi*200/(2*1000);N,wn=cheb1ord(wc,wt,0.8,20,s);B,A=cheby1(N,0.8,wn,high,s);num,dem=bilinear(B,A,1000);h,w=freqz(num,dem);f=w/pi*500;plot(f,20*log10(abs(h);axis(0,500,-300,10);grid;xlabel(频率 /HZ);ylabel(幅度 /dB);（2）双线性变换法：Wp=0.2*1000*pi;Ws=0.3*1000*pi;Rp=1;At=25;Fs=1000;Ts=1/Fs;wp=2*tan(Wp/2)/Ts;ws=2*tan(Ws/2)/Ts; N,wn=buttord(Wp,Ws,Rp,At,s); fprintf(滤波器阶数N=%.0fn,N);b,a=butter(N,wn,s);bz,az=bilinear(b,a,Fs);disp(分子系数ba); fprintf(%.4e,bz); fprintf(n); disp(分母系数aa); fprintf(%.4e,bz); fprintf(n);w=linspace(0,pi,1024);h=freqz(bz,az,w);plot(w/pi,20*log10(abs(h);axis(0 1 -350 50);grid;xlabel(频率 );ylabel(幅度(dB);脉冲响应不变法：wp=0.2*1000*pi;ws=0.3*1000*pi;Rp=1;At=25;Fs=1000;Ts=1/Fs;Wp=wp*Fs;Ws=ws*Fs; N,Wn=buttord(Wp,Ws,Rp,At,s);b,a=butter(N,Wn,s);bz,az=impinvar(b,a);w=linspace(0,pi,1024);h=freqz(bz,az,w);plot(w/pi,20*log10(abs(h);grid;xlabel(频率);ylabel(幅度(dB);（3）Wp=1.2*1000*pi;Ws=2*1000*pi;Rp=0.5;At=40;Fs=8000;Ts=1/Fs;wp=2*tan(Wp/2)/Ts;ws=2*tan(Ws/2)/Ts; N,wn=buttord(Wp,Ws,Rp,At,s); fprintf(滤波器阶数N=%.0fn,N);b,a=butter(N,wn,s);bz,az=bilinear(b,a,Fs);w=linspace(0,pi,1024);h=freqz(bz,az,w);plot(w/pi,20*log10(abs(h);axis(0 1 -350 50);grid;xlabel(频率/kHZ );ylabel(幅度(dB);巴特沃思滤波器阶数N=12clear;Wp=1.2*1000*pi;Ws=2*1000*pi;Rp=0.5;At=40;Fs=8000;Ts=1/Fs; N,Wn=cheb1ord(Wp,Ws,Rp,At,s); fprintf(滤波器阶数N=%.0fn,N);z0,p0,k0=cheb1ap(N,Rp);b0=k0*real(poly(z0);a0=real(poly(p0);H,W=freqs(b0,a0);plot(W*Wn,20*log10(abs(H)/max(abs(H),grid ;xlabel(频率/HZ);ylabel(幅度(dB);切比雪夫滤波器阶数N=6clear;Wp=1.2*1000*pi;Ws=2*1000*pi;Rp=0.5;At=40;Fs=8000;Ts=1/Fs; N,Wn=ellipord(Wp,Ws,Rp,At,s); fprintf(滤波器阶数N=%.0fn,N);b,a=ellip(N,Rp,At,Wn,s); omega=0:200:18000*pi;h=freqs(b,a,omega);gain=20*log10(abs(h)plot(omega/(2*pi),gain);axis(0 9000 -100 10);grid;xlabel(频率/HZ );ylabel(幅度(dB);椭圆滤波器阶数N=4（4）w1=2*30000*tan(2*pi*2000/(2*30000);w2=2*30000*tan(2*pi*3000/(2*30000);wr=2*30000*tan(2*pi*6000/(2*30000);Rp=3;At=5;Fs=30000;N,wn=buttord(w1 w2,1 wr,3,5,s); B,A=butter(N,wn,s);num1,den1=impinvar(B,A,Fs);h1,w=freqz(num1,den1); B,A=butter(N,wn,s);num2,den2=bilinear(B,A,Fs);h2,w=freqz(