1、数字信号处理与 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 单位抽样序列 生成函数 impseq.m)(0nx,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
2、0 0n=-5:5;x1=3*impseq(2,-4,5)-impseq(-4,-5,4) %起点到终点长度要一致x1 =0 -1 0 0 0 0 3 0 0 04.2 单位阶跃序列 生成函数 stepseq.m)(nux,n=stepseq(no,ns,nf) %序列的起点为 ns，终点为 nf，在 n=n0 点处生成一个单位阶跃4.3 两个信号相加的生成函数 sigadd.my,n=sigadd(x1,n1,x2,n2)4.4 两个信号相乘的生成函数 sigmult.my,n=sigmult(x1,n1,x2,n2)4.5 序列移位 y(n)=x(n-n0)的生成函数 sigshift.m
3、y,n=sigshift(x,m,n0)4.6 序列翻褶 y(n)=x(-n)的生成函数 sigfold.my,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)(zzX122534)(zzX )
4、()(21zXY程序：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 23 19 15n =-3 -2 -1 0 1 2因此 2123 5937)( zzzY计算 ， 得 Z 反变换121)9.0().() zzzX9.0b=1;a=poly(0.9 0.9 -0.9);r,p,k=residuez(b,a)结果：r =0.25000.50000.2500p =0.90000.9000-0.9000k =因此得到 ,9.0125)9.01(5.20)( 12zzzzX 9.0相应的 )(.)()5
5、nunununn5. H,w=freqzB,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. 傅里叶变换信号时域和频域变换连续对应着非周期，离散对应着周期
6、。一个域的离散必然导致另一个域的周期延拓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 是向量，则
7、此傅里叶变换即变成一维傅里叶变换 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
8、和 idct：离散余弦正变换和离散余弦逆变换Y=dct(x,N)完成如下变换，N 的默认值为 length(x).k=0, 1, ，N-1nnkxKY1)12(cos)(2)(12. fftshiftY=fftshift(x)用来重新排列 X=fft(x)的输出，当 X 为向量时，把 X 的左右两半进行交换，从而将零频分量移至频谱的中心；如果 X 为二维傅里叶变换的结果，它同时将 X 的左右和上下部分进行交换。13. fftfiltY=fftfilt(b,x)采用重叠相加法 FFT 对信号向量 x 快速滤波，得到输出序列向量 y，向量 b 为 FIR 滤波的单位脉冲相应，h(n)=b(n+1)
9、,n=0,1,length(b)-1.Y=fftfilt(b,x,N)自动选取 FFT 长度 NF=2nextpow2(N), 输入数据 x 分段长度 M=NF-length(b)+1, 其中nextpow2(N)函数求的一个整数，满足 2(nextpow2(N)-1)= abs(A)L=pow2(nextpow2(M+N-1) M, N 分别为序列延拓周期14. 用 DFT 进行谱分析时，必须将序列阶段为长度为 N 的有限长序列。造成频谱泄露和谱间干扰。泄露使频谱变得模糊，分辨率降低；旁瓣引起不同分量间的干扰。截断效应无法完全消除，可以加宽窗和缓慢截断。矩形窗比海明窗的频率分辨率高（泄漏小）
10、 ，但谱间干扰大，因此海明窗是以牺牲分辨率来换取谱间干扰的降低。栅栏效应是只能在离散点的地方看到真实的像，其余频谱被遮挡。为减少栅栏效应，可以在时域数据末端增加一些零点，使周期内点数增加，但不改变原有数据，即增加频域抽样点数 N，频域抽样为 ，这样必然使谱线更密，这样原来看不到的谱分量就可能看到了。k215. 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如： %直接型5432161873)( zzzHb=1，-3，11，27，
11、18；a=16，12，2，-4 ，-1%级联型)81.09.1)(5.0(42364) 2zzzHb0=4;B=1,1,0;1,-1.4142136,1;A=1.-0.5.0;1.0.9.0.81%并联型21213)( zzzC=0;B=-14,-12;24,26;A=1,-2,3;1,-1,1delta=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(b
12、0,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,R
13、s,s)格式(1)对应数字滤波器, 式中 Wp,Ws 分别是通带和阻带的截止频率，实际上它们是归一化频率，其值在 01 之间，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,
14、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)的分子
15、、分母多项式的系数向量；在格式（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等文件都包含进去，
16、从而使设计过程更简捷。其调用格式是(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是通带截止频率，范围在 01 之间，1 对应抽样频率的一半。若 Wn 是标量，则格式(1)用来设计低通数字滤波器，若 Wn 是 1 2 的向量，则格式(1)用来设计数字带通滤波器；格式（2）用来设计数字高通滤波器；格式（3）用来设计数字带阻滤波器，显然，这时的 Wn是 1 2 的向量：
17、格式（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模拟巴特沃兹滤波器低通滤波器设计实例 mo
18、nibashi.m其中，巴特沃兹模拟滤波器的设计子程序设计和非归一化巴特沃兹模拟低通滤波器原型子程序 afd_butt.m计算系统函数的幅度响应和相位响应子程序 freqs_m.m直接形式转换成级联形式子程序 sdir2cas.m模拟切比雪夫 1 型低通滤波器设计实例 moniqieshi.m（2 型.m 文件中改为 2）其中，切比雪夫 1 型模拟滤波器的设计子程序设计和非归一化切比雪夫 1 型模拟低通滤波器原型子程序 afd_chb1.m计算系统函数的幅度响应和相位响应子程序 freqs_m.m直接形式转换成级联形式子程序 sdir2cas.m29. 用脉冲响应不变法设计低通数字滤波器, 模
19、拟滤波器采用切比雪夫 1 型滤波器原型模拟滤波器采用切比雪夫 1 型滤波器原型程序 shuziqieshi.m脉冲响应不变法子程序 imp_invr数字滤波器响应子程序 freqz_m直接型转换成并联型子程序 dir2par.m比较两个含同样标量元素,但(可能) 有不同下标的复数对及其相应留数向量子程序cplxcomp.m用脉冲响应不变法设计数字滤波器，模拟滤波器采用巴特沃兹滤波器原型 shuzibashi.m30. 模拟低通滤波器变换为数字高通滤波器 Alp2Dhp.m模拟低通滤波器变换为数字带通滤波器 Alp2Dbp.m模拟低通滤波器变换为数字带阻滤波器 Alp2Dbs.m31. 数字低通
20、到数字滤波器的设计为了从低通数字滤波原型得到新数字滤波器的有理函数，必须实现有理代换，通常采用zmapping.m 函数实现。32. 用 zmapping 函数实现高通滤波器实例 zmappinghp.m33. 一个含有高频噪声的信号保存在数据文件 noisbloc 中，试用脉冲响应不变法设计低通滤波器对其进行除噪。lpchuzao.m34. 无限长单位冲击响应（IIR）数字滤波器的优点是可以利用模拟滤波器设计的结果，而模拟滤波器的设计有大量图表可查，方便简单。但是它也有明显的缺点，就是相位的非线性，在图像处理以及数据传输等要求信道具有线性相位特性的场合，IIR 滤波器就不太适用了。有限单位冲
21、击响应（FIR）数字滤波器则可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。此外，FIR 滤波器的单位抽样相应是有限长的，因而 FIR 滤波器一定是稳定的。再有，只要经过一定的延时，任何非因果有限长序列都能变成因果的有限长序列，因而总能用因果系统来实现。最后，FIR 滤波器由于单位冲击响应是有限长的，因而可以用快速傅里叶变换（FFT）算法来过滤信号，从而可大大提高运算效率。但是，要取得很好的衰减特性，FIR 滤波器 H(z)的阶次比 IIR 滤波器的要高。IIR 滤波器设计中的各种变换法对 FIR 滤波器设计是不适用的，这是因为那里是利用有理分式的系统函数，而 FIR 滤波器的系统函
22、数只是 的多项式。1z因为最感兴趣的是具有线性相位的 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 是通带截止频率，其值在 01 之间，1 对应抽样频率的一半；b 是设计好的滤波器系数 h(n).对于第一种格式，若 Wn 是一标量，则可用来设计低通滤波器；若 Wn 是 1 2 的向量，
23、则用来设计带通滤波器；若 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 是频率向量，其值在 01 之间，M 是与 F 相对应的所希望的幅频响应。不指定窗函数的类型时，将自动选择汉
24、明窗。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 是频率向量，其值在 01 之间；A 是对应 F 的各频段上的理想幅频响应；W 是各频段上的加权向量。值得注意的是，若 b 的长度为偶数，涉及高通和带阻滤波器时有可能出现错误，因此最好保证 b 的长度为奇数，
25、即 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 结合使用
26、。值得说明是，remeord 给出的阶次 N 有可能偏低，这是适当增加 N 即可；另外，若 N 为奇数，就可令其加 1，使其变为偶数，这样 b 的长度为奇数。39. sgolay.m本文件用来设计 Savitzky-Golay 平滑滤波器。其调用格式是 b=sgolay（k,f）式中，k 是多项式的阶次，f 是拟合的双边点数。要求 kf，且 f 为奇数。40. firls.m本文件用最小平方法设计线性相位 FIR 数字滤波器。可设计任意给定的理想幅频特性。41. fircls.m本文件用带约束的最小平方法设计线性相位 FIR 数字滤波器。可设计任意给定的理想幅频特性。42. fircls1.m本文件用带约束的最小平方法设计线性相位 FIR 低通和高通滤波器。可设计任意给定的理想幅频特性。43. firrcos.m本文件用来设计低通线性相位 FIR 数字滤波器，其过渡带为余弦函数形式。44. 低通滤波器用矩形窗和汉明窗观察其频谱响应 lpchuang.m45. 设计多带滤波器 duodaifilter.m46. 利用切比雪夫最佳一致逼近法设计一低通滤波器 qiebijinlp.m47. 利用切比雪夫最佳一致逼近法设计一多阻带陷波器 qiebijinduobs
