Matlab幅频曲线和滤波器设计

---

前言少叙，下面开始正题。

一.离散数字信号的表示

n=-3:5;subplot(221);x1=(n==0);stem(n,x1,'.');title('单位冲击');axis([-4,4,-1,2]);grid on;subplot(222);x2=[n>=0];stem(n,x2,'.');title('单位阶跃');axis([-4,4,-1,2]);grid;subplot(223);x3=(n>=0&n<=2);stem(n,x3,'.');title('矩形序列');axis([-4,4,-1,2]);grid;%%%%%   %   用&&会产生错误提示：|| 和 && 运算符的操作数必须能够转换为逻辑标量值%   这说明A&&B，A和B不能是矩阵，只能是标量，而我们的n为1*9 double型矩阵，故只能用&subplot(224);x4=exp(1./2.*n);stem(n,x4,'.');title('指数序列');axis([-4,4,-1,2]);grid;

二.时间卷积和频率卷积

1.线性卷积

nx=[0:45];x1=[nx>=0];x2=[(nx-10)>=0];x=x1-x2;        %矩形方波信号nh=[0:45];h=(0.9).^nh;y=conv(x,h);ny=[0:60];subplot(311);stem(nx,x);grid;axis([-1 46 -0.2 1.2]);title('x');subplot(312);stem(nh,h);grid;axis([-1 46 -0.2 1.2]);title('h');subplot(313);stem(ny,y(1:61));grid;              %       卷积之后的长度是两个序列的长度之和减一。故y是一个1*91的矩阵，y(1:61)是提取前61个元素，即提取矩阵从第1列到第61列axis([-1 61 -2 8]);title('y=x*h');

要特别注意： C=conv（A,B）;计算两个有限长序列的卷积，卷积结果的长度为N+M-1。conv函数默认A和B表示的两个序列都是从0开始，且结果序列也是从0开始。当两个序列不是从0开始时，必须对conv函数加以扩展。
已知{[x(n);nx=xs:xf]},{[h(n);nh=hs:hf]},则y(n)=x(n)*h(n),且y(n)的起点和终点分别是ys=xs+hs,yf=xf+hf,下面写出通用卷积函数convu();

function [y,ny]=conv(h,nh,x,nx);ys=nx(1)+nh(1);yf=nx(end)+nh(end);y=conv(h,x);ny=ys:yf;end;%%%%%%下面对这个函数进行注释：convu()是通用卷积函数%   h和x都是有限长序列，y为卷积结果序列向量%   nh和nx分别是h和x的位置向量，ny是y的位置向量%   前两个end表示最后一个元素的下标

2.循环卷积

当循环卷积区间长度L大于等于N+M-1（即线性卷积的长度）时，循环卷积就等于线性卷积。

以下摘自：严素清等人的数字信号处理实训 (论文)说明书《用matlab实现两信号的卷积 》，详见百度文库.

function y=myconv(x1,x2,N)x1=input('x1=');x2=input('x2=');N=input('N=');if length(x1)>N    error('N必须大于或等于x1的长度');endif length(x2)>N    error('N必须大于或等于x2的长度'); endx1=[x1,zeros(1,N-length(x1))];x2=[x2,zeros(1,N-length(x2))];x3=circel(x2);y=x1*x3;        %   x1:1*N,x3:N*N.stem('y');xlabel('n');ylabel('y');grid on;title('循环卷积结果y')

当然，循环函数circel();
流程图如下：

function v=circlel(y)  %把一个序列转换成循环矩阵N=length(y);v=zeros(N,N);        %建立一个全为0的n行n列的矩阵if(rem(N,2)==1)     %这个if-else计算出原序列需要以N为周期            count=N/2+1;  %扩展变换后第一次转化时的次数else    count=N/2;   %偶数次的第一次变换次数endfor j=2:(count+1)     %把原序列变换成初始循环序列    mk=y(N-j+2);    y(N-j+2)=y(j);    y(j)=mk;endfor i=1:N            %把循环的起始序列转换成矩阵算法    k=1;    for j=1:N        v(i,j)=y(j);    end    L=y(N);    for k=N:-1:2           y(k)=y(k-1);    end    y(1)=L;end

有

3.时域卷积定理和频域卷积定理

三.绘制幅相特性曲线

