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

%---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)/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)

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;

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); %绘制x

stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on; subplot(232); %绘制h

stem(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、

%求序列的自相关函数 clear N=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)/N

nx=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两个信号 %产生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=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;%时域相乘，频域卷积 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*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).*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; end

Y=Y/1000;

subplot(212);plot(f(1:N/2),Y(1:N/2));title('正弦加白噪声后与w时域相乘的FFT 1000次平均') xlabel('Hz')

x(t) w(t)10.50-0.5-10246x(t)频谱3020810121000510152025Hz3035404550

正弦加白噪声后与w时域相乘10.50-0.5024681012正弦加白噪声后与w时域相乘的FFT10500510152025Hz3035404550

正弦加白噪声的FFT 1000次平均1816141210052530354045Hz正弦加白噪声后与w时域相乘的FFT 1000次平均10152050654320510152025Hz3035404550

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*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)420-2-40100200300400500600x(t)加窗之后的傅里叶变换0.80.60.40.200510152025Hz3035404550 %分段函数----------------------------------------------- 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)为分段函数10.50-0.5-10100200300x(t)频谱分析0.20.150.10.0500510152025Hz3035404550400500600

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/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(4*45:4*60));title('N=1024');grid on;

xlabel('Hz')

N=256N=10241506001004005020001.522.501.52HzN=2561501005000246Hz

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; 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')

2.5

% 在数据末补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')

不补零82005Hz补7N个零1082005Hz补29N个零10补N个零101055005Hz10005Hz10

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: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)100-105000102030fft(x)405060700-0.5100-10-0.4-0.3-0.2-0.10x(n)实部0.10.20.30.40.5010203040506070

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*h z=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 卷积')

x(t)105020100201000102030405060708090100fftfilt 叠接相加法01020304050conv 卷积607080901000102030405060708090100

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),'.');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')

未补零8200102030405060708090Hz补N个零82001020304050Hz60708090100

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)

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)

stem(f,X(1:N/2),'.');grid on; xlabel('Hz')

1.5x 10-1410.500204060Hz820050100Hz

9、

%周期延拓 clear all;

820020406080Hz

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*(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.50-0.5-10246810121501005000510152025Hz3035404550

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)');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)100005000000.51HzFFT x(n)1.522.51000050000-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5

11、

clear;

%补零对频谱分析的影响

f1=10.8;f2=11.75;f3=12.55;fs=40;w=2*pi/fs; N=;

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个零403020100051015Hz补7N个零2040302010005补3N个零1015Hz补15N个零204030201000510Hz15204030201000510Hz1520

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),'.');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('改进升余弦窗')

矩形窗10.5005101520Hz平顶窗1.510.5005101520Hz

13、

%频谱分析，平均

升余弦窗10.5005101520Hz改进升余弦窗10.5005101520Hz

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: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);

subplot(212);plot(f(1:N/2),Y(1:N/2));title('平均后') axis([0,50,0,60]); xlabel('Hz')

平均前60402000510152025Hz平均后303540455060402000510152025Hz3035404550

1高通 2带阻 3低通 4带通 5切比雪夫 带通

1、1高通，双线性变换法，巴特沃斯（两种方法）

% to design 高通high-pass DF with s=2/Ts[(z-1)/(z+1)]

%-给出：通带下限频率，阻带上限频率，通带衰减，阻带衰减，采样频率----------------------------------------------- %---------------------- clear all

wp=.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是模拟低通滤波器的阶次,巴特沃斯滤波器阶数的选择，（最小阶数的选择）

% 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,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个频率等分点。 figure(1);

plot(w,20*log10(abs(h)));grid;

ylabel(' High-pass DF');%高通滤波器

%直接设计高通滤波器，把上一种方法的buttord，buttap，lp2hp全包括进去 % Directly to design H(z) by butter.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-pass DF');

0-50-100 High-pass DF-150-200-250-300-35000.050.10.150.20.250.30.350.40.450.5

1、2高通，双线性变换法，巴特沃斯和切比雪夫（两种方法）

（1）% 高通,双线性Z变换法，巴特沃斯,数字滤波器 s=2/Ts[(z-1)/(z+1)] clear all %

%给定参数

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对应的角频率 %设计模拟滤波器

% 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，巴特沃斯滤波器阶数的选择 % 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);grid on; ylabel(' 相频特性')

%直接设计法

% Directly to design H(z) by butter.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-pass DF');title('直接设计法') subplot(212);

plot(w1,Hphase);title('直接设计法') grid on;

直接设计法200 High-pass DF0-200-400050100150200250300350400450500直接设计法50-5-10050100150200250300350400450500

