《数字信号处理胡广书例题作业程序.doc》由会员分享，可在线阅读，更多相关《数字信号处理胡广书例题作业程序.doc（48页珍藏版）》请在三一文库上搜索。
1、1、%-filter求卷积，B(Z)/A(Z)=H(Z),已知B(Z)和A(Z)，求y(n)=x(n)*h(n)-clear;x=ones(100);t=1:100;b=.001836,.007344,.011016,.007374,.001836;a=1,-3.0544,3.8291,-2.2925,.55075;% y=filter(b,a,x);% 求所给系统的输出，本例实际上是求所给系统的阶跃响应；plot(t,x,r.,t,y,k-);grid on;ylabel(x(n) and y(n)xlabel(n)1、%-filter求卷积，B(Z)/A(Z)=H(Z),已知B(Z)和A(
2、Z)，求y(n)=x(n)*h(n)-clear;x=ones(100);t=1:100;b=.001836,.007344,.011016,.007374,.001836;a=1,-3.0544,3.8291,-2.2925,.55075;% y=filter(b,a,x);% 求所给系统的输出，本例实际上是求所给系统的阶跃响应；plot(t,x,r.,t,y,k-);grid on;ylabel(x(n) and y(n)xlabel(n)第一章 产生信号，求卷积和自相关函数1、%信号产生n=0:100;%工频f0=50;A=220;fs=400;x1=A*sin(2*pi*f0*n/fs
3、);subplot(321);plot(n,x1);xlabel(n);ylabel(x1(n) ;grid on;%率减正弦f0=2;A=2;alf=0.5;fs=16;x2=A*exp(-alf*n/fs).*sin(2*pi*f0*n/fs);subplot(323);plot(n,x2);xlabel(n);ylabel(x2(n) ;grid on;%谐波信号f0=5;A1=1.0;A2=0.5;A3=0.2;fs=100;x3=A1*sin(2*pi*f0*n/fs)+A2*sin(2*pi*2*f0*n/fs)+A3*sin(2*pi*3*f0*n/fs);subplot(322
4、);plot(n,x3);xlabel(n);ylabel(x3(n) ;grid on;%哈明窗f0=10;fs=1000;x4=0.54-0.46*cos(2*pi*f0*n/fs);subplot(324);plot(n,x4);xlabel(n);ylabel(x4(n) ;grid on;%采样n=-50:50;f0=10;fs=400;w=2*pi*f0*n/fs;x5=sinc(w);subplot(325);plot(n,x5);xlabel(n);ylabel(x5(n) ;grid on;2、% 产生均匀分布的白噪信号，使均值为0，功率为p%-clear;p=0.01;N=
5、50000;u=rand(1,N);u=u-mean(u);a=sqrt(12*p); u1=u*a;power_u1=dot(u1,u1)/N subplot(211)plot(u1(1:200);grid on;ylabel(u(n) xlabel(n) 3、% 产生高斯分布的白噪信号，使功率为p，并观察数据分布的直方图%-clear;p=0.1;N=500000;u=randn(1,N);a=sqrt(p);u=u*a; power_u=var(u);subplot(211)plot(u(1:200);grid on;ylabel(u(n);xlabel(n) subplot(212)h
6、ist(u,50);grid on; ylabel(histogram of u(n); 4、% 产生一 sinc 函数；%-clear;n=200;stept=4*pi/n;t=-2*pi:stept:2*pi;y=sinc(t);plot(t,y,t,zeros(size(t);ylabel(sinc(t);xlabel(t=-2*pi2*pi);grid on;5、% 产生一 chirp 信号；% chirp(T0,F0,T1,F1):% T0: 信号的开始时间； F0:信号在T0时的瞬时频率，单位为Hz；% T1: 信号的结束时间； F1:信号在T1时的瞬时频率，单位为Hz；%-cle
7、ar; t=0:0.001:1; x=chirp(t,0,1,125);plot(t,x);ylabel(x(t) xlabel(t) 6、% 计算两个序列的线性卷积；%-clear;N=5;%第一个序列的长度M=6;%第二个序列的长度L=N+M-1;x=1,2,3,4,5;h=6,2,3,6,4,2;y=conv(x,h);nx=0:N-1;nh=0:M-1;ny=0:L-1;subplot(231);%绘制xstem(nx,x,.k);xlabel(n);ylabel(x(n);grid on;subplot(232);%绘制hstem(nh,h,.k);xlabel(n);ylabel(
8、h(n);grid on;subplot(233);%绘制卷积stem(ny,y,.k);xlabel(n);ylabel(y(n);grid on; 7、% 求两个序列的互相关函数，或一个序列的自相关函数；%-clear;N=500;p1=1;p2=0.1;f=1/8;Mlag=50;%自相关的单边长度u=randn(1,N);n=0:N-1;s=sin(2*pi*f*n);% 混有高斯白噪的正弦信号的自相关u1=u*sqrt(p1);%高斯白噪声x1=u1(1:N)+s;%混合信号rx1=xcorr(x1,Mlag,biased);%自相关，无偏估计subplot(221);plot(x1
9、(1:Mlag);title(信号 x1);xlabel(n);ylabel(x1(n);grid on;subplot(223);plot(-Mlag:Mlag),rx1);title(x1自相关);grid on;xlabel(m);ylabel(rx1(m);% 高斯白噪功率由原来的p1减少为p2，再观察混合信号的自相关u2=u*sqrt(p2);%改变高斯白噪声x2=u2(1:N)+s;%新的混合信号rx2=xcorr(x2,Mlag,biased);subplot(222);plot(x2(1:Mlag);title(信号 x2);xlabel(n);ylabel(x2(n);gri
10、d on;subplot(224);plot(-Mlag:Mlag),rx2);title(x2自相关);grid on;xlabel(m);ylabel(rx2(m); 8、%求序列的自相关函数clearN=500;Mlag=50;%单边长度nx=0:N-1;x=exp(-nx*0.1);rx=xcorr(x,Mlag,biased);nrx=-Mlag:Mlag;%自相关序列的程度subplot(211);plot(nx,x);xlabel(n);ylabel(x(n) ;grid on;subplot(212);plot(nrx,rx);xlabel(n);ylabel(rx(n) ;g
11、rid on;9、%正弦加白噪声，自相关p=0.1;N=5000;Mlag=100;u=rand(1,N);u=u-mean(u);a=sqrt(12*p);u=u*a;power_u=dot(u,u)/Nnx=1:1000;x=1.414*sin(nx*pi/16.0);x1=x(1:1000)+5*u(1:1000);rx=xcorr(x1,Mlag,biased);nrx=-Mlag:Mlag;subplot(211);plot(x1(1:200);xlabel(n);ylabel(x(n) ;grid on;subplot(212);plot(nrx,rx);xlabel(n);yla
12、bel(rx(n) ;grid on;1、%产生信号，求卷积，FFT，求平均clear all; N=1024;%采样点数 fs=100.0;%采样频率 alf1=-1.0; f1=5; alf2=-1.5; f2=8; alf3=-0.7; f3=10; %产生x和w两个信号 %产生x u=rand(1,N); u=u-mean(u);%均值为0的白噪声 t=0:1/fs:(N-1)/fs; x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3*t); x=
13、x+u; x=x/max(x); %产生w alf4=-1.0; w=1.0*exp(alf4*t); %x=x-mean(x); figure(1); subplot(211); plot(t,x,t,w,r);title(x(t) w(t) % 应用FFT 求频谱； f=0:fs/N:fs/N*(N-1); X=fft(x,N); X=abs(X);% X=20*log10(X); subplot(212); plot(f(1:N/2),X(1:N/2);title(x(t)频谱) % stem(f,X,.);grid on; xlabel(Hz) %求卷积-y=x.*w;%时域相乘，频域
14、卷积figure(2);subplot(211);plot(t,y);title(正弦加白噪声后与w时域相乘)Y=fft(y,N); Y=abs(Y);subplot(212);plot(f(1:N/2),Y(1:N/2);title(正弦加白噪声后与w时域相乘的FFT)xlabel(Hz)xlabel(Hz)%平均1000次- figure(3); Y=zeros(1,N); u=rand(1,1000*N);u=u-mean(u);%零均值白噪声 for i=0:999 x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*
15、f2*t)+1.0*exp(alf3*t).*sin(2*pi*f3*t); x=x+10*u(1+i*N:i*N+N); x=x/max(x); X=fft(x,N); X=abs(X); Y=Y+X; end Y=Y/1000; subplot(211);plot(f(1:N/2),Y(1:N/2);title(正弦加白噪声的FFT 1000次平均) xlabel(Hz) Y=zeros(1,N); for i=0:999 x=1.0*exp(alf1*t).*sin(2*pi*f1*t)+1.0*exp(alf2*t).*sin(2*pi*f2*t)+1.0*exp(alf3*t).*s
16、in(2*pi*f3*t); x=x+10*u(1+i*N:i*N+N); x=x/max(x); x=x.*w; X=fft(x,N); X=abs(X); Y=Y+X;end Y=Y/1000; subplot(212);plot(f(1:N/2),Y(1:N/2);title(正弦加白噪声后与w时域相乘的FFT 1000次平均)xlabel(Hz)2、clear all;%三个正弦信号相加，分段函数，进行频谱分析 % 产生三个正弦相加的函数； N=512; f0=10;fs=100.0; t=0:N-1; x=1.0*sin(2*pi*f0*t/fs)+1.0*sin(2*pi*2*f0
17、*t/fs)+1.0*sin(2*pi*3*f0*t/fs); subplot(211); plot(t(1:N),x(1:N);title(x(t) %加窗 w=1-1.93*cos(2*pi*t/N)+1.29*cos(4*pi*t/N)-0.388*cos(6*pi*t/N)+0.0322*cos(8*pi*t/N); %w=1.0-cos(2*pi*t/N); x=w.*x;%加窗等于时域点乘 % 应用FFT 求频谱； f=0:fs/N:fs/N*(N-1); X=fft(x,N);%先点乘再进行傅里叶变换 X=abs(X)/N; subplot(212); plot(f(1:N/2)
18、,X(1:N/2);title(x(t)加窗之后的傅里叶变换) xlabel(Hz) %分段函数- M=170;L=N-2*M; x(1:M)=1.0*sin(2*pi*f0*t(1:M)/fs); x(M+1:2*M)=1.0*sin(2*pi*2*f0*t(1:M)/fs); x(2*M+1:N)=1.0*sin(2*pi*3*f0*t(1:L)/fs); figure; subplot(211); plot(t(1:N),x(1:N);title(x(t)为分段函数) %w=1-1.93*cos(2*pi*t/M)+1.29*cos(4*pi*t/M)-0.388*cos(6*pi*t/
19、M)+0.0322*cos(8*pi*t/M); %x(1:M)=w(1:M).*x(1:M); X=fft(x,N); X=abs(X)/N; subplot(212); plot(f(1:N/2),X(1:N/2);title(x(t)频谱分析) xlabel(Hz)3、%采样长度不同对FFT的影响-clear all;% 观察数据长度N的变化对DTFT分辨率的影响 f1=2;f2=2.02;f3=2.07;fs=10; w=2*pi/fs; N=256;%N=256 x=sin(w*f1*(0:N-1)+sin(w*f2*(0:N-1)+sin(w*f3*(0:N-1); f=0:fs/
20、N:fs/2-1/N; X=fft(x); X=abs(X); subplot(221); plot(f(45:60),X(45:60);title(N=256);grid on; subplot(223) plot(f(1:N/2),X(1:N/2);title(N=256);grid on; xlabel(Hz)% N=N*4;%N=1024 x=sin(w*f1*(0:N-1)+sin(w*f2*(0:N-1)+sin(w*f3*(0:N-1); f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X); subplot(222) plot(f(45*4:4*60),X
21、(4*45:4*60);title(N=1024);grid on; xlabel(Hz)4、%补零的影响-clear;% 计算长度为N的原始信号的DTFT f1=2.67;f2=3.75;f3=6.75;fs=20;w=2*pi/fs; N=16; x=sin(w*f1*(0:N-1)+sin(w*f2*(0:N-1)+pi/2)+sin(w*f3*(0:N-1); f=0:fs/N:fs/2-1/N; X=fft(x); X=abs(X); f=fs/N*(0:N/2-1); subplot(221) stem(f,X(1:N/2),.);title(不补零);grid on; xlabe
22、l(Hz) % 在数据末补N个零 x(N:2*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:N-1)/(2*N); subplot(222) stem(f,X(1:N),.);title(补N个零);grid on; xlabel(Hz) % 在数据末补7*N个零 x(N:8*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:4*N-1)/(8*N); subplot(223) stem(f,X(1:4*N),.);title(补7N个零);grid on; xlabel(Hz) % 在数据末补29*N个零 x(N:30*N-1)=0; X=ff
23、t(x); X=abs(X); f=fs*(0:15*N-1)/(30*N); subplot(224) plot(f,X(1:15*N);title(补29N个零);grid on; xlabel(Hz)5、%x(n)是两个正弦信号和一个白噪声相加，FFT和IFFT- clear all; % 产生两个正弦加白噪声； N=256; f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+randn(1,N); % 应用FFT 求频谱； subplot(3,1,1); plot(x(1:
24、N/4);title(x(n) f=-0.5:1/N:0.5-1/N; X=fft(x); y=ifft(X); subplot(3,1,2); plot(f,fftshift(abs(X);title(fft(x) subplot(3,1,3); plot(real(x(1:N/4); title(x(n)实部)6、%fftfilt和conv比较n=0:2;h=1./2.n;for i=1:51 x(i)=(i-1)/5;end for i=52:100 x(i)=20-(i-1)/5;end y=fftfilt(h,x); % y=x*hz=conv(h,x);subplot(311)ti
25、tle(x(t)hold;plot(x);subplot(312)title(fftfilt 叠接相加法)hold;plot(y);subplot(313)plot(z);axis(0,100,0,20);title (conv 卷积)7、clear;%补零后对频谱的影响% 计算长度为N的原始信号的DFT f0=50;%信号频率 fs=200;%采样频率 N=16;%抽样点数 w=2*pi/fs;%数字角频率 x=sin(w*f0*(0:N-1); X=fft(x); X=abs(X); f=fs/N*(0:N/2-1);%N/2的数据 subplot(211) stem(f,X(1:N/2)
26、,.);title(未补零);grid on; xlabel(Hz) % 在数据末补N个零 x(N:2*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:N-1)/(2*N); subplot(212) stem(f,X(1:N),.);title(补N个零);grid on; xlabel(Hz)8、%抽样频率不同的影响% 计算长度为N的原始信号的DFT% fs=100 f0=50;fs=100;w=2*pi/fs; N=16; x=sin(w*f0*(0:N-1); X=fft(x); X=abs(X); f=fs/N*(0:N/2-1); subplot(221)
27、stem(f,X(1:N/2),.);grid on; xlabel(Hz) % fs=150 f0=50;fs=150;w=2*pi/fs; N=16; x=sin(w*f0*(0:N-1); X=fft(x); X=abs(X); f=fs/N*(0:N/2-1); subplot(222) stem(f,X(1:N/2),.);grid on; xlabel(Hz) % fs=200 f0=50;fs=200;w=2*pi/fs; N=16; x=sin(w*f0*(0:N-1); X=fft(x); X=abs(X); f=fs/N*(0:N/2-1); subplot(223) st
28、em(f,X(1:N/2),.);grid on; xlabel(Hz)9、%周期延拓clear all; M=128;N=1024; alf=-3.0; f0=10;fs=100.0; %u=rand(1,5000);u=u-mean(u); t=0:1/fs:(M-1)/fs; x1=1.0*exp(alf*t).*sin(2*pi*f0*t); t=0:1/fs:(N-1)/fs; for i=0:7 x(M*i+1:M*i+M)=x1(1:M); end x=x-mean(x); subplot(211); plot(t,x); % 应用FFT 求频谱； f=0:fs/N:fs/N*(
29、N-1); X=fft(x,N); X=abs(X); %X=20*log10(X); subplot(212); plot(f(1:N/2),X(1:N/2); % stem(f,X,.);grid on; xlabel(Hz)10、%正弦加白噪声并FFT频谱分析 % 产生4096点的两个正弦加白噪声； N=4096; f1=2.1;f2=2.2;fs=5; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+10*randn(1,N); % 应用FFT 求频谱； f=fs/N*(0:N/2-1); X=fft(x,
30、N); X=abs(X); subplot(211); stem(f,X(1:N/2),.);title(x(n);grid on; xlabel(Hz) f=-0.5:1/N:0.5-1/N; subplot(212); stem(f,fftshift(X),.);title(FFT x(n);grid on;11、clear;%补零对频谱分析的影响 f1=10.8;f2=11.75;f3=12.55;fs=40;w=2*pi/fs; N=64; x=sin(w*f1*(0:N-1)+sin(w*f2*(0:N-1)+pi/2)+sin(w*f3*(0:N-1); X=fft(x); X=a
31、bs(X); f=fs/N*(0:N/2-1); subplot(221) stem(f,X(1:N/2),.);title(补0个零);grid on; xlabel(Hz) % 在数据末补3N个零 x(N:4*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:2*N-1)/(4*N); subplot(222) stem(f,X(1:2*N),.);title(补3N个零);grid on; xlabel(Hz) % 在数据末补7*N个零 x(N:8*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:4*N-1)/(8*N); subplot(2
32、23) stem(f,X(1:4*N),.);title(补7N个零);grid on; xlabel(Hz) % 在数据末补15*N个零 x(N:16*N-1)=0; X=fft(x); X=abs(X); f=fs*(0:8*N-1)/(16*N); subplot(224) plot(f,X(1:8*N);title(补15N个零);grid on; xlabel(Hz) 12、%窗函数的影响clear; %窗函数对频谱的影响 f1=10.8;f2=11.75;f3=12.55; fs=40; w=2*pi/fs; N=1024;t=0:N-1; f=fs/N*(0:N/2-1);%矩形
33、窗 x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); X=fft(x); X=2*abs(X)/N; subplot(221) stem(f,X(1:N/2),.);grid on; xlabel(Hz) title(矩形窗) %升余弦窗 x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=0.5-0.5*cos(2*pi*t/N); x=wt.*x; X=fft(x); X=2*2*abs(X)/N; subplot(222) stem(f,X(1:N/2),.);grid on; xlabel(Hz) title(升余弦窗)%平顶
34、窗 x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=1-1.93*cos(2*pi*t/N)+1.29*cos(4*pi*t/N)-0.388*cos(6*pi*t/N)+0.0322*cos(8*pi*t/N); x=wt.*x; X=fft(x); X=2*abs(X)/N; subplot(223) stem(f,X(1:N/2),.);grid on; xlabel(Hz) title(平顶窗) %改进升余弦窗 x=sin(w*f1*t)+sin(w*f2*t)+sin(w*f3*t); wt=0.54-0.46*cos(2*pi*t/N); x=w
35、t.*x; X=fft(x); X=1.852*2*abs(X)/N; subplot(224) stem(f,X(1:N/2),.);grid on; xlabel(Hz) title(改进升余弦窗)13、%频谱分析，平均clear all;%产生正弦加白噪声信号，求频谱，平均 % 产生两个正弦加白噪声； N=1024; f0=10;fs=100.0; u=rand(1,N);u=u-mean(u); t=0:1/fs:(N-1)/fs; x=1.0*sin(2*pi*f0*t)+0.1*sin(2*pi*2*f0*t)+10.0*u; % 应用FFT 求频谱； f=0:fs/N:fs/N*
36、(N-1); X=fft(x,N); X=abs(X); X=20*log10(X); subplot(211); plot(f(1:N/2),X(1:N/2);title(平均前) % stem(f,X,.);grid on; xlabel(Hz) Y=zeros(1,N); u=rand(1,1000*N);u=u-mean(u); for i=0:999 x=1.0*sin(2*pi*f0*t)+0.1*sin(2*pi*2*f0*t)+10*u(1+i*N:i*N+N); X=fft(x,N); X=abs(X); Y=Y+X;end Y=Y/1000; Y=20*log10(Y);s
37、ubplot(212);plot(f(1:N/2),Y(1:N/2);title(平均后)axis(0,50,0,60);xlabel(Hz)1高通 2带阻 3低通 4带通 5切比雪夫 带通1、1高通，双线性变换法，巴特沃斯（两种方法）% to design 高通high-pass DF with s=2/Ts(z-1)/(z+1)%-给出：通带下限频率，阻带上限频率，通带衰减，阻带衰减，采样频率-%-clear allwp=.8*pi;%通带下限频率ws=.44*pi;%阻带上限频率rp=3;%通带衰减rs=20;%阻带衰减Fs=2000;%采样频率%设计模拟滤波器% Firstly to
38、finish frequency prewarping;wap=2*Fs*tan(wp/2);%通带截止频率，模拟角频率，预畸变was=2*Fs*tan(ws/2);%阻带截止频率，模拟角频率，预畸变n,wn=buttord(wap,was,rp,rs,s);%求取模拟低通滤波器阶数，n是模拟低通滤波器的阶次,巴特沃斯滤波器阶数的选择，（最小阶数的选择）% Note: s!z,p,k=buttap(n);%设计模拟低通滤波器，极点，零点，增益b,a=zp2tf(z,p,k);%零极点增益模型转换为传递函数模型bt,at=lp2hp(b,a,wap);%模拟低通滤波器转换成高通滤波器，G(p)转
39、换成H（s）,wap为低通的截止频率%模拟转换成数字% Note: z=(2/ts)(z-1)/(z+1);ts=1,that is 2fs=1,fs=0.5;bz,az=bilinear(bt,at,Fs);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z），bz,az分别是H（Z）的分子分母多项式的系数h,w=freqz(bz,az,256,1);%求离散系统频响特性，H包含了离散系统频响在 0pi范围内N个频率等分点的值，w则包含了范围内N个频率等分点。figure(1);plot(w,20*log10(abs(h);grid;ylabel( High-pass DF);%高通
40、滤波器%直接设计高通滤波器，把上一种方法的buttord，buttap，lp2hp全包括进去% Directly to design H(z) by butter.mn,wn=buttord(wp/pi,ws/pi,rp,rs);%要知道阶数N和3dB截止频率矢量wn，buttord()计算这些参数b,a=butter(n,wp/pi,high);%设计高通滤波器h1,w1=freqz(b,a,256,1);%计算由n指定的频率点上，数字滤波器H(Z)的频率响应h1=20*log10(abs(h1);figure(2);plot(w,h1);grid;ylabel( High-pass DF);1、2高通，双线性变换法，巴特沃斯和切比雪夫（两种方法）（1）% 高通,双线性Z变换法，巴特沃斯,数字滤波器 s=2/Ts(z-1)/(z+1)clear all%给定参数fp=400;%通带下限频率fs=300;%阻带上限频率Fs=1000;%采样频率rp