频率响应函数：H(e^jw),即系统的单位脉冲响应h(n)的傅里叶变换
系统函数：H(z)，………..（同上）……………的z变换
两者有关系的条件是H(z)的收敛域包括单位圆|z|=1,也即是单位圆上的z变换是傅里叶变换。

1.求H(z)零极点，画出H(e^jw)幅频、相频特性曲线

b=[1 0.3 -0.1];a=[1 -0.9 0.25 -0.06];p=roots(a);z=roots(b);      %   roots();求的是方程组的根subplot(221);zplane(z,p);[h,w]=freqz(b,a,100);magh=abs(h);phah=angle(h);      %   求频响的幅值和相位subplot(222);plot(w/pi,magh);grid;title('magnitude');subplot(223);plot(w/pi,phah);grid;title('pha');

代码解释：

zplane(B,A);其中B.A分别是H(z)的分子和分母的系数向量，表示绘制出系统函数的零极点图.
[h,w]=freqz(B,A,M);计算出M个频率点上的频率响应，存放在h向量中.且freqz()自动将这M个频率点均匀设置在频率范围[0,pi]上.
注：freqz是z变换和FT用的，而freqs是拉氏变换用的。

2.函数传输函数求输出信号

系统函数同上，求单位冲击和单位阶跃经过系统输出的信号.

b=[1 0.3 -0.1];a=[1 -0.9 0.25 -0.06];n=[-10:30];x=[n==0];h1=filter(b,a,x);subplot(121);stem(n,h1);axis([-5,30,-0.5 1.5]);grid on;title('responce of impluse');x=[0<=n];h2=filter(b,a,x);subplot(122);stem(n,h2);axis([-5,30,-0.5 5]);grid on;title('responce of stemp');%%  filter(b,a,x);已知输入信号和传输函数，求输出信号

四.（数字）滤波器设计

1.巴特沃斯滤波器

a. 试给出巴特沃斯模拟滤波器的2.5.10.20.50阶的归一化幅频图。

clear;for q=1:5;    N=[2,5,10,20,50];   %   画5条条件一样的线，当然for语句比较好使    u=['b','r','g','m','k'];    N1=N(q);u1=u(q);    [z,p,k]=buttap(N1);    [b,a]=zp2tf(z,p,k);    n=0:0.02:2;    [h,w]=freqs(b,a,n);    H=(abs(h).^2);       %   归一化是幅度平方函数|Ha(jw)^2|=1/(1+(w/wc)^(2*N))    plot(w,H,u1);    hold on;endaxis([0,2,-0.2,1.2]);title('butterwoth');legend('N=2','N=5','N=10','N=20','N=50');grid on;hold off;

、

解释：
[z,p,k]=buttap(N);计算N阶巴特沃斯归一化（3DB截止频率Wc=1）的模拟低筒滤波器系统函数的零极点和增益因子。返回长度为N的列向量z和p,k表示滤波器增益。
若转化成一般式，可调用[b,a]=zp2tf(z,p,k);表示将系统函数的零极点转化为系统函数一般形式的系数，b和a分别是分子和分母多项式系数向量。
如图，上面为零极点形式，下面为一般式。

b.
已知通带截止频率850Hz，阻带截止频率2000Hz，波纹3DB，衰减25DB，采样频率10kHz，求数字低通巴特沃斯滤波器阶数，归一化转折频率和幅频相频特性曲线。

clc;wp=850;ws=2000;wc=5000;       wp=wp/wc;ws=ws/wc;      [N,wn]=buttord(wp,ws,3,25)[b,a]=butter(N,wn)freqz(b,a,1024,10000);axis([0,5000,-100,10]);

解释：
[N,wc]=buttord(wp,ws,Rp,As);用于计算巴特沃斯数字滤波器的阶数N和3db截止频率wc,wp.ws.wc均为归一化值，这就要求0<=wp<=1,0<=ws<=1,其中1表示数字频率pi，对应于模拟频率Fs/2,Fs位采样频率。Rp,As为通带最大衰减和阻带最小衰减。
[b,a]=butter(N,wc)；计算分子和分母多项式（一般式）的系数向量b和a。

注：若求时间序列信号的响应则

%%% x定义;y=filte(b,a,x);         %同上文

2.切比雪夫滤波器（略）

3.IIR滤波器

4.FIR滤波器