《数字信号处理相关matlab文件整理》

hhppan

已于 2023-11-04 21:24:13 修改

阅读量1k

点赞数 3

分类专栏：数字信号处理文章标签：信号处理 matlab

于 2023-10-31 21:52:59 首次发布

本文链接：https://blog.csdn.net/weixin_50889820/article/details/134149620

版权

数字信号处理专栏收录该内容

1 篇文章 0 订阅

订阅专栏

第一章离散时间信号与离散时间系统

1、rand函数

rand函数可以用来产生均值为0.5，幅度在0-1之间均匀分布的伪随机数，一般用它来近似均匀分布的白噪声信号u(n)。白噪声信号频谱幅度在整个频率范围内都一样，且功率等于方差。

例1、产生均匀分布白噪声的代码以及波形（结果如图1.1所示）

clear;
N=50000;%u(n)的长度
u=rand(1,N);%调用rand，得到均匀分布的伪随机数u(n)
u_mean=mean(u);%用matlab中的m文件求u(n)的均值
power_u=var(u);%用matlab中的m文件求u(n)的方差，方差即为功率
subplot(211);%绘图空间划分成两行一列，该图在第一行
plot(u(1:100));%绘制u(n)图像
grid on; %加网格
ylabel('(u(n)');%y轴标记
subplot(212);%该图在第二行
histogram(u,50);%绘制直方图，50表示绘制50个等距bin的数量
grid on;
ylabel('histogram of u(n)');

图1.1 均匀分布的白噪声信号

例2、产生一个均匀分布、均值为0、功率为0.01的白噪声信号u(n)。（结果如图1.2所示）

由上例可知，rand(N)给出的u(n)均值为0.5，功率约为1/12，

均值变为0可以通过u(n)-均值来实现

功率变为p=0.01，则需要使u(n)=a*u(n),原来功率为p1其中a=sqrt(p/p1)

                                       -------用功率表达式即可证明
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;
%试验u1(n)的功率是否满足要求，dot是Matlab的m文件，用来实现两个向量的内积，该句等效于var(u1)
plot(u1(1:100));
grid on;

图1.2 均匀分布、均值为0、功率为0.01的白噪声信号

2、randn函数

该文件可以产生均值为0、方差为1、服从高斯分布的白噪声信号u(n)，使用方法与rand函数相似。

例。产生均值为0、方差为1、服从高斯分布的白噪声信号u(n)。（结果如图1.3所示）

%randn可以用来产生均值为0，方差为1的服从正态分布的白噪声信号u(n)
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:100));
ylabel('u(n)');
subplot(212);
histogram(u,50);
ylabel('u(n)的直方图');

图1.3 randn函数产生的白噪声

3、sinc函数

本文件可用来产生sinc函数，对于连续信号，sinc函数定义

t代表时间，对于离散信号，相应的sinc函数定义为：

例、产生一sinc信号并显示其波形。（结果如图1.4所示）

clear;
n=200;%时间轴分点数
stept=4*pi/n;%时间轴的范围是从-2Π-2Π，stept是步长
t=-2*pi:stept:2*pi;
y=sinc(t);
plot(t,y,t,zeros(size(t)));%后面两个参数代表在y=0的位置画一条（-2Π，2Π）的直线
ylabel('sinc(t)');
xlabel('t');
grid on;

图1.4 sinc函数图像

4、chirp函数

本文件可以产生一种叫做chirp的信号，该信号基本形式为：

文件的调用格式为，其中T是横轴的时间范围向量，F0是起始频率，F1是在T1时刻所具有的频率。缺省时，F0=0,T1=1,F1=100

例、产生一chirp信号，初始频率为0，结尾频率为125。（结果如图1.5所示）

t=0:0.001:1;
y=chirp(t,0,1,125);
plot(t,y);
ylabel('x(t)');
xlabel('t');

图1.5 chirp信号

5、conv函数

该文件用来实现两个离散序列的线性卷积，调用格式为：

若X长度为N，h长度为M，则y长度为N+M-1（可以使用图像直观地证明）。

例、令x(n)={1,2,3,4,5},h(n)={6,2,3,6,4,2},y(n)=x(n)*h(n),画出y(n)。（如图1.6所示）

