测试信号处理（第3版）程序代码汇 第2-10章

例2.4参考代码：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%a.利用conv函数T=0.1;%%%时间步长t1=0:T:10;%%%时间序列f1=cos(t1);%%%信号f1t2=t1;f2=exp(t2)+exp(-2*t2);%%%信号f2f=T*conv(f1,f2);%%%计算卷积k0=t1(1)+t2(1);%%%卷积输出序列的起始k3=length(f1)+length(f2)-2;t=k0:T:(k0+T*k3);%%%卷积结果对应的时间向量%subplot(3,1,1);%%%绘制信号f1plot(t1,f1,'linewidth',2);title('f1(t)');subplot(3,1,2);%%%绘制信号f2plot(t2,f2,'linewidth',2);title('f2(t)');subplot(3,1,3);%%%绘制卷积结果plot(t,f,'linewidth',2);title('convolutionoff1(t)andf2(t)');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%b.编程实现T=0.1;%%%时间步长t1=0:T:10;%%%时间序列f1=cos(t1);%%%信号f1t2=t1;f2=exp(t2)+exp(-2*t2);%%%信号f2lf1=length(f1);%%%信号f1长度lf2=length(f2);%%%信号f2长度fork=1:lf1+lf2-1y(k)=0;%%%y赋初始值forii=max(1,k-(lf2-1)):min(k,lf1)y(k)=y(k)+f1(ii)*f2(k-ii+1);%%%信号相乘和求和endyzsappr(k)=T*y(k);%%%用乘和加运算来近似积分运算endt0=t1(1)+t2(1);%%%卷积输出序列起点t3=lf1+lf2-2;t=t0:T:(t0+t3*T);%%%卷积输出对应的时间序列subplot(3,1,1);%%%绘制信号f1的波形plot(t1,f1,'linewidth',2);title('f1(t)');subplot(3,1,2);%%%绘制信号f2的波形plot(t2,f2,'linewidth',2);title('f2(t)');subplot(3,1,3);%%%绘制卷积结果的波形plot(t,f,'linewidth',2);title('convolutionoff1(t)andf2(t)');例3.1%例3.1求解周期信号傅里叶级数的MATLAB程序%先求出信号傅里叶级数的系数symsTEtnx;pi=sym('pi');%创建符号对象a0=2/T*int(-E,t,-T/2,0)+2/T*int(E,t,0,T/2),A0=a0/2;%计算系数an=2/T*int(-E*cos(2*pi*n*t/T),t,-T/2,0)+2/T*int(E*cos(2*pi*n*t/T),t,0,T/2)bn=2/T*int(-E*sin(2*pi*n*t/T),t,-T/2,0)+2/T*int(E*sin(2*pi*n*t/T),t,0,T/2)%求出E=1时各谐波分量的幅度symsEAncbn=2*E*c;E=1;c=-(cos(n*pi)-1)/(n*pi);bn=subs(bn,{sym('E'),sym('c')},{E,c})%符号变量置换%画出频谱图A=[00000000000];%存储11个谐波分量幅度的数组n=1:1:11;forn=1:11bn=2*(-cos(n*pi)+1)/(n*pi);A(n)=double(vpa(bn));%对任意精度的符号类数据进行规范An=A(n)endx1=1;x2=11;x=x1:x2;stem(x,A,'r','filled')%画出11个谐波分量的幅度谱例3.3%产生周期为2π的方波信号t=0:1/10:10;y=square(2*pi*(0.05*pi)*t);subplot(2,2,1);plot(t,y);xlabel('(a)方波波形')axis([010-1.51.5])%基波信号t=0:1/10:10;y=sin(t);subplot(2,2,2);plot(t,y);xlabel('(b)基波波形')%基波+3次谐波合成的波形t=0:1/10:10;y=sin(t)+sin(3*t)/3;subplot(2,2,3);plot(t,y);xlabel('(c)基波+3次谐波')%基波+3次+5次+7次+9次谐波合成的波形t=0:1/10:10;y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;subplot(2,2,4);plot(t,y);xlabel('(d)基波+3次+5次+7次+9次谐波') %最高谐波阶次为19时的吉伯斯现象t=0:31/1000:5;y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:19,x=x+sin(k*t)/k;y((k+1)/2,:)=x;end;pause,plot(y(1:9,:)),pause,mesh(y),pauseclc%最高谐波阶次为45时的吉伯斯现象y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:45,x=x+sin(k*t)/k;y((k+1)/2,:)=x;end;pause,plot(y(1:9,:)),pause,mesh(y),pauseclc例3.10%例3.10中傅里叶变换的MATLAB程序symsEtaw;f=E*Heaviside(t+a)-E*Heaviside(t-a)F=fourier(f);F1=simplify(F)%得到F的简化形式E=1;a=1;F2=subs(F1,{sym('E'),sym('a')},{Ea})%将F1中的变量替换为新变量ezplot(abs(F2),[-3*pi,3*pi])title('F(Ω)')xlabel('Ω')f2=ifourier(F2)运行结果为F=(E*(sin(a*w)+cos(a*w)*1i))/w-(E*(cos(a*w)*1i-sin(a*w)))/w例3.14%例3.14高斯信号的频谱symstwsigma=sym('sigma','positive');f=exp(-t^2/(2*sigma^2))/(sqrt(2*pi)*sigma);F=fourier(f,t,w);F1=simplify(F);digits(4);F2=vpa(F1)sigma=1;F3=subs(F2,'sigma',sigma)ezplot(F3,[-2*pi,2*pi]),gridontitle('F(Ω)');xlabel('Ω')例3.15%信号最高频率wm＝2ws（ws:采样频率）临界采样的情况wm=1;%信号最高频率（带宽）wc=wm;%预处理用低通滤波器截止频率Ts=pi/wm;%采样间隔ws=2*pi/Ts;%变为采样角频率n=-100:100;%采样点数nTs=n*Ts;f=sinc(nTs/pi);%对信号Sa(t)的抽样信号Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));%信号的重构t1=-10:0.5:10;f1=sinc(t1/pi);%Sa(t)的连续波形表示%画出抽样信号与抽样信号的重构波形subplot(2,1,1);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的临界抽样');subplot(2,1,2);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由抽样信号重构sa(t)');grid%欠抽样情况，抽样频率ws与信号最高频率wm之间，wm>2ws.wm=1;wc=1.1*wm;Ts=2*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Ts;f=sinc(nTs/pi);Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));%求解重构信号与原信号的误差t1=-10:0.5:10;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的欠抽样信号');subplot(312);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由欠抽样信号重构的sa(t)');grid;subplot(313);plot(t,error);xlabel('t');ylabel('error(t)');title('欠抽样信号与原信号的误差error(t)')%过抽样情况，wm<2ws.wm=1;wc=1.1*wm;Ts=0.5*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Ts;f=sinc(nTs/pi);Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));t1=-10:0.5:10;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的采样信号');subplot(312);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由过抽样信号重构的sa(t)');grid;subplot(313);plot(t,error);xlabel('t');ylabel('error(t)');title('过抽样信号与原信号的误差error(t)')例4.17%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%symst;%%%定义符号变量f=cos(2*t);%%%信号表达式F=laplace(f);%%%调用MATLAB函数对x做拉普拉斯变换fl=cos(2*(t-1));%%%信号xl表达式Fl=laplace(fl);%%%对xl做拉普拉斯变换例4.18%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%symss;%%%定义符号变量F=5*(s+2)*(s+5)/(s*(s+1)*(s+3));%%%信号表达式f=ilaplace(F);%%%逆拉普拉斯变换例4.19%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%A=[83-21];%%%分子系数B=[10-7-6];%%%分母系数z=roots(A);%%%解得零点p=roots(B);%%%解得极点x=max(abs([z;p]));%%%取零极点的最大值x=x+0.1;y=x;figure;%%%绘制零极点图holdon;%%%使多次绘图同时保留plot([-xx],[00],'--');%%%绘制实轴（或x轴）plot([00],[-yy],'--');%%%绘制虚轴（或y轴）plot(real(z),imag(z),'bo',real(p),imag(p),'kx');xlabel('RealPart');ylabel('ImaginaryPart');axis([-xx-yy]);例4.20%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%symst;%%%定义符号变量f=exp(-2*t)*heaviside(t);%%%信号表达式F=laplace(f);%%%对信号做拉普拉斯变换%绘制零极点图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%A=[1];%%%多项式分子B=[12];%%%多项式分母z=roots(A);%%%解得零点p=roots(B);%%%解得极点x=max(abs([z;p]));%%%取零极点的最大值x=x+0.1;y=x;figure;%%%绘制零极点图holdon;plot([-xx],[00],'--');%%%绘制实轴plot([00],[-yy],'--');%%%绘制虚轴plot(real(z),imag(z),'bo',real(p),imag(p),'kx');xlabel('RealPart');ylabel('ImaginaryPart');axis([-xx-yy]);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%length=100;%%%定义仿真的长度xs=zeros(length,length);%%%仿真区域初始化delta_t=0.1;%%%时间步长formm=1:length%%%对仿真区域每一个点fornn=1:lengthdelta_xs=0;%%%保存积分运算的中间结果s=(mm/10-5)+1j*(nn/10-5)%%%s的取值范围为[-5-5i,5+5i]xs(mm.nn)=1./(s+2);%%%信号的拉普拉斯变换ifmm<30%%%相当于sigma<-2fortt=0:delta_t:50delta_xs=delta_xs+exp(-(2+s)*tt);%%%对delta_xs累加endxs(mm,nn)=delta_xs*delta_t;%%%对xt进行积分endendendxs=log10(abs(xs));figure;mesh(1:length,1:length,xs);set(gca,'xtick',[0102030405060708090100]);set(gca,'xticklabel',{'-5','-4','-3','-2','-1','0','1','2','3','4','5'});set(gca,'ytick',[0102030405060708090100]);set(gca,'yticklabel',{'-5','-4','-3','-2','-1','0','1','2','3','4','5'});xlabel('j\omega');ylabel('sigma');zlabel('log10(X(s))');例5.3%MATLABPROGRAMExample5.3b=input('差分方程左端系数b=[b(1),…,b(nb)]=');a=input('差分方程右端系数a=[a(1),...,a(na)]=');x=input('输入信号序列x=');nb=length(b);na=length(a);nx=length(x);s=['现时刻前',int2str(nb-1),'点ym的值=[ym(1),...,ym(na-1)='];ym=zeros(1,nb+nx-1);%先预设ym序列长度为nb+nx,并初始化为0ym(1:nb-1)=input(s);%对ym序列的前nb个点赋初值xm=[zeros(1,nb-1),x];%给序列xm赋初值forn=nb:nb+nx-1%n以ym的起点为初值ys=ym(n-1:-1:n-nb+1);%得出ysxs=xm(n-1:-1:n-na);%得出xsym(n)=(a*xs'-b(2:nb)*ys')/a(1);%递推求ymend%把ym左移nb-1位,求得yy=ym(nb:nb+nx-1);stem(y),gridline([0,60],[0,0])%给出起点、终点坐标，画横坐标例6.2%x(n)的图形表示n=-5:5;x=(-0.95).^n;subplot(3,1,1);stem(n,x);title('序列x(n)的图形');xlabel('n');ylabel('x(n)')%x(n)的傅里叶变换k=-200:200;w=(pi/100)*k;%横坐标(频率轴)的示值缩小比例系数X=x*(exp(-j*pi/100)).^(n'*k);magX=abs(X);angX=angle(X);subplot(3,1,2);plot(w/pi,magX);gridon;title('x(n)的幅度谱');axis([-2,2,0,15]);xlabel('频率ω(单位:π)');ylabel('|X|')subplot(3,1,3);plot(w/pi,angX/pi);gridon;title('x(n)的相位谱');axis([-2,2,-1,1]);xlabel('频率ω(单位:π)');ylabel('φ/π')例6.4%构成并画出受噪声污染的信号clear,randn('state',0);%正态随机数发生器置0t=linspace(0,10,512);%生成均匀采样数组y=6*cos(6*t)-4*sin(12*t)+6*randn(size(t));%生成波形和%标准差为6，标准偏差为0的正态随机噪声subplot(311);plot(t,y)xlabel('t');%标注坐标轴ylabel('信号值f(t)');%计算并画出受噪声污染信号的幅度谱Y=fft(y);%计算信号的DFTTs=t(2)-t(1)%计算信号的采样周期Ws=2*pi/Ts;%计算信号的采样频率Wn=Ws/2w=linspace(0,Wn,length(t)/2);%频率轴的频率间隔Ya=abs(Y(1:length(t)/2));%频谱幅度subplot(312);plot(w,Ya)xlabel('\omega')%标注坐标轴ylabel('频谱幅值F(\omega)')ii=find(w<=20);%在有效信号频率段寻找数组中非零元素序号subplot(313);plot(w(ii),Ya(ii))xlabel('\omega')%标注坐标轴ylabel('局部放大谱幅F(\omega)')grid例6.5clear%序列的表示a=ones(1,15);a([1,2,3,4,5])=0;%序列a(n)b=ones(1,12);b([1,2,3])=0;%序列b(n)%序列直接卷积c=conv(a,b);%直接卷积运算%应用FFT及DFT性质，进行序列的快速卷积M=32;%加长序列符合快速卷积要求AF=fft(a,M);%序列a(n)FFT得AFBF=fft(b,M);%序列b(n)FFT得BFCF=AF.*BF;%AF和BF相乘得CFcc=real(ifft(CF));%CF进行IFFT%求直接卷积和快速卷积两种计算结果之间的误差nn=0:(M-1);c(M)=0;error=c-cc;%图形表示出序列、直接卷积和快速卷积以及两种卷积结果的误差subplot(411);stem([1:15],a),ylabel('a(n)');subplot(412);stem([1:12],b),ylabel('b(n)');subplot(413),stem(nn,c,'fill'),grid,axis([0,31,0,9])ylabel('a(n)*b(n)')subplot(414),stem(nn,error,'fill'),axis([0,31,-1,1])xlabel('n'),ylabel('误差')例6.6N=64;T=5/N;t=0:T:5;x=exp(-1.5*t);%连续时间信号x(t)Xx=T*fft(x,N);%x(t)的近似频谱n=0:1:63;x0=1/T.*ifft(Xx,N);%重构x(n)%分别画出x(t)波形、幅谱X(k)和x(n)subplot(311);plot(t,x);xlabel('t');ylabel('x(t)')subplot(312);stem(n,Xx);xlabel('k'),ylabel('X(k)')subplot(313)stem(n,x0);xlabel('n'),ylabel('x(n)')例6.7%example6.7MATLABPROGRAMN=256;%采样点数dt=0.0004%采样时间间隔wn=500;%二阶系统的固有频率seta=0.47;%系统阻尼比dw=2*pi/(N*dt);%频率间隔a=wn^2;b=[1,2*seta*wn,a];t=[0:dt:(N-1)*dt];c=step(a,b,t);%求二阶系统的阶跃响应w=[0:dw:(N-1)*dw];[mag,phase]=bode(a,b,w);%计算二阶系统理论频率特性ycw=fft(c);%求系统阶跃响应瞬态分量的傅里叶变换Re=real(ycw);%ycw的实部Im=imag(ycw);%ycw的虚部fori=1:NRw(i,1)=1-Im(i,1)*(i-1)*dw*dt;Iw(i,1)=Re(i,1)*(i-1)*dw*dt;end%计算频率特性的实部和虚部分量ffw=Rw+Iw*sqrt(-1);Aw=abs(ffw)%系统幅频特性semilogx(w,20*log10(mag),'r-')%理论幅频特性曲线axis([100,10000,-30,10])text(600,12,'对数幅频特性')ylabel('A(ω)(dB)');holdonsemilogx(w,20*log10(Aw))%由阶跃响应求得的幅频特性曲线axis([100,10000,-30,10])xlabel('ω(1/s)');gridon例7.11%MATLABPROGRAMofExample7.11b=[0,1];a=[3,-4,1];[R,p,k]=residuez(b,a)例7.12%例7.12中求解z变换的MATLAB程序symsaznTf=a^(n*T)F=factor(ztrans(f))%做因式分解处理%若T=1,a=1求其z变换symsknwzF=ztrans(1^n)例7.13symsznabF=k*z^2/((z-a)*(z-b))f=iztrans(F)%若k=1,a=1,b=0.5时,求其z反变换symsnzF=z^2/((z-1)*(z-0.5))f=iztrans(F)例7.18%Exampe7.18MATLABProgram%画出单位响应h(n)clear;n=0:7;hn=(0.9).^n;subplot(3,1,1);stem(n,hn,'k')xlabel('n');ylabel('h(n)');%画出频率响应（幅频和相频）k=0:200;w=(pi/100)*k;Hew=hn*(exp(-j*pi/100).^(n'*k));mag_H=abs(Hew);ang_H=angle(Hew);subplot(3,1,2);plot(w/pi,mag_H,'k');xlabel('w');ylabel('H(w)');subplot(3,1,3);plot(w/pi,ang_H/pi,'k');xlabel('w');ylabel('\Phi');例8.4%例8.4的MATLAB程序n=0:0.01:2;fori=1:4switchicase1N=2;case2N=5;case3N=10;case4N=20;end[z,p,k]=buttap(N);%巴特沃思滤波器原型设计函数，n:阶数；%z,p,k：滤波器零点、极点和增益[b,a]=zp2tf(z,p,k);%零极点增益转换为传递函数[H,w]=freqs(b,a,n);%模拟滤波器频率响应函数magH2=(abs(H)).^2;%幅度平方函数holdonplot(w,magH2);axis([0201]);endxlabel('w/wc');ylabel('|H(jw)|^2');grid例8.5%绘制切比雪夫低通滤波器幅度平方函数的NATLAB程序n=0:0.01:2;fori=1:4switchicase1N=2;case2N=4;case3N=6;case4N=8;endRp=1;%滤波器通带波纹系数[z,p,k]=cheb1ap(N,Rp);%设计切比雪夫模拟低通滤波器原型函数,%z,p,k分别为滤波器的零点、极点和增益[b,a]=zp2tf(z,p,k);%零点、极点和增益转换为传递函数[H,w]=freqs(b,a,n);%模拟滤波器的频率响应magH2=(abs(H)).^2;%幅度平方函数%Outputposplot=['22',num2str(i)];%定义字符串变量subplot(posplot)plot(w,magH2,'k');ylabel('|H(jw)^2');title(['N=',num2str(N)]);gridend例8.7%设计切比雪夫模拟低通滤波器的MATLAB程序wp=200*pi;ws=300*pi;Rp=1;%通带波纹Rs=16;%阻带衰减%计算滤波器阶数ebs=sqrt(10^(Rp/10)-1);%波纹系数A=10^(Rs/20);%A为参变量Wc=wp;%截止频率Wr=ws/wp;g=sqrt(A*A-1)/ebs;%g为参变量N1=log10(g+sqrt(g*g-1))/log10(Wr+sqrt(Wr*Wr-1));%滤波器阶数计算N=ceil(N1);%N应取整例8.8%设计巴特沃思带通滤波器MATLAB程序%滤波器设计指标wp=[20003000]*2*pi;ws=[15003500]*2*pi;Rp=1;Rs=100;%计算阶数和截止频率[N,Wn]=buttord(wp,ws,Rp,Rs,'s');Fc=Wn/(2*pi);%计算滤波器传递函数多项式系数[b,a]=butter(N,Wn,'s');%画出滤波器幅频特性w=linspace(1,4000,1000)*2*pi;%生成线性等间隔的向量H=freqs(b,a,w);magH=abs(H);phaH=unwrap(angle(H));plot(w/(2*pi),20*log10(magH),'k');xlabel('频率(Hz)');ylabel('幅度(dB)');title('巴特沃思模拟滤波器')gridon例8.9%设计切比雪夫带阻滤波器MATLAB程序%滤波器设计指标wp=[20003000]*2*pi;ws=[15003500]*2*pi;Rp=1;Rs=60;%计算阶数和截止频率[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');%计算滤波器传递函数多项式系数[b,a]=cheby1(N,Rp,Wn,'stop','s');%画出滤波器幅频特性w=linspace(1,4000,1000)*2*pi;%生成线性等间隔的向量H=freqs(b,a,w);%相位展开magH=abs(H);phaH=unwrap(angle(H));plot(w/(2*pi),20*log10(magH),'k');xlabel('频率(Hz)');ylabel('幅度(dB)');title('切比雪夫模拟滤波器')gridon例8.15%例8.15中的MATLAB程序%定义滤波器的性能指标wp=2000*2*pi;ws=4000*2*pi;Rp=3;Rs=15;Fs=20000;Nn=128;%给出滤波器的序列点数%计算滤波器阶数和截止频率[N,Wn]=buttord(wp,ws,Rp,Rs,'s');%设计模拟滤波器原型[z,p,k]=buttap(N);[Bap,Aap]=zp2tf(z,p,k);%将系统函数的零极点形式变为分子分母多项式形式[b,a]=lp2lp(Bap,Aap,Wn);%低通至低通模拟滤波器的变换%冲激响应不变法实现模拟数字滤波器的变换[bz,az]=impinvar(b,a,Fs);freqz(bz,az,Nn,Fs)例8.19%例8.19双线性变换法设计数字低通滤波器MATLAB程序%给定滤波器指标wp=0.2*pi;ws=0.3*pi;Rp=3;Rs=15;Ts=0.001;Nn=128;%数字频率模拟频率非线性变换Wp=(2/Ts)*tan(wp/2);Ws=(2/Ts)*tan(ws/2);%计算滤波器阶次和截止频率[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s');%设计模拟原型[z,p,k]=buttap(N);[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,Wn);%双线性变换法设计数字滤波器[bz,az]=bilinear(b,a,1/Ts);freqz(bz,az,Nn,1/Ts)例8.22%给出滤波器设计要求Fs=10000;wp=[20003000]*2/Fs;ws=[15003500]*2/Fs;Rp=1;Rs=20;Nn=128;%计算模拟原型的阶次和截止频率[N,Wn]=cheb1ord(wp,ws,Rp,Rs);[b,a]=cheby1(N,Rp,Wn);%得到频响特性freqz(b,a,Nn,Fs)例8.23%例8.23中设计巴特沃思高通数字滤波器MATLAB程序%给定滤波器设计指标Fs=5000;wp=2000*2/Fs;ws=1500*2/Fs;Rp=1;Rs=20;Nn=128;%计算滤波器阶次和截止频率[N,Wn]=buttord(wp,ws,Rp,Rs);%设计数字滤波器[b,a]=butter(N,Wn,'high');%数字滤波器频响特性freqz(b,a,Nn,Fs)例8.25%滤波器设计要求wn=[0.280.58];N=50;%使用FIR1类函数计算并求出滤波器特性b=fir1(2*N,wn);freqz(b,1,512)例8.27%例8.27中设计FIR滤波器MATLAB程序N=20;alpha=(N-1)/2;l=0:N-1;wl=(2*pi/N)*l;Hrs=[1,1,1,zeros(1,15),1,1];Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1;angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];H=Hrs.*exp(j*angH);h=real(ifft(H,N));[db,mag,pha,grd,w]=freqz_m(h,1);[Hr,ww,a,L]=Hr_Type2(h);subplot(1,1,1)subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);axis([0,1,-0.1,1.1]);title('频率样本:N=20')xlabel('频率(单位:pi)');ylabel('Hr(k)')set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,2);stem(l,h);axis([-1,N,-0.1,0.3])title('单位抽样响应');ylabel('h(n)');text(N+1,-0.1,'n')subplot(2,2,3);plot(ww/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.2,1.2]);title('振幅响应')xlabel('频率(单位:pi)');ylabel('Hr(w)')set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,4);plot(w/pi,db);axis([0,1,-200,10]);gridtitle('幅度响应');xlabel('频率(单位:pi)');ylabel('分贝数');set(gca,'XTickMode','Manual','XTick',[0;0.2;0.3;1]);set(gca,'YTickMode','Manual','YTick',[-16;0]);set(gca,'YTickLabelMode','manual','YTickLabels',['16';'0'])%function[db,mag,pha,grd,w]=freqz_m(b,a);%------------------------------------%[db,mag,pha,grd,w]=freqz_m(b,a);%db：0-π区间内的相对振幅(dB)%mag：0-π区间内的绝对振幅%pha：0-π区间内的相移%grd：0-π区间内的群延迟%w：0-π区间内的501个抽样点频率%b：系统函数H(z)中分子多项式系数(对FIRb=h)%a：系统函数H(z)中分母多项式系数(对FIR:a=[1])%[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);grd=grpdelay(b,a,w);扩展函数Hr_Type2的MATLAB程序为function[Hr,w,b,L]=Hr_Type2(h);%-----------------------------------------------------------%[Hr,w,b,L]=Hr_Type2(h)%Hr:频响振幅%w：0-π区间内计算Hr的频率点%b：低通滤波器的系数%L：Hr的阶次%h：低通滤波器的单位抽样响应%M=length(h);L=M/2;b=2*[h(L:-1:1)];n=[1:1:L];n=n-0.5;w=[0:1:500]'*pi/500;Hr=cos(w*n)*b';例8.28%例8.28求解三种结构形式及转换的滤波系数MATLAB程序%b=[1,-3,11,-27,18];a=[16,12,2,-4,-1];[b0,B,A]=dir2cas(b,a);[C,B1,A1]=dir2par(b,a);N=24;n=1:N+1;formatlong;delta=impseq(0,0,N);h1=filter(b,a,delta);h2=casfiltr(b0,B,A,delta);h3=parfiltr(C,B1,A1,delta);figure(1)subplot(3,1,1)plot(n,h1)title('直接型')subplot(3,1,2)plot(n,h2)title('级联型')subplot(3,1,3)plot(n,h3)title('并联型')function[b,a]=dir2cas(b0,B,A);%将直接型结构变成级联型结构,求系统函数的系数%b0为常数;B、A分别为级联型系统函数的分子、分母%b、a分别为直接型系统函数的分子、分母%计算增益系数b0b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;%补零,使b和a等长M=length(b);N=length(a);ifN>Mb=[b,zeros(1,N-M)];elseifM>Na=[a,zeros(1,M-N)];N=M;elseNM=0;end%计算多项式的根K=floor(N/2);B=zeros(K,3);A=zeros(K,3);ifK*2==Nb=[b0];a=[a,0];end%用函数cplxpair把根以共轭复根对的次序排列,并用poly函数转换成二阶多项式broots=cplxpair(roots(b));aroots=cplxpair(roots(a));fori=1:2:2*KBrow=broots(i:1:i+1,:);Brow=real(poly(Brow));B(fix((i+1)/2),:)=Brow;Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));A(fix((i+1)/2),:)=Arow;end%直接型到并联型的转换%--------------------------------------%function[C,B,A]=dir2par(b,a)%C=当length(b)>=length(a)时的多项式部分%B=包含各系数bk的实系数矩阵%A=包含各ar系数的实系数矩阵%b=直接型的分子多项式系数%a=直接型的分母多项式系数function[C,B,A]=dir2par(b,a)M=length(b);N=length(a);[r1,p1,C]=residuez(b,a);p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);ifK*2==N;%Neven,orderofA(z)odd,onefactorisfirstorderfori=1:2:N-2Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix((i+1)/2),:)=real(Brow);A(fix((i+1)/2),:)=real(Arow);end[Brow,Arow]=residuez(r(N-1),p(N-1),[]);B(K,:)=[real(Brow)0];A(K,:)=[real(Arow)0];elsefori=1:2:N-1Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix((i+1)/2),:)=real(Brow);A(fix((i+1)/2),:)=real(Arow);endendfunctiony=casfiltr(b0,B,A,x);%IIR%-----------------------------------------------%y=casfiltr(b0,B,A,x);%y=输出序列%b0=级联型的增益系数%B=包含各系数bk的实系数矩阵%A=包含各ar系数的实系数矩阵%x=输入序列%[K,L]=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=x;fori=1:1:Kw(i+1,:)=filter(B(i,:),A(i,:),w(i,:));endy=b0*w(K+1,:);functiony=parfiltr(C,B,A,x);%IIR滤波器的并联型实现%----------------------------------------%[y]=parfiltr(C,B,A,x);%y=输出序列%C=当M>=N时(FIR)的多项式部分%B=包含各bk系数的实系数矩阵%A=包含各ar系数的实系数矩阵%x=输入序列%[K,L]=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=filter(C,1,x);fori=1:1:Kw(i+1,:)=filter(B(i,:),A(i,:),x);endy=sum(w);function[x,n]=impseq(n0,n1,n2)%产生单位抽样序列n=[n1:n2];x=[(n-n0)==0];functionI=cplxcomp(p1,p2)%比较两个包含同样标量元素,但可能有不同下标的复数对,返回其数组下标%这一函数必须用在cplxpair函数之后,以便对频率极点矢量及其响应的留数矢量重新排序p2=cplxpair(p1)I=[];forj=1:1:length(p2)fori=1:1:length(p1)if(abs((p1(i))-p2(j))<0.0001)I=[I,i];endendendI=I';例8.29%例8.29中计算FIR滤波器级联型结构的MATLAB程序b=[1,0,0,0,16+1/16,0,0,0,1];[b0,B,A]=dir2cas(b,1)例9.6%例9.6中周期图法计算信号功率谱的Matlab程序clfFs=2000;%情况1:数据长度N1=256N1=256;N1fft=256;n1=0:N1-1;t1=n1/Fs;f1=50;f2=100;xn1=sin(2*pi*f1*t1)+2*sin(2*pi*f2*t1)+randn(1,N1);Pxx1=10*log10(abs(fft(xn1,N1fft).^2)/N1);f1=(0:length(Pxx1)-1)*Fs/length(Pxx1);subplot(2,1,1)plot(f1,Pxx1)ylabel(′功率谱(dB)′);title(′数据长度N1=256′)grid%情况2:数据长度N2=1024N2=1024;N2fft=1024;n2=0:N2-1;t2=n2/Fs;f1=50;f2=100;xn2=sin(2*pi*f1*t2)+2*sin(2*pi*f2*t2)+randn(1,N2);Pxx2=10*log10(abs(fft(xn2,N2fft).^2)/N2);f2=(0:length(Pxx2)-1)*Fs/length(Pxx2);subplot(2,1,2)plot(f2,Pxx2)xlabel(‘频率(Hz)’);ylabel(‘功率谱(dB)’);title(‘数据长度N2=1024’)例9.8%基于功率谱估计的系统辨识%产生输入信号N=1024;nFFT=256;window=hanning(256);noverlap=128;dflag=′none′;Fs=1000;n=0:N-1;t=n/Fs;randn(′state′,0);xn=sin(2*pi*50*t)+randn(1,N);%构建系统：一个滤波器特性h=ones(1,10)/10;%计算系统输出yn=filter(h,1,xn);%进行系统函数估计[Txy,f]=tfe(xn,yn,nFFT,Fs,window,noverlap,dflag);%计算系统频率响应H=freqz(h,1,f,Fs);%比较响应结果subplot(2,1,1);plot(f,abs(H));xlabel(′频率（Hz）′);ylabel(′幅度′);title(′频率响应′);axis([050001]);gridsubplot(2,1,2);plot(f,abs(Txy));xlabel(′频率（Hz）′);ylabel(′幅度′);title(′频率响应′);axis([050001])grid例9.9%example9.9MatlabProgramN=1000;n=0:N-1;Fs=500;t=n/Fs;Lag=100;%相关信号的最大延迟量x1=sin(2*pi*10*t)+0.6*randn(1,length(t));%含白噪声的正弦信号x1[c,lags]=xcorr(x1,Lag,′unbiased′);%无偏自相关函数的计算subplot(2,2,1)%画出x1曲线plot(t,x1);xlabel(′t′);ylabel(′x1(t)′);title(′含白噪声的正弦信号x1′);gridsubplot(2,2,2);%画x1自相关曲线plot(lags/Fs,c);xlabel(′t′);ylabel(′Rxx1(t)′);title(′x1的自相关函数′);gridx2=randn(1,length(t));%发生白噪声x2\[c,lags\]=xcorr(x2,Lag,′unbiased′);%白噪声x2的无偏自相关函数subplot(2,2,3)%画x2曲线plot(t,x2);xlabel(′t′);ylabel(′x2(t)′);title(′白噪声x2′);gridsubplot(2,2,4)%画x2自相关曲线plot(lags/Fs,c)xlabel(′t′);ylabel(′Rxx2(t)′);title(′x2的自相关曲线′);grid例9.10%采用白噪声输入，进行系统辨识h=fir1(30,0.2,boxcar(31));x=randn(16384,1);y=filter(h,1,x);tfe(x,y,1024,[],[],512)例10.3t=0:0.01:2.10;%抽样频率为100Hzs1=sin(2*pi*30*t);%产生频率为30Hz的正弦信号s2=s1+0.5*[zeros(1,50)s1(1:161)];%加上回声，其幅度是原信号的一半%时域上延迟50个抽样周期即0.5秒c=cceps(s2);%用函数cceps求出倒谱subplot(121),plot(t,s2)%画出倒谱图xlabel('\fontsize{14}t')ylabel('\fontsize{14}\fontname{courier}信号+噪声的幅值')subplot(122),plot(t,c)xlabel('\fontsize{14}t')ylabel('\fontsize{14}\fontname{courier}倒谱')例10.4%%例10.4，利用FFT对信号进行STFT分析A=2;w0=10*pi*0.00001;%产生一个被分析的信号M=20000;n=0:M-1;x=A*sin(n.^2*w0);%%进行短时傅里叶变换L=200;P=L/2;N=256;Fs=10000;%%汉宁窗的窗口长度L，做FFT运算w=0.5*(1-cos(2*pi*(0:L-1)/(L-1)));Q=fix((M-P)/(L-P));forq=0:Q-1x0=x(q*(L-P)+1:q*(L-P)+L).*w;X(q+1,:)=fft(x0,N);endtn=((0:Q-1)*(L-P))/Fs;fk=(0:N/2)*Fs/N;zhX=X