数字信号处理实验指导书-程序.doc

10数字信号处理实验指导书数字信号处理实验指导书课程组扬州大学信息工程学院2013年10月扬州大学信息工程学院2 信号采样与重建 - matlab仿真2.4 实验内容实验内容一MATLAB参考程序和仿真内容%*%f余弦信号的频率% M基2 FFT幂次数 N=2M为采样点数，这样取值是为了便于作基2的FFT分析%2. 采样频率Fs%*%function samples(f,Fs,M)N=2M; % fft点数=取样总点数Ts=1/Fs; % 取样时间间隔T=N*Ts; % 取样总时间=取样总点数*取样时间间隔n=0:N-1;t=n*Ts;Xn=cos(2*f*pi*t);subplot(2,1,1);stem(t,Xn);axis(0 T 1.1*min(Xn) 1.1*max(Xn);xlabel(t -);ylabel(Xn);Xk=abs(fft(Xn,N);subplot(2,1,2);stem(n,Xk);axis(0 N 1.1*min(Xk) 1.1*max(Xk);xlabel(frequency -);ylabel(!Xk!);%*%实验内容二 MATLAB参考程序和仿真内容% 信号取样与重构% k = 1 临界取样 0k1 欠取样display(Please input the value of k);k = input(k = );wm = 1; %信号带宽Ts = k*pi/wm; %取样间隔ws = 2*pi/Ts; %取样角频率wc = ws/2; %为简便，滤波器截止频率wc=ws/2n = -10:10; %定义取样点的数量m = fix(length(n)/2); %单边取样周期数,fix函数的功能为取整nTs = n*Ts; %计算每个取样点dt = 0.05;t = -m*Ts:dt:m*Ts;%Sa(t)信号取样f = sinc(nTs/pi); %计算Sa(t)在取样点处的函数值，得到取样信号%Sa(t)信号重构fa = f*Ts*wc/pi*sinc(wc/pi)*(ones(length(nTs),1)*t-nTs*ones(1,length(t); %Sa(t)包络线f2=sinc(t/pi);subplot(311);stem(nTs,f); %绘制取样信号hold onplot(t,f2,r:); %绘制包络线xlabel(t);title(对sa(t)信号进行取样);axis(-m*Ts m*Ts -0.5 1.2);hold offsubplot(312);h1=plot(t,fa); %绘制重构信号hold onfor i=-m*Ts:Ts:m*Ts ft = sinc(i/pi)*sinc(wc/pi*(t-i); h2=plot(t,ft,m:); %绘制重构信号的各个分量 hold on yy = sinc(i/pi); yy1 = yy:yy; plot(i,yy1,o);endxlabel(t);title(重构信号);axis(-m*Ts m*Ts -0.5 1.2);legend(h1,h2,重构信号,分信号);hold offsubplot(313);%计算原信号和重构信号之间的误差error=abs(fa-f2);plot(t,error);axis(-m*Ts m*Ts min(error) max(error)+0.1*max(error);xlabel(t);ylabel(error(t); %End4 离散傅里叶变换及其快速算法 matlab仿真4.4 实验内容MATLAB参考程序和仿真内容%*%functionx=ffts(mode,M)Nfft=2M;x=zeros(1,Nfft); %定义一个长度为Nfft的一维全0数组if mode=1 for n=0:Nfft-1 x(n+1)=sin(2*pi*n/Nfft); end end %定义一个长度为Nfft的单周期正弦序列if mode=2 for n=0:Nfft-1 x(n+1)=sin(4*pi*n/Nfft); end end %定义一个长度为Nfft的双周期正弦序列if mode=3 for n=0:Nfft/2-1 x(n+1)=sin(4*pi*n/Nfft); end end %定义一个长度为Nfft/2的正弦序列，后面一半为0序列。if mode=4 for n=0:Nfft-1 x(n+1)=square(2*pi*n/Nfft); end endif mode=5 for n=0:Nfft-1 x(n+1)=square(2*pi*n/Nfft); end endif mode=6 for n=0:Nfft/2-1 x(n+1)=square(4*pi*n/Nfft); end endn=0:Nfft-1;subplot(2,1,1);stem(n,x);axis(0 Nfft-1 1.1*min(x) 1.1*max(x);xlabel(Points-);ylabel(x(n);X=abs(fft(x,Nfft);subplot(2,1,2);stem(n,X);axis(0 Nfft-1 1.1*min(X) 1.1*max(X);xlabel(frequency-);ylabel(!X(k)!);%*%实验内容二MATLAB参考程序和仿真内容k=8; n1=0:1:19; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel(t/T);ylabel(x(n); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel(k);ylabel(X(k); n2=0:1:15; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel(t/T);ylabel(x(n); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel(k);ylabel(X(k);实验内容三 相关分析MATLAB参考程序和仿真内容x=1 3 -1 1 2 3 3 1; y=2 1 -1 1 2 0 -1 3; k=length(x); xk=fft(x,2*k);yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k); m=(-k+1):(k-1); stem(m,rm) xlabel(m); ylabel(幅度);MATLAB参考程序和仿真内容x=2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2; y=0 0 0 0 2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2; k=length(y); e=rand(1,k)-0.5; y=y+e; xk=fft(x,2*k); yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k); m=(-k+1):(k-1); stem(m,rm)xlabel(m); ylabel(幅度); 6 IIR滤波器的设计方法 matlab仿真6.4 实验内容实验内容一MATLAB参考程序和仿真内容%*%mode: 1-巴特沃斯低通；2-巴特沃斯高通；3-巴特沃斯带通；4-巴特沃斯带阻% 5-契比雪夫低通；6-契比雪夫高通；7-契比雪夫带通；8-契比雪夫带阻%fp1,fp2： 通带截止频率，当高通或低通时只有fp1有效%fs1, fs2： 阻带截止频率，当高通或低通时只有fs1有效%rp: 通带波纹系数%as: 阻带衰减系数%sample: 采样率%h: 返回设计好的滤波器系数%*%functionb,a=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)wp1=2*fp1/sample;wp2=2*fp2/sample;ws1=2*fs1/sample;ws2=2*fs2/sample;%得到巴特沃斯滤波器的最小阶数N和3bd频率wnif mode3N,wn=buttord(wp1,ws1,rp,as);elseif mode5N,wn=buttord(wp1 wp2,ws1 ws2,rp,as);%得到契比雪夫滤波器的最小阶数N和3bd频率wnelseif mode);phase=angle(freq_response);subplot(3,1,2);plot(f,phase);grid; %相频特性axis(0 sample/2 1.1*min(phase) 1.1*max(phase);ylabel(Phase);xlabel(Frequency-);h=impz(b,a,32); %32点的单位函数响应t=1:32;subplot(3,1,3);stem(t,h);grid;axis(0 32 1.2*min(h) 1.1*max(h);ylabel(h(n);xlabel(n-);%*%实验内容二1 数字低通滤波器设计MATLAB参考程序和仿真内容 B,A=butter(3,2*pi*1000,s); num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,s); num2,den2=bilinear(B,A,4000); h2,w=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),-.,f,abs(h2),-); grid; xlabel(频率/Hz ) ylabel(幅值/dB) 2 数字高通滤波器设计MATLAB参考程序和仿真内容 wc=2*1000*tan(2*pi*400/(2*1000); wt=2*1000*tan(2*pi*317/(2*1000); N,wn=cheb1ord(wc,wt,0.5,19,s); B,A=cheby1(N,0.5,wn,high,s); num,den=bilinear(B,A,1000); h,w=freqz(num,den); f=w/pi*500; plot(f,20*log10(abs(h); axis(0,500,-80,10); grid; xlabel() ylabel(幅度/dB) 3 数字带通滤波器设计MATLAB参考程序和仿真内容w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2,0 wr,3,10,s); B,A=butter(N,wn,s); num,den=bilinear(B,A,400); h,w=freqz(num,den); f=w/pi*200; plot(f,20*log10(abs(h); axis(40,160,-30,10); grid; xlabel(频率/kHz) ylabel(幅度/dB) 4 数字带阻滤波器设计MATLAB参考程序和仿真内容w1=2*1000*tan(2*pi*95/(2*1000); w2=2*1000*tan(2*pi*105/(2*1000); B,A=butter(1,w1, w2, stop,s); num,den=bilinear(B,A,1000); h,w=freqz(num,den); f=w/pi*500; plot(f,20*log10(abs(h); axis(50,150,-30,10); grid; xlabel(频率/kHz) ylabel(幅度/dB) 8 FIR滤波器的设计方法 - 窗函数法matlab仿真8.4 实验内容MATLAB参考程序和仿真内容%*%mode: 模式（1-高通；2-低通；3-带通；4-带阻）%n: 阶数，加窗的点数为阶数加1%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率%fs: 带通和带阻时指示上限频率%window:加窗（1-矩形窗；2-三角窗；3-巴特利特窗；4-汉明窗；% 5-汉宁窗；6-布莱克曼窗；7-凯泽窗；8-契比雪夫窗）%r: 代表加chebyshev窗的r值和加kaiser窗时的beta值%sample: 采样率%h: 返回设计好的FIR滤波器系数%*%mode: 模式（1-高通；2-低通；3-带通；4-带阻）%n: 阶数，加窗的点数为阶数加1%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率%fs: 带通和带阻时指示上限频率%window:加窗（1-矩形窗；2-三角窗； 3-巴特利特窗；4-汉明窗；% 5-汉宁窗；6-布莱克曼窗；7-凯泽窗； 8-契比雪夫窗）%r: 代表加chebyshev窗的r值和加kaiser窗时的beta值%sample: 采样率%h: 返回设计好的FIR滤波器系数%*%function h=usefir1(mode,n,fp,fs,window,r,sample)if window=1 w=boxcar(n+1);endif window=2 w=triang(n+1);endif window=3 w=bartlett(n+1);endif window=4 w=hamming(n+1);endif window=5 w=hanning(n+1);endif window=6 w=blackman(n+1);endif window=7 w=kaiser(n+1);endif window=8 w=chebwin(n+1);endwp=2*fp/sample;ws=2*fs/sample;if mode=1 h=fir1(n,wp,high,w);endif mode=2 h=fir1(n,wp,low,w);endif mode=3 h=fir1(n,wp,ws,w);endif mode=4 h=fir1(n,wp,ws,stop,w);end%set(gcf,menubar,menubar);m=0:n;subplot(3,1,1);stem(m,h);grid; %单位函数响应axis(0 n 1.1*min(h) 1.1*max(h);ylabel(h(n);xlabel(Points-);freq_response=freqz(h,1);Magnitude=20*log10(abs(freq_response);M=0:511;f=M*sample/(2*511);subplot(3,1,2);plot(f,Magnitude);grid; %幅频特性axis(0 sample/2 1.1*min(Magnitude) 1.1*max(Magnitude);ylabel(Magnitude);xlabel(Frequency-);phase=angle(freq_response);subplot(3,1,3);plot(f,phase);grid; %相频特性axis(0 sample/2 1.1*min(phase) 1.1*max(phase);ylabel(Phase);xlabel(Frequency-);%*%10 FIR滤波器的设计方法 - 频率采样法matlab仿真10.4 实验内容MATLAB参考程序和仿真内容%*%n： 阶数，加窗时的点数为阶数加1%fp： 高通和低通时指示截止频率，带通和带阻时指示下限频率%fs： 带通和带阻时指示上限频率% Window：加窗（1-矩形窗；2-三角窗；3-巴特利特窗；4-汉明窗；% 5-汉宁窗；6-布莱克曼窗；7-凯泽窗；8-契比雪夫窗）%r 代表加chebyshev窗的r值和加kaiser窗时的beta值%sample： 采样率%h： 返回设计好的FIR滤波器系数%*%function h=usefir2(n,f,m, window,r,sample)freq=2*f/sample;if window=1 w=boxcar(n+1);e