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

实用文档/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-');gridon;ylabel('x(n)andy(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-');gridon;ylabel('x(n)andy(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)');gridon;%率减正弦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)');gridon;%谐波信号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)');gridon;%哈明窗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)');gridon;%采样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)');gridon;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));gridon;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));gridon;ylabel('u(n)');xlabel('n')subplot(212)hist(u,50);gridon;ylabel('histogramofu(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');gridon;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;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)');gridon;subplot(232);%绘制hstem(nh,h,'.k');xlabel('n');ylabel('h(n)');gridon;subplot(233);%绘制卷积stem(ny,y,'.k');xlabel('n');ylabel('y(n)');gridon;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)');gridon;subplot(223);plot((-Mlag:Mlag),rx1);title('x1自相关');gridon;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)');gridon;subplot(224);plot((-Mlag:Mlag),rx2);title('x2自相关');gridon;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)');gridon;subplot(212);plot(nrx,rx);xlabel('n');ylabel('rx(n)');gridon;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)');gridon;subplot(212);plot(nrx,rx);xlabel('n');ylabel('rx(n)');gridon;1、%产生信号，求卷积，FFT，求平均clearall;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,'.');gridon;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);%零均值白噪声fori=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('正弦加白噪声的FFT1000次平均')xlabel('Hz')Y=zeros(1,N);fori=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时域相乘的FFT1000次平均')xlabel('Hz')2、clearall;%三个正弦信号相加，分段函数，进行频谱分析%产生三个正弦相加的函数；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')%分段函数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')3、%采样长度不同对FFT的影响clearall;%观察数据长度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');gridon;subplot(223)plot(f(1:N/2),X(1:N/2));title('N=256');gridon;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');gridon;xlabel('Hz')4、%补零的影响clear;%计算长度为N的原始信号的DTFTf1=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('不补零');gridon;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个零');gridon;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个零');gridon;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个零');gridon;xlabel('Hz')5、%x(n)是两个正弦信号和一个白噪声相加，FFT和IFFTclearall;%产生两个正弦加白噪声；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)实部')6、%fftfilt和conv比较n=0:2;h=1./2.^n;fori=1:51x(i)=(i-1)/5;endfori=52:100x(i)=20-(i-1)/5;endy=fftfilt(h,x);%y=x*hz=conv(h,x);subplot(311)title('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的原始信号的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的数据subplot(211)stem(f,X(1:N/2),'.');title('未补零');gridon;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个零');gridon;xlabel('Hz')8、%抽样频率不同的影响%计算长度为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),'.');gridon;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),'.');gridon;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),'.');gridon;xlabel('Hz')9、%周期延拓clearall;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];fori=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,'.');gridon;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,N);X=abs(X);subplot(211);stem(f,X(1:N/2),'.');title('x(n)');gridon;xlabel('Hz')f=-0.5:1/N:0.5-1/N;subplot(212);stem(f,fftshift(X),'.');title('FFTx(n)');gridon;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个零');gridon;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个零');gridon;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个零');gridon;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个零');gridon;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);%矩形窗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),'.');gridon;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),'.');gridon;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),'.');gridon;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),'.');gridon;xlabel('Hz')title('改进升余弦窗')13、%频谱分析，平均clearall;%产生正弦加白噪声信号，求频谱，平均%产生两个正弦加白噪声；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,'.');gridon;xlabel('Hz')Y=zeros(1,N);u=rand(1,1000*N);u=u-mean(u);fori=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')1高通2带阻3低通4带通5切比雪夫带通1、1高通，双线性变换法，巴特沃斯（两种方法）%todesign高通high-passDFwiths=2/Ts[(z-1)/(z+1)]%-给出：通带下限频率，阻带上限频率，通带衰减，阻带衰减，采样频率%clearallwp=.8*pi;%通带下限频率ws=.44*pi;%阻带上限频率rp=3;%通带衰减rs=20;%阻带衰减Fs=2000;%采样频率%设计模拟滤波器%Firstlytofinishfrequencyprewarping;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)转换成H（s）,wap为低通的截止频率%模拟转换成数字%Note:z=(2/ts)(z-1)/(z+1);ts=1,thatis2fs=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个频率等分点。figure(1);plot(w,20*log10(abs(h)));grid;ylabel('High-passDF');%高通滤波器%直接设计高通滤波器，把上一种方法的buttord，buttap，lp2hp全包括进去%DirectlytodesignH(z)bybutter.m[n,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-passDF');1、2高通，双线性变换法，巴特沃斯和切比雪夫（两种方法）（1）%高通,双线性Z变换法，巴特沃斯,数字滤波器s=2/Ts[(z-1)/(z+1)]clearall%%给定参数fp=400;%通带下限频率fs=300;%阻带上限频率Fs=1000;%采样频率rp=3;%通带衰减rs=35;%阻带衰减%频率转化%%得到数字角频率wp=2*pi*fp/Fs;%得到数字角频率，2*pi对应Fs，fp对应的角频率ws=2*pi*fs/Fs;%得到数字角频率，2*pi对应Fs，fp对应的角频率%设计模拟滤波器%Firstlytofinishfrequencyprewarping;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);%模拟低通滤波器转换成高通滤波器%双线性变换法模拟转换成数字[bz,az]=bilinear(bt,at,Fs);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z）[h,w]=freqz(bz,az,256,Fs);%求离散系统频响特性Hphase=angle(h);%相位Hphase=unwrap(Hphase);%解卷绕，使相位在pi处不发生跳变，从而反应出真实的相位变化figure(1)subplot(211);plot(w,20*log10(abs(h)));grid;%对数幅频ylabel('幅频特性')subplot(212);plot(w,Hphase);gridon;ylabel('相频特性')%直接设计法%DirectlytodesignH(z)bybutter.m[n,wn]=buttord(wp/pi,ws/pi,rp,rs);[b,a]=butter(n,wp/pi,'high');[h1,w1]=freqz(b,a,256,Fs);Hphase=angle(h1);Hphase=unwrap(Hphase);h1=20*log10(abs(h1));figure(2)subplot(211);plot(w1,h1);grid;ylabel('High-passDF');title('直接设计法')subplot(212);plot(w1,Hphase);title('直接设计法')gridon;（2）%高通滤波器,双线性Z变换法，切比雪夫滤波器high-passDFwiths=2/Ts[(z-1)/(z+1)]clearall%%初始参数fp=400;%通带下限频率fs=300;%阻带上限频率Fs=1000;%采样频率rp=3;%通带衰减rs=35;%阻带衰减%%归一化频率，得到数字角频率wp=2*pi*fp/Fs;%得到数字角频率，2*pi对应Fs，fp对应的角频率ws=2*pi*fs/Fs;%得到数字角频率，2*pi对应Fs，fp对应的角频率%设计模拟滤波器%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp/2);%通带下限频率，模拟角频率，预畸变was=2*Fs*tan(ws/2);%阻带上限频率，模拟角频率，预畸变%模拟滤波器设计[n,wn]=cheb1ord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数n，切比雪夫滤波器阶数的选择%Note:'s'!%就这一句和巴特沃斯滤波器不同[z,p,k]=cheb1ap(n,rp);%设计模拟低通滤波器，极点，零点，增益[b,a]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型[bt,at]=lp2hp(b,a,wap);%模拟低通滤波器转换成高通滤波器%双线性变换法模拟转换成数字[bz,az]=bilinear(bt,at,Fs);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z）[h,w]=freqz(bz,az,256,Fs);%求离散系统频响特性Hphase=angle(h);%相位Hphase=unwrap(Hphase);%解卷绕，使相位在pi处不发生跳变，从而反应出真实的相位变化figure(1)subplot(211);plot(w,20*log10(abs(h)));grid;%对数幅频ylabel('幅频特性')subplot(212);plot(w,Hphase);ylabel('相频特性');gridon;%直接设计法%DirectlytodesignH(z)bycheby1[n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs);[bz1,az1]=cheby1(n,rp,wp/pi,'high');[h1,w1]=freqz(bz1,az1,256,Fs);Hphase=angle(h1);Hphase=unwrap(Hphase);h1=20*log10(abs(h1));figure(2)subplot(211);plot(w1,h1);title('直接设计法');gridylabel('High-passDF')subplot(212);plot(w1,Hphase);gridon;2、1带阻，巴特沃斯，双线性变换%带阻，陷波TodesignIIRButteworthbandstopDFbyanalog-lowpass,%给出：通带上下限频率，阻带上下限频率，抽样频率，通带衰减，阻带衰减clearall;%基本参数fp=[95105];%通带下限和上限频率,数字角频率fs=[99101];%阻带的上边和下边频率，数字角频率%wp=[.19*pi0.21*pi];ws=[.198*pi0.202*pi];Fs=1000;%抽样频率rp=3;%通带衰减，单位dBrs=14;%阻带衰减，单位dB%频率转换wp=fp*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率ws=fs*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率%wp=fp*2*pi*Fs;%数字角频率向模拟角频率转换%ws=fs*2*pi*Fs;%数字角频率向模拟角频率转换%设计模拟滤波器%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp./2);%预畸变,wp=(2/T)*tan(wp/200)=2*Fs*tan(wp*Fs/2)，模拟角频率was=2*Fs*tan(ws./2);%预畸变,wp=(2/T)*tan(wp/200)=2*Fs*tan(wp*Fs/2)，模拟角频率[n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数n%Note:'s'!%滤波器设计[z,p,k]=buttap(n);%设计模拟低通滤波器，极点，零点，增益[b,a]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型bw=wap(2)-wap(1);%通带带宽w0=sqrt(wap(1)*wap(2));%阻带中心频率[bt,at]=lp2bs(b,a,w0,bw);%从低通到带阻的转换,中心频率w0，带宽bw%%Note:z=(2/ts)(z-1)/(z+1);[bz1,az1]=bilinear(bt,at,Fs);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z），bz,az分别是H（Z）的分子分母多项式的系数[h,w]=freqz(bz1,az1,256,Fs);%求离散系统频响特性figure(1);plot(w,20*log10(abs(h)));title('对数幅频')%直接设计，双线性变换，巴特沃斯%DirectlytodesignH(z)bybutter.m[n,wn]=buttord(wp/pi,ws/pi,rp,rs);%要知道阶数N和3dB截止频率矢量wn，buttord()计算这些参数[b,a]=butter(n,wp/pi,'stop');%设计带阻滤波器[h1,w1]=freqz(b,a,256,Fs);%离散系统频响特性figure(2);subplot(211)plot(w1,20*log10(abs(h1)));title('对数幅频')%分贝subplot(212);plot(w1,abs(h1));title('幅值')%幅值2、2带阻，巴特沃斯，切比雪夫，双线性变换巴特沃斯%带阻，陷波，双线性变换法，巴特沃斯数字，带阻滤波器%clear;%初始参数fp=[4456];fs=[4753];rp=3;rs=50;Fs=1000;wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;%%求出阶次；[n,wn]=buttord(wp/pi,ws/pi,rp,rs);%直接设计Butterworth带阻滤波器；[b,a]=butter(n,wp/pi,'stop');[h,w]=freqz(b,a,256,Fs);Hphase=angle(h);Hphase=unwrap(Hphase);h=20*log10(abs(h));subplot(211);plot(w,h);gridon;ylabel('BandstopDF')xlabel('Hz')subplot(212);plot(w,Hphase);gridon;（2）%双线性变换法，切比雪夫I型，带阻数字滤波器%clearall;%初始参数f1=44;%通带下限截止频率f3=56;%通带上限截止频率fsl=47;%阻带下限频率fsh=53;%阻带上限频率rp=3;%通带衰减rs=35;%阻带衰减Fs=1000;%采样频率%%频率归一化wp1=2*pi*f1/Fs;%频率归一化，数字角频率，通带下限截止频率wp3=2*pi*f3/Fs;%通带上限截止频率wsl=2*pi*fsl/Fs;%阻带下限频率wsh=2*pi*fsh/Fs;%阻带上限频率wp=[wp1wp3];%通带截止频率ws=[wslwsh];%阻带截止频率%%直接设计切比雪夫滤波器；[n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs);%求取数字带通滤波器阶数n，切比雪夫I型滤波器阶数的选择%wp，ws为1X2向量，求出的是带通和带阻滤波器[bz1,az1]=cheby1(n,rp,wp/pi,'stop');%设计切比雪夫滤波器，stop代表带阻滤波器[h,w]=freqz(bz1,az1,256,Fs);%求离散系统频响特性Hphase=angle(h);%相位Hphase=unwrap(Hphase);%解卷绕h=20*log10(abs(h));%对数幅频subplot(211);plot(w,h);gridon;ylabel('幅频特性')subplot(212);plot(w,Hphase);gridon;ylabel('相频特性')3、1低通，双线性变换，巴特沃斯%低通，双线性变换，巴特沃斯%totestbuttord,lp2lp,bilinear;todesignLow-passDFwiths=(z-1)/(z+1)%通带上限频率，阻带下限频率，通带衰减，阻带衰减，采样频率clearall;%初始参数fp=100;%通带上限频率，数字角频率fs=300;%阻带下限频率，数字角频率Fs=1000;%采样频率rp=3;%通带允许最大衰减rs=20;%阻带应达到最小衰减%频率转换wp=2*pi*fp/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率ws=2*pi*fs/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率Fs=Fs/Fs;%letFs=1%Firstlytofinishfrequencyprewarping;wap=tan(wp/2);%预畸变,wp=2*Fs*tan(wp*Fs/2)，Fs=1，模拟角频率，数字转换成模拟was=tan(ws/2);%预畸变,wp=2*Fs*tan(wp*Fs/2)，Fs=1，模拟角频率，数字转换成模拟[n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数n%Note:'s'!%滤波器设计，双线性变换[z,p,k]=buttap(n);%设计模拟低通滤波器，极点，零点，增益[bp,ap]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型[bs,as]=lp2lp(bp,ap,wap);%从低通到低通的转换,wap为截止频率%Note:s=(2/Ts)(z-1)/(z+1);Ts=1,thatis2fs=1,fs=0.5;[bz,az]=bilinear(bs,as,Fs/2);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z），bz,az分别是H（Z）的分子分母多项式的系数[h,w]=freqz(bz,az,256,Fs*1000);%plot(w,abs(h));gridon;%直接设计，巴特沃斯，双线性变换wp=.2*pi;ws=.6*pi;Fs=1000;rp=3;rs=20;%Firstlytofinishfrequencyprewarping;[n,wn]=buttord(wp/pi,ws/pi,rp,rs);[bz,az]=butter(n,wp/pi);%[bz1,az1]=butter(n,wn);%[h,w]=freqz(bz,az,128,Fs);[h1,w1]=freqz(bz1,az1,128,Fs);plot(w,abs(h),w1,abs(h1),'g.');gridon;%todesignLow-passDFwiths=2/Ts[(z-1)/(z+1)]%clearall;wp=.2*pi;ws=.6*pi;Fs=1000;rp=3;rs=20;%%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp/2);was=2*Fs*tan(ws/2);[n,wn]=buttord(wap,was,rp,rs,'s');%Note:'s'![z,p,k]=buttap(n);[bp,ap]=zp2tf(z,p,k);[bs,as]=lp2lp(bp,ap,wap);w1=[0:499]*2*pi;h1=freqs(bs,as,w1);%Note:z=(2/ts)(z-1)/(z+1);[bz,az]=bilinear(bs,as,Fs);[h2,w2]=freqz(bz,az,500,Fs);plot(w1/2/pi,abs(h1),w2,abs(h2),'k');gridon;3、2低通，双线性变换，巴特沃斯（三种方法）%totestbuttord,lp2lp,bilinear;todesignLow-passDFwiths=2/Ts[(z-1)/(z+1)]%letFs=1,s=(z-1)/(z+1)%低通，双线性Z变换法，巴特沃斯，数字滤波器%给定参数：通带衰减，阻带衰减，通带截止频率，阻带截止频率，采样频率clearall;%给定初始参数Fs=100000;%抽样频率rp=3;%通带衰减rs=30;%阻带衰减%频率归一化wp=0.2*pi;%通带截止频率ws=0.5*pi;%阻带截止频率Fs=Fs/Fs;%letFs=1%求取模拟低通滤波器阶数n%Firstlytofinishfrequencyprewarping;wap=tan(wp/2);%通带截止频率，模拟角频率，预畸变was=tan(ws/2);%阻带截止频率，模拟角频率，预畸变%模拟滤波器设计[n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟低通滤波器阶数n%Note:'s'!%设计模拟低通滤波器[z,p,k]=buttap(n);%设计模拟低通滤波器[bp,ap]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型[bs,as]=lp2lp(bp,ap,wap);%模拟低通滤波器转换成低通滤波器%Note:s=(2/Ts)(z-1)/(z+1);Ts=1,thatis2fs=1,fs=0.5;%双线性变换法，数字滤波器%由模拟滤波器H（s）得到数字滤波器H（Z）[bz,az]=bilinear(bs,as,Fs/2);%由模拟滤波器H（s）得到数字滤波器H（Z），实现双线性变换[h,w]=freqz(bz,az,256,Fs*100000);%数字滤波器H(Z)的频率响应plot(w,abs(h));gridon;%todesignLow-passDFwiths=2/Ts[(z-1)/(z+1)]%clearall;wp=.2*pi;%通带截止频率ws=.5*pi;%阻带截止频率Fs=100000;%抽样频率rp=3;%通带衰减rs=30;%阻带衰减%频率归一化%Firstlytofinishfrequencyprewarping;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);%设计模拟低通滤波器[bp,ap]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型[bs,as]=lp2lp(bp,ap,wap);%模拟低通滤波器转换成低通滤波器%Note:z=(2/ts)(z-1)/(z+1);[bz,az]=bilinear(bs,as,Fs);%由模拟滤波器H（s）得到数字滤波器H（Z），实现双线性变换[h2,w2]=freqz(bz,az,500,Fs);%数字滤波器H(Z)的频率响应%[h2,w2]=freqz(bz,az,256,Fs);%数字滤波器H(Z)的频率响应plot(w2,abs(h2));gridon;%直接设计法wp=.2*pi;ws=.5*pi;Fs=100000;rp=3;rs=30;%Firstlytofinishfrequencyprewarping;[n,wn]=buttord(wp/pi,ws/pi,rp,rs);[bz,az]=butter(n,wp/pi);%[bz1,az1]=butter(n,wn);%[h,w]=freqz(bz,az,128,Fs);[h1,w1]=freqz(bz1,az1,128,Fs);plot(w,abs(h),w1,abs(h1),'g.');gridon;4、1带通，双线性变换，巴特沃斯，直接和间接设计%带通，双线性变换，巴特沃斯，todesignaButterworthBandpassdigitalfilter.%给出：通带截止频率，阻带截止频率，通带衰减，阻带衰减，采样频率%--clearall;%初始参数fp=[300400];%通带范围fs=[200500];%阻带边缘频率rp=3;%通带带边衰减rs=18;%阻带带边衰减Fs=2000;%采样频率%频率转换wp=fp*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率ws=fs*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率%设计模拟滤波器%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp./2);%预畸变,wp=(2/T)*tan(wp/200)=2*Fs*tan(wp*Fs/2)，模拟角频率was=2*Fs*tan(ws./2);%预畸变,wp=(2/T)*tan(wp/200)=2*Fs*tan(wp*Fs/2)，模拟角频率%模拟滤波器设计[n,wn]=buttord(wap,was,rp,rs,'s');%求取模拟带通滤波器阶数n%Note:'s'![z,p,k]=buttap(n);%设计模拟带通滤波器，极点，零点，增益[bp,ap]=zp2tf(z,p,k);%零－极点增益模型转换为传递函数模型%bw=wap(2)-wap(1);%通带带宽w0=sqrt(wap(1)*wap(2));%通带中心频率[bs,as]=lp2bp(bp,ap,w0,bw);%从低通到带通的转换,中心频率w0，带宽bw%求模拟低通滤波器的频响特性[h1,w1]=freqs(bp,ap);%求模拟低通滤波器的频响特性,不出图形figure(1);subplot(311);plot(w1,abs(h1));grid;%幅频特性，出图形ylabel('模拟低通G(p)');%求模拟带通滤波器的频响特性w2=[0:Fs/2-1]*2*pi;h2=freqs(bs,as,w2);%求模拟带通滤波器的频响特性,不出图形subplot(312);plot(w2,abs(h2));grid;%幅频特性，出图形ylabel('模拟带通幅值G(p)');subplot(313);plot(w2/2/pi,20*log10(abs(h2)));ylabel('模拟带通对数幅频G(p)');%数字滤波器设计，双线性变换法%Note:z=(2/Ts)(z-1)/(z+1);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z）[bz1,az1]=bilinear(bs,as,Fs);%实现双线性变换，由模拟滤波器H（s）得到数字滤波器H（Z）[h3,w3]=freqz(bz1,az1,1000,Fs);%求离散系统频响特性figure(2);plot(w2/2/pi,20*log10(abs(h2)),w3,20*log10(abs(h3)));grid;ylabel('BandpassAFandDF');xlabel('Hz');%直接设计fp=[300400];fs=[200500];rp=3;rs=18;Fs=2000;wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;%不需要再预畸变%求出阶次；[n,wn]=buttord(wp/pi,ws/pi,rp,rs);%再设计Butterworth带通滤波器；[b,a]=butter(n,wp/pi)[h,w]=freqz(b,a,256,Fs);h=20*log10(abs(h));plot(w,h);grid;ylabel('BandpassDF')xlabel('Hz')4、2带通，双线性变换，巴特沃斯，切比雪夫，直接设计（1）%带通,双线性变换，巴特沃斯数字滤波器%clear;%初始参数fp=[300400];%通带范围fs=[200500];%阻带边缘频率rp=3;%通带带边衰减rs=40;%阻带带边衰减Fs=2000;%采样频率%频率归一化wp=fp*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率ws=fs*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率%设计巴特沃斯数字滤波器%求出阶次；%直接设计法[n,wn]=buttord(wp/pi,ws/pi,rp,rs);%求取模拟带通滤波器阶数n%再设计Butterworth数字带通滤波器；[b,a]=butter(n,wp/pi);%设计带通数字滤波器，已经包含bilinear双线性变换[h,w]=freqz(b,a,256,Fs);%求数字带通滤波器的频响特性,不出图形%相位和幅频特性Hphase=angle(h);%相位Hphase=unwrap(Hphase);%解卷绕，使相位在pi处不发生跳变，从而反应出真实的相位变化h=20*log10(abs(h));%对数幅频subplot(211);plot(w,h);gridon;ylabel('幅频特性');title('直接设计法')xlabel('Hz')subplot(212);plot(w,Hphase);ylabel('相频特性');gridon;（2）%带通,双线性变换，切比雪夫I型数字滤波器%clearall;%初始参数f1=300;%通带下限频率f3=400;%通带上限频率fsl=200;%下阻带的上限频率fsh=500;%上阻带的下限频率rp=3;%通带衰减rs=40;%阻带衰减Fs=2000;%采样频率%%归一化频率wp1=2*pi*f1/Fs;%归一化角频率wp3=2*pi*f3/Fs;%归一化角频率wsl=2*pi*fsl/Fs;%归一化角频率wsh=2*pi*fsh/Fs;%归一化角频率wp=[wp1wp3];ws=[wslwsh];%%直接设计切比雪夫滤波器[n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs);%求取数字带通滤波器阶数n，切比雪夫I型滤波器阶数的选择[bz1,az1]=cheby1(n,rp,wp/pi);%设计切比雪夫滤波器[h,w]=freqz(bz1,az1,256,Fs);%求离散系统频响特性Hphase=angle(h);%相位Hphase=unwrap(Hphase);%解卷绕，使相位在pi处不发生跳变，从而反应出真实的相位变化%h=20*log10(abs(h));h=abs(h);subplot(211);plot(w,h);gridon;%幅频ylabel('幅频')subplot(212);plot(w,Hphase);gridon;%相频ylabel('相频')5、切比雪夫滤波器，双线性变换，带通%切比雪夫I型，带通，双线性变换todesignaChebyshev-IBandpassDF.%clearall;%给定参数，频率值f1=300;f3=500;fsl=200;fsh=600;rp=0.1;rs=30;Fs=2000;%%频率转换wp1=2*pi*f1/Fs;%频率归一化wp3=2*pi*f3/Fs;wsl=2*pi*fsl/Fs;wsh=2*pi*fsh/Fs;wp=[wp1wp3];ws=[wslwsh];%%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp./2);%预畸变was=2*Fs*tan(ws./2);%模拟滤波器设计[n,wn]=cheb1ord(wap,was,rp,rs,'s');%Note:'s'![z,p,k]=cheb1ap(n,rp);[bp,ap]=zp2tf(z,p,k);bw=wap(2)-wap(1);w0=sqrt(wap(1)*wap(2));[bs,as]=lp2bp(bp,ap,w0,bw);%%双线性变换，数字滤波器设计%Note:z=(2/Ts)(z-1)/(z+1);[bz1,az1]=bilinear(bs,as,Fs)[h,w]=freqz(bz1,az1,256,Fs);plot(w,20*log10(abs(h)));gridon;%直接设计法%初始参数f1=300;f3=500;fsl=200;fsh=600;rp=0.1;rs=30;Fs=2000;%%归一化频率wp1=2*pi*f1/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsl/Fs;wsh=2*pi*fsh/Fs;wp=[wp1wp3];ws=[wslwsh];%%设计切比雪夫滤波器，双线性变换[n,wn]=cheb1ord(ws/pi,wp