N=5;M=6;L=N+M-1;
x=[1,2,3,4,5];
nx=0:N-1;
h=[6,2,3,6,4,2];
nh=0:M-1;
y=conv(x,h);
ny=0:L-1;
%stem(Y) 将数据序列Y从x轴到数据值按照茎状形式画出，以圆圈终止。如果Y是一个矩阵，则将其每一列按照分隔方式画出。
%stem(X,Y)在X的指定点处画出数据序列Y. 
subplot(231);stem(nx,x,'.');xlabel('x(n)');grid on;
subplot(232);stem(nh,h,'.');xlabel('h(n)');grid on;
subplot(233);stem(ny,y,'.');xlabel('y(n)');grid on;

图1.6 卷积结果

6、xcorr函数

该文件用来求两个信号的互相关或者一个信号的自相关，调用格式为：

格式（1）是求两个序列x,y的互相关，若x,y长度都为N，则 $r_{xy}$ 的长度为2N-1，如果x,y长度不相等，短的一方补零。

格式（2）是求序列x的自相关，Mlag表示rx的单边长度，总的长度为，Mlag+1，flag是定标标志，flag=biased表示是有偏估计，需要将rx(m)都除以N，若flag=unbiased，则表示是无偏估计，需将rx(m)都除以（N-abs(m))。

例、求白噪声和正弦信号自相关函数。（结果如图1.7所示）

clear;
N=500;p1=1;p2=0.1;f=1/8;Mlag=60;
u=randn(1,N);
u1=u*sqrt(p2);
n=[0:N-1];
s=sin(2*pi*f*n);
x1=u(1:N)+s;
rxl=xcorr(x1,Mlag,'biased');
u2=u*sqrt(p2);x2=u2(1:N)+s;
rx2=xcorr(x2,Mlag,'biased');
subplot(211);
plot(x1);
subplot(212);
plot(rx2);

图1.7 自相关函数

第二章 Z变换及离散时间系统分析

1、filter.m文件

例、求该系统的阶跃响应，即求系统对阶跃输入的输出。

 x=ones(100);%x(n)=1,即产生阶跃信号
 t=1:100;
 b=[0.001836,0.007344,0.011016,0.007374,0.001836];
 a=[1,-3.0544,3.8291,-2.2925,0.55075];
 y=filter(b,a,x);
 plot(t,x,'g.',t,y,'k-');%实现将x(n)和y(n)画在同一个图上
 xlabel('x');
 ylabel('x(n),y(n)');

图2.1 系统的阶跃响应

2、impz.m文件

例、求该系统的单位抽样响应h(n)。

b=[0.001836,0.007344,0.011016,0.007374,0.001836];
a=[1,-3.0544,3.8291,-2.2925,0.55075];
[h,t]=impz(b,a,40);
stem(t,h,'.');
grid on;
xlabel('n');
ylabel('h(n)');

图2.2 系统的单位抽样响应

3、freqz.m文件

例、求该系统的频率响应，即幅频和相频响应。

b=[0.001836,0.007344,0.011016,0.007374,0.001836];
a=[1,-3.0544,3.8291,-2.2925,0.55075];
[H,w]=freqz(b,a,256,'whole',1);
%256代表频率轴的分点数，1是抽样频率 
Hr=abs(H);%求出H的模，即幅频响应
Hphase=angle(H);
Hphase=unwrap(Hphase);%求相频响应并解卷绕
subplot(211);
plot(w,Hr);ylabel('幅频');xlabel('w/2pi');
subplot(212);
plot(w,Hphase);ylabel('相位');xlabel('w/2pi')

图2.3 频率响应

4、zplane.m文件

本文件可以用来显示离散系统的极零图，调用格式为zplane(z,p)或zplane(b,z)，前者是在已知系统零点的列向量z和极点的列向量p的情况下画出极零图，后者是在仅已知B(z)和A(z)的情况下画出极零图。

例、显示该系统的极零图以及下面FIR系统的极零图

b=[0.001836,0.007344,0.011016,0.007374,0.001836];
a=[1,-3.0544,3.8291,-2.2925,0.55075];
subplot(221);zplane(b,a);
b=[1,-1.7,1.53,-0.684];a=1;%没有分母，所以a=1
subplot(222);zplane(b,a);

图2.4 不同系统极零点图

5、residuez.m文件

clear;
b=[1.7,-1.69,.39];%形成分子多项式向量
a=[1,-1.7,0.8,-.1];
[r,p,k]=residuez(b,a);%做z变换，求出向量r,p,k
disp(r);disp(p);
[b1,a1]=residuez(r,p,k);%反过来，由求出的r,p,k求多项式向量b,a
[r,p,k]=residuez(a,b);%交换分子分母多项式，再求相应r,p,k;

6、其它有关极零点的文件