（2）

% 高通滤波器,双线性Z变换法，切比雪夫滤波器high-pass DF with s=2/Ts[(z-1)/(z+1)] clear all %

%初始参数

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对应的角频率 %设计模拟滤波器

% Firstly to finish frequency prewarping;

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('相频特性');grid on;

%直接设计法

% Directly to design H(z) by cheby1 [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('直接设计法');grid ylabel(' High-pass DF')

subplot(212);plot(w1,Hphase);grid on;

直接设计法0 High-pass DF-100-200-300-4000501001502002503003504004505000-2-4-6-8050100150200250300350400450500

2、1带阻，巴特沃斯，双线性变换

% 带阻，陷波 To design IIR Butteworth bandstop DF by analog-lowpass,

% 给出：通带上下限频率，阻带上下限频率，抽样频率，通带衰减，阻带衰减--------------------------------- clear all;

%基本参数

fp=[95 105];%通带下限和上限频率,数字角频率 fs=[99 101];%阻带的上边和下边频率，数字角频率 %wp=[.19*pi 0.21*pi];ws=[.198*pi 0.202*pi]; Fs=1000;%抽样频率

rp=3;%通带衰减，单位dB rs=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;%数字角频率向模拟角频率转换 %设计模拟滤波器

% Firstly to finish frequency prewarping;

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('对数幅频')

%直接设计，双线性变换，巴特沃斯 % Directly to design H(z) by butter.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('幅值')%幅值

对数幅频对数幅频102000-20-40-10-60-20050100150200250幅值3003504004505001.5-301-400.5-500501001502002503003504004505000050100150200250300350400450500

2、2带阻，巴特沃斯，切比雪夫，双线性变换

（1）巴特沃斯

%带阻，陷波，双线性变换法，巴特沃斯数字，带阻滤波器 % -------------------------------------------------------------------------

clear;

%初始参数

fp=[44 56];fs=[47 53]; 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);grid on; ylabel('Bandstop DF') xlabel(' Hz') subplot(212);

plot(w,Hphase);grid on;

（2）

% 双线性变换法，切比雪夫I型，带阻数字滤波器

%---------------------------------------------------------------------------- clear all;

%初始参数

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=[wp1 wp3];%通带截止频率 ws=[wsl wsh];%阻带截止频率 %

% 直接设计切比雪夫滤波器；

[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);grid on; ylabel('幅频特性') subplot(212);

plot(w,Hphase);grid on; ylabel('相频特性')

0-20幅频特性-40-60-80-1000501001502002503003504004505000相频特性-5-10-15050100150200250300350400450500

3、1低通，双线性变换，巴特沃斯

% 低通，双线性变换，巴特沃斯

%to test buttord,lp2lp,bilinear ;to design Low-pass DF with s=(z-1)/(z+1) %通带上限频率，阻带下限频率，通带衰减，阻带衰减，采样频率--------------------------------------

clear all;

%初始参数

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; % let Fs=1

% Firstly to finish frequency prewarping ;

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,that is 2fs=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));grid on;

%直接设计，巴特沃斯，双线性变换 wp=.2*pi; ws=.6*pi; Fs=1000; rp=3;rs=20;

% Firstly to finish frequency prewarping; [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.');grid on;

1.41.210.80.60.40.20

% to design Low-pass DF with s=2/Ts[(z-1)/(z+1)]

%----------------------------------------------------------------------------- clear all;

wp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20; %

% 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'); % 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');grid on;

0501001502002503003504004505003、2低通，双线性变换，巴特沃斯（三种方法）

% to test buttord,lp2lp,bilinear ;to design Low-pass DF with s=2/Ts[(z-1)/(z+1)] %let Fs=1 ,s=(z-1)/(z+1)

%低通，双线性Z变换法，巴特沃斯，数字滤波器

%给定参数：通带衰减，阻带衰减，通带截止频率，阻带截止频率，采样频率------------------------------------ clear all;

%给定初始参数

Fs=100000;%抽样频率

rp=3;%通带衰减 rs=30;%阻带衰减 %频率归一化

wp=0.2*pi;%通带截止频率 ws=0.5*pi;%阻带截止频率 Fs=Fs/Fs; % let Fs=1 %求取模拟低通滤波器阶数n

% Firstly to finish frequency prewarping ;

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,that is 2fs=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));grid on;

% to design Low-pass DF with s=2/Ts[(z-1)/(z+1)]

%----------------------------------------------------------------------------- clear all;

wp=.2*pi;%通带截止频率 ws=.5*pi;%阻带截止频率 Fs=100000;%抽样频率 rp=3;%通带衰减 rs=30;%阻带衰减 %频率归一化

