数字信号处理(胡广书例题作业程序)

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(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);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);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=50000;u=rand(1,N);u=u-mean(u);a=sqrt(12*p);u1=u*a;power_u1=dot(u1,u1)/Nsubplot(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);plot(u(1:200));grid on;ylabel('u(n)');xlabel('n')subplot(212)hist(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*pi~2*pi');grid on;5、% 产生一chirp 信号；% chirp(T0,F0,T1,F1):% T0: 信号的开始时间；F0:信号在T0时的瞬时频率，单位为Hz；% T1: 信号的结束时间；F1:信号在T1时的瞬时频率，单位为Hz；%-----------------------------------------------------------------clear;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;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('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(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)');grid 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)') ;grid 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');ylabel('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两个信号%产生xu=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=x+u;x=x/max(x);%产生walf4=-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;%时域相乘，频域卷积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:999x=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=x+10*u(1+i*N:i*N+N); x=x/max(x);X=fft(x,N);X=abs(X);Y=Y+X;endY=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:999x=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=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;endY=Y/1000;subplot(212);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声后与w时域相乘的FFT 1000次平均') xlabel('Hz')24681012-1-0.500.51x(t) w(t)05101520253035404550102030x(t)频谱Hz024681012-0.50.51正弦加白噪声后与w 时域相乘05101520253035404550510正弦加白噪声后与w 时域相乘的FFTHz正弦加白噪声的FFT 1000次平均05101520253035404550Hz正弦加白噪声后与w时域相乘的FFT 1000次平均05101520253035404550Hz2、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*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),X(1:N/2));title('x(t)加窗之后的傅里叶变换')xlabel('Hz')x(t)0100200300400500600x(t)加窗之后的傅里叶变换05101520253035404550Hz%分段函数-----------------------------------------------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/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')x(t)为分段函数0100200300400500600x(t)频谱分析05101520253035404550Hz3、%采样长度不同对FFT的影响----------------------------------------------------clear all;% 观察数据长度N的变化对DTFT分辨率的影响f1=2;f2=2.02;f3=2.07;fs=10;w=2*pi/fs;N=256;%N=256x=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(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=1024x=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(4*45:4*60));title('N=1024');grid on;xlabel('Hz')1.52 2.5N=256N=256Hz1.52 2.5020*******N=1024Hz4、%补零的影响------------------------------------------------------------------- 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; xlabel('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=fft(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')510不补零Hz 0510补N 个零Hz 0510补7N 个零Hz510510补29N 个零Hz5、%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: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)实部')x(n)010203040506070fft(x)-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5x(n)实部0102030405060706、%fftfilt和conv比较n=0:2;h=1./2.^n;for i=1:51x(i)=(i-1)/5;endfor i=52:100x(i)=20-(i-1)/5;endy=fftfilt(h,x); % y=x*hz=conv(h,x);title('x(t)')hold;plot(x);subplot(312)title('f ftfilt 叠接相加法')hold;plot(y);subplot(313)plot(z);axis([0,100,0,20]);title ('conv 卷积')x(t)0102030405060708090100fftfilt 叠接相加法0102030405060708090100conv 卷积01020304050607080901007、clear;%补零后对频谱的影响% 计算长度为N的原始信号的DFTf0=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的数据stem(f,X(1:N/2),'.');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')未补零Hz补N个零Hz8、%抽样频率不同的影响% 计算长度为N的原始信号的DFT% fs=100f0=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)stem(f,X(1:N/2),'.');grid on;xlabel('Hz')% fs=150f0=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=200f0=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)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')0204060-14HzHzHz9、%周期延拓 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:7x(M*i+1:M*i+M)=x1(1:M);endx=x-mean(x);subplot(211);plot(t,x);% 应用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));% stem(f,X,'.');grid on;xlabel('Hz')024681012Hz10、%正弦加白噪声并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,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;x(n)HzFFT x(n)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=abs(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(223)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')补0个零Hz 05101520补3N 个零Hz 05101520补7N 个零Hz5101520010203040补15N 个零Hz12、%窗函数的影响 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); %矩形窗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('升余弦窗') %平顶窗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=wt.*x; X=fft(x);X=1.852*2*abs(X)/N; subplot(224)stem(f,X(1:N/2),'.');grid on; xlabel('Hz')title('改进升余弦窗')Hz 矩形窗Hz升余弦窗Hz平顶窗Hz改进升余弦窗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*(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:999x=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;endY=Y/1000;Y=20*log10(Y);subplot(212);plot(f(1:N/2),Y(1:N/2));title('平均后')axis([0,50,0,60]);xlabel('Hz')5101520253035404550204060平均前Hz 05101520253035404550204060平均后Hz1高通 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 finish frequency prewarping;wap=2*Fs*tan(wp/2);%通带截止频率，模拟角频率，预畸变 was=2*Fs*tan(ws/2);%阻带截止频率，模拟角频率，预畸变 [n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数，n 是模拟低通滤波器的阶次,巴特沃斯滤波器阶数的选择，（最小阶数的选择）[z,p,k]=buttap(n);%设计模拟低通滤波器，极点，零点，增益 [b,a]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型 [bt,at]=lp2hp(b,a,wap);%模拟低通滤波器转换成高通滤波器，G(p)转换成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 包含了离散系统频响在 0~pi 范围内N 个频率等分点的值，w 则包含了范围内N 个频率等分点。