clear;
B=[1,-3.3,7.25,-6.7,3,-0.8];
L=length(B);
A=zeros(1,L);
A(1)=1;
[Z,P,K]=tf2zp(B,A);
%tf2zp函数是找到分子为B，分母为A的多项式对应的极点和零点，极点序列为P,零点序列为Z，K为增益
sort(Z);
%sort函数将序列进行排序
[b,a]=zp2tf(Z,P,K);
%zp2tf函数是求出零点为Z，极点为P，增益为K所代表的多项式，b为分子，a为分母，与上面的B,A相对应
Z1=roots(B);
%roots函数可以求出由B组成的多项式的根
poly(Z1);
%poly函数将分解的根再综合成多项式

7、系统转移函数以及极零点与二阶子系统之间的转换函数

clear;
B=[0.0201,0,-0.0402,0,0.0202];
A=[1,-1.637,2.237,-1.307,0.641];
[sos,G]=tf2sos(B,A);%将转移函数分解为二阶子系统的级联；
[B,A]=sos2tf(sos,G);%由级联的二阶子系统重构转移函数
[z,p,k]=tf2zp(B,A);%求转移函数的极零点
[sos,G]=zp2sos(z,p,k);%由系统的极零点得到级联的二阶子系统
[z,p,k]=sos2zp(sos,G);%由级联的二阶子系统得到系统的极零点

第三章离散时间信号的傅里叶变换

1、fftfilt.m文件

例、令x(n)是正弦信号加白噪声信号，长度为500，h(n)是设计出来的低通FIR滤波器，长度为11，试用fftfilt实现长序列的卷积。

clear;
h=fir1(10,0.3,hanning(11));%设计低通滤波器，得到h(n)，长度为11，Wn是截止频率
N=500;
p=0.05;
f=1/16;
u=randn(1,N)*sqrt(p);%产生白噪声信号
s=sin(2*pi*f*[0:N-1]);
x=u(1,N)+s;
y=fftfilt(h,x);%实现叠加相加法滤波
subplot(211);plot(x);%画出原信号
subplot(212);plot(y);%画出滤波后的信号

2、hilbet.m文件

该文件用来计算x(n)的希尔伯特变换，调用格式为y=hilbert(x)，y的实部就是x(n),虚部就是x(n)的希尔伯特变换x~(n)；

clear all;
N=25;
f=1/16;
x=sin(2*pi*f*[0:N-1]);
y=hilbert(x);
% 求 x 的希尔伯特变换；
subplot(221)
stem(x,'.');
hold on;
plot(zeros(size(x)));xlabel('n');ylabel('x(n)'); 
subplot(222)
stem(imag(y),'.');
% imag(y) is the Hilbert transform of x(n);
hold on;
plot(zeros(size(x)));xlabel('n');ylabel('x~(n)');

第四章快速傅里叶变换

1、fft(x).m文件

该文件调用格式为X=fft(x)或者X=fft(x,N);第一种调用格式，如果x长度不是2的整次幂，实现的是慢速的非2的整数次幂的变换；对后者，N应该是2的整次幂，如果x的长度小于N，则补零，如果超过N，则舍弃N之后的数据，ifft是快速求反傅里叶变换的函数，与该函数调用方式相同。

例、令x(n)是两个正弦信号及其白噪声的叠加，对其做频谱分析

clear;
n=(1:256);
%产生两个正弦信号
x1=5*sin(2*pi*0.1*n);
x2=3*sin(2*pi*0.2*n);
N=256;
u=randn(1,N);%产生白噪声信号
x=x1+x2+u;
subplot(311);
plot(x(1:N/4));
xlabel('n');
ylabel('原信号');
%fft
y=fft(x);
f=-0.5:1/N:0.5-1/N;
subplot(312);
%fftshift的作用是看到（-fs,fs)的频率区间，方便查看是否产生混叠
plot(f,abs(fftshift(y)));
xlabel('f');
ylabel('幅度值');
%ifft
x1=ifft(y);
subplot(313);
plot(x1(1:N/4));
xlabel('n');
ylabel('恢复出的信号');

2、czt.m文件

其调用格式维X=czt(x,M,W,A),式中x是待变换的时域信号x(n),其长度设为N，M是变换的长度，W确定变换的步长，A确定变换的起点，若M=N，A=1,则CZT变为DFT。

例、求三个正弦信号相加的fft和czt

n=(1:128);
x1=sin(2*pi*8/40*n);
x2=sin(2*pi*8.22/40*n);
x3=sin(2*pi*9/40*n);
x=x1+x2+x3;
 