% 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 % 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));grid on;

%直接设计法 wp=.2*pi; ws=.5*pi; Fs=100000; rp=3;rs=30;

% Firstly to finish frequency prewarping; [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.');grid on;

s=(z-1)/(z+1)1.41.210.80.60.40.2000.511.522.533.544.5x 1054

直接设计法s=2/Ts[(z-1)/(z+1)]10.90.80.70.60.50.40.30.20.1000.511.522.533.544.5x 105410.90.80.70.60.50.40.30.20.1000.511.522.533.544.5x 1054

4、1带通，双线性变换，巴特沃斯，直接和间接设计

% 带通，双线性变换，巴特沃斯，to design a Butterworth Bandpass digital filter. % 给出：通带截止频率，阻带截止频率，通带衰减，阻带衰减，采样频率------------------------------------ % -- clear all;

%初始参数

fp=[300 400];%通带范围 fs=[200 500];%阻带边缘频率 rp=3;%通带带边衰减 rs=18;%阻带带边衰减 Fs=2000;%采样频率

%频率转换

wp=fp*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率 ws=fs*2*pi/Fs;%得到角频率，2*pi对应Fs，fp对应的角频率，数字角频率 %设计模拟滤波器

% Firstly to finish frequency prewarping;

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('Bandpass AF and DF'); xlabel(' Hz');

%直接设计

fp=[300 400];fs=[200 500]; 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('Bandpass DF') xlabel(' Hz')

1 模拟低通 G(p)0.5010123456710 模拟带通幅值 G(p)0.500-100-200-300-40001000200030004000500060007000 模拟带通对数幅频 G(p)010020030040050060070080090010000-50-100Bandpass AF and DF-150-200-250-300-3500100200300400500 Hz6007008009001000

0-20直接设计Bandpass DF-40-60-80-100-1200100200300400500 Hz6007008009001000

4、2带通，双线性变换，巴特沃斯，切比雪夫，直接设计

（1）% 带通,双线性变换，巴特沃斯数字滤波器

% ------------------------------------------------------------------------- clear;

%初始参数

fp=[300 400];%通带范围 fs=[200 500];%阻带边缘频率 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);grid on;

ylabel('幅频特性');title('直接设计法') xlabel(' Hz') subplot(212); plot(w,Hphase);

ylabel('相频特性');grid on;

（2）

% 带通,双线性变换，切比雪夫I型数字滤波器

%---------------------------------------------------------------------------- clear all;

%初始参数

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=[wp1 wp3]; ws=[wsl wsh]; %

% 直接设计切比雪夫滤波器

[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);grid on;%幅频 ylabel('幅频') subplot(212);

plot(w,Hphase);grid on;%相频 ylabel('相频')

切比雪夫10.8幅频0.60.40.20010020030040050060070080090010000-5相频-10-1501002003004005006007008009001000

5、切比雪夫滤波器，双线性变换，带通

% 切比雪夫I型，带通，双线性变换to design a Chebyshev-I Bandpass DF. %---------------------------------------------------------------------------- clear all;

%给定参数，频率值 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=[wp1 wp3];ws=[wsl wsh]; %

% Firstly to finish frequency prewarping; 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)));grid on;

%直接设计法

%初始参数 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=[wp1 wp3];ws=[wsl wsh]; %

% 设计切比雪夫滤波器，双线性变换 [n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi) [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h)); plot(w,h);grid on;

0-50-100-150-200-250-30001002003004005006007008009001000

直接设计0-50-100-150-200-250-300010020030040050060070080090010001低通 2多阻带

1、切比雪夫最佳一致逼近，低通滤波器

% 切比雪夫最佳一致逼近，低通滤波器

%------------------------------------------------------------------------- clear all;

%给定参数 f=[0 0.6 0.7 1];

% 给定频率轴分点；通带边缘频率0.6*pi,阻带边缘频率0.7*pi A=[1 1 0 0];

% 给定在这些频率分点上理想的幅频响应； %代表低通 weigh=[10 1];

% 给定在这些频率分点上的加权；

%设计滤波器

b=remez(32,f,A,weigh);

%32阶，f频率向量，A幅频响应，weight权重,b是设计出的滤波器的系数 % 设计出切比雪夫最佳一致逼近滤波器；

%频率响应

[h,w]=freqz(b,1,256,1);

%系统频率响应b/a=h，a=1

%256为频率轴的分点，指定计算频率的范围从0到1， h=abs(h);

h=20*log10(h); figure(1)

stem(b,'.');title('系数');grid; figure(2)

