MATLAB常用函数

MATLAB常用函数
数字信号处理与MATLAB 实现
1． n1=[ns:nf]；
x1=[zeros(1,n0-ns),1,zeros (1,nf-n0)]； %单位抽样序列的产生
2. subplot(2,2,4) 画2行2列的第4个图
3. stem(n,x) %输出离散序列，(plot 连续)
4. 编写子程序可调用
4.1 单位抽样序列)(0n n -δ生成函数impseq.m
[x,m]=impseq(n0,ns,nf); %序列的起点为ns ，终点为nf ，在n=n0点处生成一个单位脉冲
n=[-5:5];x1=3*impseq(2,-5,5)-impseq(-4,-5,5)
x1 =
0 -1 0 0 0 0 0 3 0 0 0
n=[-5:5];x1=3*impseq(2,-4,5)-impseq(-4,-5,4) %起点到终点长度要一致
x1 =
0 -1 0 0 0 0 3 0 0 0
4.2 单位阶跃序列)(0n n u -生成函数stepseq.m
[x,n]=stepseq(no,ns,nf) %序列的起点为ns ，终点为nf ，在n=n0点处生成一个单位阶跃
4.3 两个信号相加的生成函数sigadd.m
[y,n]=sigadd(x1,n1,x2,n2)
4.4 两个信号相乘的生成函数sigmult.m
[y,n]=sigmult(x1,n1,x2,n2)
4.5 序列移位y(n)=x(n-n0)的生成函数sigshift.m
[y,n]=sigshift(x,m,n0)
4.6 序列翻褶y(n)=x(-n)的生成函数sigfold.m
[y,n]=sigfold(x,n)
4.7 evenodd.m 函数可以将任一给定的序列x(n)分解为xe(n)和
xo(n)两部分
[xe,xo,m]=evenodd(x,n)
4.8 序列从负值开始的卷积conv_m, conv 默认从0开始
function [y,ny]=conv_m(x,nx,h,nh)
有{x(n):nx1≤n ≤nx2},{h(n):nh1≤n ≤nh2}, 卷积结果序列为
{y(n):nx1+nh1≤n ≤nx2+nh2}
例. 设1132)(-++=z z z X ,1225342)(-+++=z z z z X ,求)()()(21z X z X z Y += 程序：
x1=[1,2,3];n1=-1:1;
x2=[2,4,3,5];n2=-2:1;
[y,n]=conv_m(x1,n1,x2,n2)
结果：
y =
2 8 17 2
3 19 15
n =
-3 -2 -1 0 1 2
因此21231519231782)(--+++++=z z
z z z z Y 计算1121))9.01()9.01(()(---+-=z z z X ，9.0>z 得Z 反变换
b=1;a=poly([0.9 0.9 -0.9]);
[r,p,k]=residuez(b,a)
结果：
r =
0.2500
0.5000
0.2500
p =
0.9000
0.9000
-0.9000
k =
[] 因此得到,9.0125.0)9.01(5.09.0125.0)(1
211---++-+-=z z z z X 9.0>z 相应的)()9.0(25.0)1()9.0)(1(9
5)()9.0(25.0)(1n u n u n n u n X n n n -++++=+ 5. [H,w]=freqz[B,A,M]
计算出M 个频率点上的频率响应，存放于H 向量中，M 个频率存放在w 向量中，freqz 函数自动将这M 个频率点均匀设置在频率范围[0，π]之间。

若缺省w 和M 时，函数自动选取512个频率点计算。

不带输出向量的freqz 函数将自动绘制幅频和相频曲线。

也可[H,w]=freqz(B,A);
plot(w/pi,abs(H)); 绘出幅频特性
6. zplane(z,p)
绘制出列向量z 中的零点（以符号o 表示）和列向量p 中的极点（以符号x 表示）以及参考单位圆。

7. 傅里叶变换信号时域和频域变换
连续对应着非周期，离散对应着周期。

一个域的离散必然导致另一个域的周期延拓
8. fft 和ifft ：一维快速傅里叶变换和逆傅里叶变换
X=fft(x,N)
采用FFT 算法计算序列向量X 的N 点DFT ，缺省N 时，fft 函数自动按X 的长度计算DFT 。

当N 为2的整数次幂时，fft 按基2算法计算，否则用混合基算法。

ifft 的调用格式类似。

9. fft2和ifft2: 二维快速正傅里叶变换和逆傅里叶变换
（1）Y=fft2(x)
数据二维傅里叶变换参数X 是向量fft2(x)相当于fft(fft(x)’)’,即先对X 的列做一维傅里叶变换，然后再对变换结果的行作一维傅里叶变换；若X 是向量，则此傅里叶变换即变成一维傅里叶变换fft; 若X 是矩阵，则是计算该矩阵的二位傅里叶变换。

（2）Y=fft2(X, M, N)
通过对X 进行补零或截断，使得X 成为（M*N ）的矩阵。

函数iff2的参数应用与函数fft2完全相同。

10. czt: 线形调频Z 变换
Y=czt(x,m,w,a)
此函数计算由z=a*w.^(-(0:m-1))定义的z 平面螺旋线上各点的z 变换，a 规定了起点，w 规定了相邻点的比例，m 规定了变换的长度，后三个变量默认值为a=1,w=exp(j*2*pi/m)及m=length(x)因此y=czt(x)就等于y=fft(x).
11. dct 和idct ：离散余弦正变换和离散余弦逆变换
Y=dct(x,N)
完成如下变换，N 的默认值为length(x).
∑=+=N n n k n n x K Y 1))12(2cos(
)(2)(π k=0, 1, …，N-1
12. fftshift
Y=fftshift(x)
用来重新排列X=fft(x)的输出，当X 为向量时，把X 的左右两半进行交换，从而将零频分量移至频谱的中心；如果X 为二维傅里叶变换的结果，它同时将X 的左右和上下部分进行交换。

13. fftfilt
Y=fftfilt(b,x)
采用重叠相加法FFT 对信号向量x 快速滤波，得到输出序列向量y ，向量 b 为FIR 滤波的单位脉冲相应，h(n)=b(n+1),n=0,1,…,length(b)-1.
Y=fftfilt(b,x,N)
自动选取FFT 长度NF=2^nextpow2(N), 输入数据x 分段长度M=NF-length(b)+1, 其中nextpow2(N)函数求的一个整数，满足2^(nextpow2(N)-1)<=N<=2^nextpow2(N)
N 缺省时，fftfilt 自动选择合适的FFT 长度NF 和对x 的分段长度M 。

p = nextpow2(A)， p 满足 2^p >= abs(A)
L=pow2(nextpow2(M+N-1)) M, N 分别为序列延拓周期
14. 用DFT 进行谱分析时，必须将序列阶段为长度为N 的有限长序列。

造成频谱泄露和谱间干扰。

泄露使频谱变得模糊，分辨率降低；旁瓣引起不同分量间的干扰。

截断效应无法完全消除，可以加宽窗和缓慢截断。

矩形窗比海明窗的频率分辨率高（泄漏小），但谱间干扰大，因此海明窗是以牺牲分辨率来换取谱间干扰的降低。

栅栏效应是只能在离散点的地方看到真实的像，其余频谱被遮挡。

为减少栅栏效应，可以在时域数据末端增加一些零点，使周期内点数增加，但不改变原有数据，即增加频域抽样点数N ，频域抽样为k N π2，这样必然使谱线更密，这样原来看不到的谱分量就可能看到了。

15. x=[one(1,5),zero(1,N-5)]; %单位阶跃信号
y = filter(b,a,x) %直接型输出信号
16. N=25;
h=impz(b,a,N) %N 次采样，b,a 为零极点多项式系数，a 不算1 如：5
43214
321421216118271131)(-----------++++++-=z z z z z z z z z z H %直接型 b=[1，-3，11，27，18]；a=[16，12，2，-4，-1] )
81.09.01)(5.01()4142136.11)(1(4)(211211------++-+-+=z z z z z z z H %级联型b0=4;B=[1,1,0;1,-1.4142136,1];A=[1.-0.5.0;1.0.9.0.81]
2
11
211126243211214)(------+--++---=z z z z z z z H %并联型C=0;B=[-14,-12;24,26];A=[1,-2,3;1,-1,1]
delta=impseq(0,0,N);
x=[one(1,5),zero(1,N-5)];
h=filter(b,a,delta); %直接型单位脉冲响应
y = filter(b,a,x); %直接型输出信号
h=casfilter(b0,B,A,delta); %级联型单位脉冲响应
h=casfilter(b0,B,A,x); %级联型输出响应
h=parfiltr(C,B,A,delta); %并联型单位脉冲响应
h=parfiltr(C,B,A,x); %并联型输出响应
[b0,B,A]=dir2cas(b,a); %直接型转换成级联型
[b,a]=cas2dir(b0,B,A); %级联型转换成直接型
[C,B,A]=dir2par(b,a); %直接型转换成并联型
[b,a]=par2dir(C,B,A); %并联型转换成直接型
17. buttord.m
用来确定数字低通或模拟低通滤波器的阶次，其调用格式分别是
(1)[N,Wn]=buttord(Wp,Ws,Rp,Rs)
(2) [N,Wn]=buttord(Wp,Ws,Rp,Rs,’s ’)
格式(1)对应数字滤波器, 式中Wp,Ws 分别是通带和阻带的截止频率，实际上它们是归一化频率，其值在0~1之间，1对应抽样频率的一半。

对低通和高通滤波器，Wp,Ws 都是标量，对带通和带阻滤波器，Wp,Ws 都是1?2的向量。

Rp,Rs 分别是通带和阻带的衰减，单位为dB 。

N 是求出的相应低通滤波器的阶次，Wn 是求出的3dB 频率，它和Wp 稍有不同。

格式（2）对应模拟滤波器，式中各个变量的含义和格式（1）相同，但Wp,Ws 及Wn 的单位为rad/s, 因此，他们实际上是频率。

18. buttap.m
用来设计模拟低通原型滤波器G(p), 其调用格式是[z,p,k]=buttap(N)
N 是与设计的低通原型滤波器的阶次，z,p 和k 分别是设计出的G(p)的极点、零点及增益。

19. lp2lp.m, lp2hp.m, lp2bp.m, lp2bs.m
以上4个文件的功能是将模拟低通原型滤波器G(p)分别转换为实际的低通、高通、带通及带阻滤波器。

其调用格式分别为
（1）[B,A]=lp2lp(b,a,Wo) 或 [B,A]=lp2hp(b,a,Wo)
（2）[B,A]=lp2bp(b,a,Wo,Bw) 或 [B,A]=lp2bs(b,a,Wo,Bw)
式中b,a 分别是模拟低通原型滤波器G(p)的分子、分母多项式的系数向量，B,A 分别是转换后的H(s)的分子、分母多项式的系数向量；在格式（1）中，Wo 是低通或高通滤波器的截止
频率；在格式（2）中Wo是带通或带阻滤波器的中心频率，Bw 是其带宽。

20. bilinear.m
实现双线性变换，即由模拟滤波器H(s)得到数字滤波器H(z)。

其调用格式是
[Bz,Az]=bilinear(B,A,Fs)式中B,A分别是H(s)的分子、分母多项式的系数向量，Bz,Az
分别是H(z)的分子、分母多项式的系数向量，Fs是抽样频率。

21. butter.m
用来直接设计巴特沃兹数字滤波器，实际上它把buttord.m, buttap.m, lp2lp.m及bilinear.m等文件都包含进去，从而使设计过程更简捷。

其调用格式是
(1) [B,A]=butter(N,Wn)
(2) [B,A]=butter(N,Wn,’high’)
(3) [B,A]=butter(N,Wn,’stop’)
(4) [B,A]=butter(N,Wn,’s’)
格式(1)~(3)用来设计数字滤波器，B、A分别是H(z)的分子、分母多项式的系数向量，Wn 是通带截止频率，范围在0~1之间，1对应抽样频率的一半。

若Wn是标量，则格式(1)用来设计低通数字滤波器，若Wn是1?2的向量，则格式(1)用来设计数字带通滤波器；格式（2）用来设计数字高通滤波器；格式（3）用来设计数字带阻滤波器，显然，这时的Wn是1?2的向量：格式（4）用来设计模拟滤波器。

22. cheb1ord.m
求切比雪夫1型滤波器的阶次。

23. cheb1ap.m
用来设计原型切比雪夫1型模拟滤波器。

24. cheby1.m
直接设计切比雪夫1型滤波器。

以上3个文件的调用格式和对应的巴特沃兹滤波器的文件类似。

25. impinvar.m
用冲激响应不变法实现ω到Ω及s到z的转换。

26. maxflat.m
设计广义巴特沃兹低通滤波器。

27. zp2tf.m
把零极点加增益的形式，改写成零极点多项式 [b,a] = zp2tf(z,p,k)
28. 低通数字滤波器设计例子lpexample.m
模拟巴特沃兹滤波器低通滤波器设计实例monibashi.m
其中，巴特沃兹模拟滤波器的设计子程序设计和非归一化巴特沃兹模拟低通滤波器原型子程序afd_butt.m
计算系统函数的幅度响应和相位响应子程序freqs_m.m
直接形式转换成级联形式子程序sdir2cas.m
模拟切比雪夫1型低通滤波器设计实例moniqieshi.m（2型.m文件中改为2）
其中，切比雪夫1型模拟滤波器的设计子程序设计和非归一化切比雪夫1型模拟低通滤波器原型子程序afd_chb1.m
计算系统函数的幅度响应和相位响应子程序freqs_m.m
直接形式转换成级联形式子程序sdir2cas.m
29. 用脉冲响应不变法设计低通数字滤波器, 模拟滤波器采用切比雪夫1型滤波器原型
模拟滤波器采用切比雪夫1型滤波器原型程序shuziqieshi.m
脉冲响应不变法子程序imp_invr
数字滤波器响应子程序freqz_m
直接型转换成并联型子程序dir2par.m
比较两个含同样标量元素,但(可能)有不同下标的复数对及其相应留数向量子程序cplxcomp.m
用脉冲响应不变法设计数字滤波器，模拟滤波器采用巴特沃兹滤波器原型shuzibashi.m 30. 模拟低通滤波器变换为数字高通滤波器
Alp2Dhp.m
模拟低通滤波器变换为数字带通滤波器Alp2Dbp.m
模拟低通滤波器变换为数字带阻滤波器Alp2Dbs.m
31. 数字低通到数字滤波器的设计
为了从低通数字滤波原型得到新数字滤波器的有理函数，必须实现有理代换，通常采用zmapping.m函数实现。

32. 用zmapping函数实现高通滤波器实例zmappinghp.m
33. 一个含有高频噪声的信号保存在数据文件noisbloc中，试用脉冲响应不变法设计低通滤波器对其进行除噪。

lpchuzao.m
34. 无限长单位冲击响应（IIR）数字滤波器的优点是可以利用模拟滤波器设计的结果，而模拟滤波器的设计有大量图表可查，方便简单。

但是它也有明显的缺点，就是相位的非线性，在图像处理以及数据传输等要求信道具有线性相位特性的场合，IIR滤波器就不太适用了。

有限单位冲击响应（FIR）数字滤波器则可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。

此外，FIR滤波器的单位抽样相应是有限长的，因而FIR滤波器一定是稳定的。

再有，只要经过一定的延时，任何非因果有限长序列都能变成因果的有限长序列，因而总能用因果系统来实现。

最后，FIR滤波器由于单位冲击响应是有限长的，因而可以用快速傅里叶变换（FFT）算法来过滤信号，从而可大大提高运算效率。

但是，要取得很好的衰减特性，FIR滤波器H(z)的阶次比IIR 滤波器的要高。

IIR滤波器设计中的各种变换法对FIR滤波器设计是不适用的，这是因为那里是利用有理分式的系统函数，而FIR滤波器的系统函数只是1-z的多项式。

因为最感兴趣的是具有线性相位的FIR滤波器，对非线性相位的FIR滤波器，一般可以用IIR滤波器来代替。

所以只讨论线性相位滤波器。

FIR数字滤波器的文件
35. fir1.m
本文件采用窗函数法设计FIR数字滤波器，其调用格式是
b=fir1(N,Wn)
b=fir1(N,Wn,’high’)
b=fir1(N,Wn,’stop’)
式中N为滤波器的阶次，因此滤波器的长度为N+1;Wn是通带截止频率，其值在0~1之间，1对应抽样频率的一半；b是设计好的滤波器系数h(n).
对于第一种格式，若Wn是一标量，则可用来设计低通滤波器；若Wn是1?2的向量，则用来设计带通滤波器；若Wn是1?L的向量，则可用来设计L带滤波器，此时，格式将变为
b=fir1(N,Wn,’DC-1’)或b=fir1(N,Wn,’DC-0’)
其中，前者保证第一个带为通带，后者保证打一个带为阻带。

第二种格式用来设计高通滤波器，第三种格式用来设计带阻滤波器。

值得注意的是，在上述所有格式中，若不指定窗函数的类型，则fir1自动选择汉明窗。

36. fir2.m
本文件采用窗函数法设计具有任意幅频特性的FIR滤波器。

其调用格式是
b=fir1(N,F,M)
其中F是频率向量，其值在0~1之间，M是与F相对应的所希望的幅频响应。

不指定窗函数的类型时，将自动选择汉明窗。

37. remez.m
本文件用来设计采用切比雪夫最佳一致逼近FIR数字滤波器。

同时，还可以用来设计希尔伯特变换器和差分器。

其调用格式是b=remez(N,F,A)
b=remez(N,F,A,W)
b=remez(N,F,A,W, ’hilbert’)
b=remez(N,F,A,’differentiator’)
其中，N是给定的滤波器的阶次；b是设计的滤波器的系数，其长度为N+1;F是频率向量，其值在0~1之间；A是对应F的各频段上的理想幅频响应；W是各频段上的加权向量。

值得注意的是，若b的长度为偶数，涉及高通和带阻滤波器时有可能出现错误，因此最好保证b的长度为奇数，即N应为偶数。

38. remexord.m
本文件采用切比雪夫一致逼近设计FIR数字滤波器时所需要的滤波器阶次。

其调用格式是[N,Fo,Ao,W]=remexord(F,A,DEV,Fs) 式中，F、A的含义同文件remez.m，是通带和阻带上的偏差；该文件输出的是符合要求的滤波器阶次N、频率向量Fo、幅度向量Ao 和加权向量W。

若设计者事先不能确定自己要设计的滤波器阶次，那么调用remexord后，就可利用这一族参数再调用remez，即[N,Fo,Ao,W]=remexord(F,A,DEV,Fs)，从而设计出所需要的滤波器。

因此，通常remez和remexord结合使用。

值得说明是，remeord给出的阶次N有可能偏低，这是适当增加N即可；另外，若N为奇数，就可令其加1，使其变为偶数，这样b 的长度为奇数。

39. sgolay.m
本文件用来设计Savitzky-Golay平滑滤波器。

其调用格式是b=sgolay（k,f）
式中，k是多项式的阶次，f是拟合的双边点数。

要求k<f，且f为奇数。

< bdsfid="275" p=""></f，且f为奇数。

<>
40. firls.m
本文件用最小平方法设计线性相位FIR数字滤波器。

可设计任意给定的理想幅频特性。

41. fircls.m
本文件用带约束的最小平方法设计线性相位FIR数字滤波器。

可设计任意给定的理想幅频特性。

42. fircls1.m
本文件用带约束的最小平方法设计线性相位FIR低通和高通滤波器。

可设计任意给定的理想幅频特性。

43. firrcos.m
本文件用来设计低通线性相位FIR数字滤波器，其过渡带为余弦函数形式。

44. 低通滤波器用矩形窗和汉明窗观察其频谱响应lpchuang.m
45. 设计多带滤波器duodaifilter.m
46. 利用切比雪夫最佳一致逼近法设计一低通滤波器qiebijinlp.m
47. 利用切比雪夫最佳一致逼近法设计一多阻带陷波器qiebijinduobs
在MATLAB中，可以用函数y=filter(p,d,x)实现差分方程的仿真，也可以用函数y=conv(x,h)计算卷积。

（1）即y=filter(p,d,x)用来实现差分方程，d表示差分方程输出y 的系数，p表示输入x的系数，而x表示输入序列。

输出结果长度数等于x的长度。

实现差分方程，先从简单的说起：
filter([1,2],1,[1,2,3,4,5])，实现y[k]=x[k]+2*x[k-1]
y[1]=x[1]+2*0=1 (x[1]之前状态都用0)
y[2]=x[2]+2*x[1]=2+2*1=4
（2）y=conv(x,h)是用来实现卷级的，对x序列和h序列进行卷积，输出的结果个数等于x的长度与h的长度之和减去1。

卷积公式：z(n)=x(n)*y(n)= ∫x(m)y(n-m)dm.
程序一：以下两个程序的结果一样
（1）h = [3 2 1 -2 1 0 -4 0 3]; % impulse response
x = [1 -2 3 -4 3 2 1]; % input sequence
y = conv(h,x);
n = 0:14;
subplot(2,1,1);。