stepf=40/128;
fs=40;
n1=0:stepf:fs/2-stepf;%确定最终画出的频谱图范围为0~fs/2，便于观察
W=exp(-j*2*pi/128);%步长，2*pi/128即为相邻抽样点的角度差
y1=czt(x,128,W,1);
subplot(311);
plot(n1,abs(y1(1:128/2)));
xlabel('f');
ylabel('|czt(x)|');
 
y2=fft(x);%做fft
subplot(312);
plot(n1,abs(y2(1:128/2)));
xlabel('f');
ylabel('|fft(x)|');
 
% 详细构造A后的czt，由于查看了局部的频谱图，所以可以清楚的看到三个正弦信号的频域叠加效果
M=60;
f0=7.2;
DELf=0.05;
A=exp(j*2*pi*f0/fs);%起始点的Z变换
W=exp(-j*2*pi*DELf/fs);%步长更新
Y3=czt(x,M,W,A);
n2=f0:DELf:f0+(M-1)*DELf;
subplot(313);plot(n2,abs(Y3));grid on;

第五章离散时间系统相位、结构与状态变量描述

1、filtfilt.m文件

该文件实现零相位滤波，零相位滤波原理为：

文件调用格式为y=filtfilt(B,A,x);B是H（z）的分子多项式，A是H（z）的分母多项式，x是待滤波的信号，y是滤波后的信号。

clear;
b=[0.06745,0.1349,0.6745];
a=[1,-1.143,0.4128];
n=(-16:48);
%为什么是cos?
x=cos(0.1*pi*n)+cos(0.2*pi*n);
subplot(211);
stem(n,x,'.');
xlabel('n');
ylabel('x(n)');
y=filtfilt(b,a,x);
subplot(212);
stem(n,y,'.');
xlabel('n');
ylabel('y(n)');

2、grepdelay.m文件

该文件用于求一个系统的群延迟，调用格式为[gd w]=grpdelay(B,A,N)或者[gd F]=grpdelay(B,A,N,Fs);gd是求出来的群延迟，w以及F是频率分点，前者单位为rad，后者单位为Hz，二者长度均为N，FS为抽样频率，单位为HZ。

例、最小相位和线性相位系统的幅频响应

clear all;
N=20;
% 给定两个系统；
a1=[1 0 -0.81];
b1=[1 -1 0.75 -0.25 0.0625];
b2=[1.0000 4.0500 8.1000 14.9956 27.7248 43.2996 51.1831 43.2996...
27.7248 14.9956 8.1000 4.0500 1.00];
a2=1;
% 分别求这两个系统的群延迟；
[gd1,w]=grpdelay(b1,a1,256,1);
[gd2,w]=grpdelay(b2,a2,256,1);
subplot(411);plot(w,gd1);grid on;
xlabel('f/Hz');
ylabel('GD for H1');
subplot(412);plot(w,gd2);grid on;
xlabel('f/Hz');
ylabel('GD for H2');
%分别求这两个系统的相频响应
[H1,W1]=freqz(b1,a1,256,'whole',1);
Hphase1=angle(H1);
Hphase1=unwrap(Hphase1);
subplot(413);
plot(W1,Hphase1);ylabel('min 相频响应');xlabel('w/2pi');
[H2,w2]=freqz(b2,a2,256,'whole',1);
Hphase2=angle(H2);
Hphase2=unwrap(Hphase2);
subplot(414);
plot(w2,Hphase2);ylabel('线性 相频响应');xlabel('w/2pi');

由此结果可以看出，线性相位系统的群延迟是一个常数，最小相位系统群延迟先减后增。

3、tf2latc.m和latc2tf.m文件

这两个文件可以实现转移函数与Lattice系数之间的转换，调用格式为

K=tf2latc(b); k=tf2latc(1,a); [k,c]=tf2latc(b,a);

第一个式子对应全零系统，第二个对应全极系统，第三个对应极零系统，latc2tf文件与tf2latc文件调用格式刚好相反。以下程序说明两个文件的用法。

clear;
b=[1 -1.7 1.53 -0.648];%
k=tf2latc(b)
% all-zeros system to LATTICE; 
b1=latc2tf(k)
% LATTICE to all-zeros system;
a=b; 
k=tf2latc(1,a)
% all-poles system to LATTICE;
a1=latc2tf(k)
% LATTICE to all-poles system;