plot(w,h);title('幅频响应');grid;

方法2

% to test remezord.m,

%-------------------------------------------------------------------------- clear all;

% 先求滤波器的阶次； f=[.6 .7]; A=[1 0]; rp=0.4916; rs=44.7;

dev=[10^(rp/20)-1 10^(-rs/20)];%通带和阻带的偏差 %

[n,f0,A0,w]=remezord(f,A,dev); %

if rem(n,2) n=n+1; end

% 再用切比雪夫最佳一致逼近设计线性相位FIR滤波器； b=remez(n,f0,A0,w);

[h,w]=freqz(b,1,256,1); h=abs(h);

h=20*log10(h); figure(1)

stem(b,'.');grid; figure(2)

plot(w,h);grid;

2、切比雪夫最佳一致逼近，高通滤波器

（1）方法1

% to use remez.m and to design high-pass FIR filter;

%高通------------------------------------------------------------------------- clear all;

f=[0 0.5 0.6 1];

% 给定频率轴分点； A=[0 0 1 1];%和低通不同

% 给定在这些频率分点上理想的幅频响应； weigh=[1 1];

% 给定在这些频率分点上的加权；

b=remez(54,f,A,weigh);

% 设计出切比雪夫最佳一致逼近滤波器； %

[h,w]=freqz(b,1,256,1); hr=abs(h);

hr=20*log10(hr);

hphase=angle(h);hphase=unwrap(hphase); figure(1)

stem(b,'.');grid; figure(2)

plot(w,hr);grid; figure(3)

plot(w,hphase);grid;

（2）方法2

% to test remezord.m and use remez.m and to design high-pass FIR filter; %高通------------------------------------------------------------------------- clear all;

% 先求滤波器的阶次； f=[.5 .6];

A=[0 1];%和低通不同 rp=0.4916; rs=44.7;

dev=[10^(-rs/20) 10^(rp/20)-1]; %

[n,f0,A0,w]=remezord(f,A,dev) %

if rem(n,2) n=n+1; end

% 再用切比雪夫最佳一致逼近设计线性相位FIR滤波器； b=remez(n,f0,A0,w);

[h,w]=freqz(b,1,256,1); hr=abs(h);

hr=20*log10(hr);

hphase=angle(h);hphase=unwrap(hphase); figure(1)

stem(b,'.');grid; figure(2)

plot(w,hr);grid; figure(3)

plot(w,hphase);grid;

2、切比雪夫最佳一致逼近，带阻滤波器，陷波

%多阻带滤波器，陷波器，滤掉50,100,150Hz的干扰---------------------------------- clear all;

% 用切比雪夫最佳一致逼近设计线性相位多带FIR滤波器； f=[0 .14 .18 .22 .26 .34 .38 .42 .46 .54 .58 .62 .66 1];

%归一化中心频率50归为0.1,100归为0.2,150归为0.3，三个频率都乘以2成为：0.2,0.4,0.6三个阻带

%[0 0.14]通带，[0.18 0.22]代表0.2阻带，[0.26 0.34]代表通带，

%[0.38 0.42]代表0.4阻带，[0.46 0.54]代表通带，[0.58 0.62]代表0.6阻带 %[0.66 1]代表通带

A=[1 1 0 0 1 1 0 0 1 1 0 0 1 1]; weigh=[8 1 8 1 8 1 8];

b=remez(,f,A,weigh);%阶，返回滤波器系数b %

[h,w]=freqz(b,1,256,1);

%系统频率响应b/a=h，a=1

%256为频率轴的分点，指定计算频率的范围从0到1， hr=abs(h); h=abs(h);

h=20*log10(h); figure(1)

stem(b,'.');grid; figure(2)

plot(w,h);grid;

方法2

% to test remezord.m,

%---------------------------------------------------------------------------- clear all;

% 先求滤波器的阶次；

f=[.14 .18 .22 .26 .34 .38 .42 .46 .54 .58 .62 .66]; A=[1 0 1 0 1 0 1];

rp=0.153;%四个通带的衰减 rs=16.92;%三个阻带的衰减

devp=10^(rp/20)-1;%通带的偏差 devs=10^(-rs/20);%阻带的偏差

dev=[devp devs devp devs devp devs devp]; [n,f0,A0,w]=remezord(f,A,dev)

%再用切比雪夫最佳一致逼近设计多带线性相位FIR滤波器； b=remez(n+1,f0,A0,w);

[h,w]=freqz(b,1,256,1); hr=abs(h); h=abs(h);

h=20*log10(h); figure(1)

stem(b,'.');grid;

figure(2)

plot(w,h);grid; 3、