b=[1 0.8 -1 -0.8]; 
[k,c]=tf2latc(b,a)
% zeros-poles system to LATTICE;
[b2,a2]=latc2tf(k,c)
% LATTICE to zeros-poles system;

4、latcfilt.m文件

该文件用来实现Lattice结构下的信号滤波，调用格式为

[y,g]=latcfilt(k,x); 对应全零系统

[y,g]=latcfilt(k,l,x);对应全极系统

[y,g]=latcfilt(k,c,x);对应极零系统

x是带滤波器的信号，y是Lattice结构正向滤波的输出，g是做反向滤波的输出。

clear;
% 给定 IIR滤波器；
B=[0.0201 0 -0.0402 0 0.0201];
A=[1 -1.637 2.237 -1.307 0.641];
% 产生信号 x；
w1=0.1*pi;w2=0.35*pi;
N=100;n=0:N-1;
x=cos(w1*n)+cos(w2*n);
%直接滤波；
y1=filter(B,A,x);
% 求出lattice 系数；
[k,c]=tf2latc(B,A);
% 用lattice 系数滤波；
[y2,g]=latcfilt(k,c,x);
subplot(411);plot(x(10:N-1));xlabel('n');ylabel('x(n)');grid on;
subplot(412);plot(y1(10:N-1));xlabel('n');ylabel('直接滤波');grid on;
subplot(413);plot(y2(10:N-1));xlabel('n');ylabel('用Lattice系数滤波');grid on;
subplot(414);plot(g(10:N-1));xlabel('n');ylabel('反向滤波');grid on;

第六章无限冲激响应数字滤波器的设计

1、设计一个巴特沃斯低通滤波器,。

用到的文件有：

1、buttored.m:该文件是用来确定数字低通或者模拟低通滤波器的阶次，调用格式为数字滤波器：[N,Wn]=buttord(Wp,Ws,Rp,Rs);

模拟滤波器：[N,Wn]=buttord(Wp,Ws,Rp,Rs,’s’);

2、Buttap.m：该文件用来设计模拟低通原型滤波器G(p)，调用格式为：

[z,p,k]=buttap(N);z,p,k分别为设计出来的G(p)的极点、零点和增益。

3、lp2lp.m,lp2hp.m,lp2bp.m,lp2bs.m这四个文件功能为将模拟低通原型滤波器G(p)转换为实际的低通、高通、带通和带阻滤波器。调用格式分别为：

[B,A]=lp2lp(ba,a,W0);[B,A]=lp2hjp(b,a,W0);

[B,A]=lp2hp(b,a,W0,BW);[B,A]=lp2bs(b,a,W0,BW);

4、bilinear.m文件，该文件实现双线性变换，由模拟滤波器H(s)得到数字滤波器H(z),调用格式为[Bz,Az]=bilinear(B,A,Fs);

clear;
fp=100;fs=300;ap=3;as=20;Fs=1000;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
%将数字滤波器的指标转换为模拟滤波器指标,oup即为通带截止频率
oup=tan(wp/2);
ous=tan(ws/2);
%因为之前已经将数字滤波器指标转化为模拟滤波器，所以该Fs暂时可以不用
Fs=Fs/Fs;
%得到滤波器的阶数
[N,Wn]=buttord(oup,ous,ap,as,"s");
%由阶数求得G(p)的参数
[z,p,k]=buttap(N);
%由已知的极零点得到系统函数的分子分母系数
[b,a]=zp2tf(z,p,k);
%再由G(p)得到相应的系统函数H(s)
[B,A]=lp2lp(b,a,oup);

%将H(s)转换为H(z)
[Bz,Az]=bilinear(B,A,Fs/2);
[h,w]=freqz(Bz,Az,256,Fs*1000);
h=abs(h);
h=20*log10(h);
plot(w,h);grid on;xlabel('w/2pi');ylabel('幅频响应');

2、直接设计数字滤波器

butter.m文件可以直接用来设计数字滤波器，使用方法为

[B,A]=butter(N,Wn);[B,A]=butter(N,Wn,’high’);

[B,A]=butter(N,Wn,’stop’);[B,A]=butter(N,Wn,’s’);

Wn是通带截止频率，范围在0-1之间，1对应抽样频率的一半，B,A分别是H（z）分子分母多项式系数，最后一个用来设计模拟滤波器。

一步步计算的方法：
clear all;

fp=[300 400];fs=[200 500];
rp=3;rs=18;
Fs=2000;
wp=fp*2*pi/Fs;
ws=fs*2*pi/Fs;
wap=2*Fs*tan(wp./2)
was=2*Fs*tan(ws./2);
%得到带通模拟滤波器的幅频响应
[n,wn]=buttord(wap,was,rp,rs,'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)
%得到模拟角频率下的频率响应
w2=[0:Fs/2-1]*2*pi;
h2=freqs(bs,as,w2);
% Note: z=(2/Ts)(z-1)/(z+1);
%将模拟滤波器转化为数字滤波器
[bz1,az1]=bilinear(bs,as,Fs)
[h3,w3]=freqz(bz1,az1,1000,Fs);
plot(w2/2/pi,20*log10(abs(h2)),w3,20*log10(abs(h3)));grid;
ylabel('Bandpass AF and DF')
xlabel(' Hz')

直接使用butter函数计算的方法：
clear;
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')

3、设计切比雪夫Ⅰ型滤波器

用到的文件有：cheblord.m：计算滤波器的阶次，与buttord文件用法一样

Cheblap.m文件：用来设计原型低通切比雪夫Ⅰ型模拟滤波器。

Cheby1.m文件：用来直接设计数字切比雪夫Ⅰ型滤波器

间接法设计：
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];
wap=2*Fs*tan(wp./2);
was=2*Fs*tan(ws./2);
%先算借此n
[n,wn]=cheb1ord(wap,was,rp,rs,'s');
%得到G(p)的零极点以及增益
[z,p,k]=cheb1ap(n,rp);
%将零极点形式转化为分子分母形式
[bp,ap]=zp2tf(z,p,k);
%得到带宽
bw=wap(2)-wap(1);
%得到通带中心频率
w0=sqrt(wap(1)*wap(2));
%将G(p)转化为G(s)
[bs,as]=lp2bp(bp,ap,w0,bw);
%进行双线性变换，将模拟滤波器转化为数字滤波器
[bz1,az1]=bilinear(bs,as,Fs)
%得到数字滤波器的频率响应
[h,w]=freqz(bz1,az1,256,Fs);
plot(w,20*log10(abs(h)));grid on;

直接法设计：
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];
% 设计切比雪夫滤波器；
[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);xlabel('f');ylabel('幅频响应');grid on;

第七章有限冲激响应数字滤波器设计

1、几种常见窗函数及其频率响应

clear;
N=128;
n=(-15:15);

%矩形窗以及其归一化之后的对数幅频响应
b1=boxcar(31);
subplot(521);
stem(n,b1,'.');
xlim([-20 20]);
[h1,w]=freqz(b1,1,256,'whole',1);
h1=20*log10(abs(h1)/abs(h1(1)));
subplot(522);
w=w-0.5;
plot(w,fftshift(h1));
xlabel('w/2pi');ylabel('矩形窗归一化频率');

%三角窗以及其归一化之后的对数幅频响应
b2=bartlett(31);
subplot(523);
stem(n,b2,'.');
xlim([-20 20]);
[h2,w]=freqz(b2,1,256,'whole',1);
h2=20*log10(abs(h2)/abs(h2(1)));
subplot(524);
w=w-0.5;
plot(w,fftshift(h2));
xlabel('w/2pi');ylabel('三角窗归一化频率');

%汉宁窗以及其归一化之后的对数幅频响应
b3=hanning(31);
subplot(525);
stem(n,b3,'.');
xlim([-20 20]);
[h3,w]=freqz(b3,1,256,'whole',1);
h3=20*log10(abs(h3)/abs(h3(1)));
subplot(526);
w=w-0.5;
plot(w,fftshift(h3));
xlabel('w/2pi');ylabel('汉宁窗归一化频率');

%汉明窗以及其归一化之后的对数幅频响应
b4=hamming(31);
subplot(527);
stem(n,b4,'.');
xlim([-20 20]);
[h4,w]=freqz(b4,1,256,'whole',1);
h4=20*log10(abs(h4)/abs(h4(1)));
subplot(528);
w=w-0.5;
plot(w,fftshift(h4));
xlabel('w/2pi');ylabel('汉明窗归一化频率');

%布莱克曼窗以及其归一化之后的对数幅频响应
b4=blackman(31);
subplot(529);
stem(n,b4,'.');
xlim([-20 20]);
[h4,w]=freqz(b4,1,256,'whole',1);
h4=20*log10(abs(h4)/abs(h4(1)));
subplot(5,2,10);
w=w-0.5;
plot(w,fftshift(h4));
xlabel('w/2pi');ylabel('布莱克曼窗归一化频率');