数字信号处理与MATLAB实现.doc

数字信号处理与MATLAB实现摘要：采用FFT算法计算序列向量X的N点DFT，缺省N时，fft函数自动按X的长度计算DFT。 . 24. cheby1.m. 直接设计切比雪夫1型滤波器。 以上3个文件的调用格式和对应的巴特沃 .关键词：算法,点,24类别：专题技术来源：牛档搜索（Niudown.COM）本文系牛档搜索（Niudown.COM）根据用户的指令自动搜索的结果，文中内涉及到的资料均来自互联网，用于学习交流经验，作品其著作权归原作者所有。不代表牛档搜索（Niudown.COM）赞成本文的内容或立场，牛档搜索（Niudown.COM）不对其付相应的法律责任！数字信号处理与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.mx,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 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.mx,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.my,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):nx1nnx2,h(n):nh1nnh2, 卷积结果序列为y(n):nx1+nh1nnx2+nh2例. 设,求程序：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因此计算，得Z反变换b=1;a=poly(0.9 0.9 -0.9);r,p,k=residuez(b,a)结果：r = 0.2500 0.5000 0.2500p = 0.9000 0.9000 -0.9000k =因此得到 相应的5. 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. 傅里叶变换信号时域和频域变换连续对应着非周期，离散对应着周期。一个域的离散必然导致另一个域的周期延拓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). k=0, 1, ，N-112. 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),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)=N= abs(A)L=pow2(nextpow2(M+N-1) M, N分别为序列延拓周期14. 用DFT进行谱分析时，必须将序列阶段为长度为N的有限长序列。造成频谱泄露和谱间干扰。泄露使频谱变得模糊，分辨率降低；旁瓣引起不同分量间的干扰。截断效应无法完全消除，可以加宽窗和缓慢截断。矩形窗比海明窗的频率分辨率高（泄漏小），但谱间干扰大，因此海明窗是以牺牲分辨率来换取谱间干扰的降低。栅栏效应是只能在离散点的地方看到真实的像，其余频谱被遮挡。为减少栅栏效应，可以在时域数据末端增加一些零点，使周期内点数增加，但不改变原有数据，即增加频域抽样点数N，频域抽样为，这样必然使谱线更密，这样原来看不到的谱分量就可能看到了。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如： %直接型b=1，-3，11，27，18；a=16，12，2，-4，-1 %级联型b0=4;B=1,1,0;1,-1.4142136,1;A=1.-0.5.0;.81 %并联型C=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(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分别是通带和阻带的截止频率，实际上它们是归一化频率，其值在01之间，1对应抽样频率的一半。对低通和高通滤波器，Wp,Ws都是标量，对带通和带阻滤波器，Wp,Ws都是12的向量。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是通带截止频率，范围在01之间，1对应抽样频率的一半。若Wn是标量，则格式(1)用来设计低通数字滤波器，若Wn是12的向量，则格式(1)用来设计数字带通滤波器；格式（2）用来设计数字高通滤波器；格式（3）用来设计数字带阻滤波器，显然，这时的Wn是12的向量：格式（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.m29. 用脉冲响应不变法设计低通数字滤波器, 模拟滤波器采用切比雪夫1型滤波器原型模拟滤波器采用切比雪夫1型滤波器原型程序shuziqieshi.m脉冲响应不变法子程序imp_invr数字滤波器响应子程序freqz_m直接型转换成并联型子程序dir2par.m比较两个含同样标量元素,但(可能)有不同下标的复数对及其相应留数向量子程序cplxcomp.m用脉冲响应不变法设计数字滤波器，模拟滤波器采用巴特沃兹滤波器原型shuzibashi.m30. 模拟低通滤波器变换为数字高通滤波器Alp2Dhp.m模拟低通滤波器变换为数字带通滤波器Alp2Dbp.m模拟低通滤波器变换为数字带阻滤波器Alp2Dbs.m31. 数字低通到数字滤波器的设计为了从低通数字滤波原型得到新数字滤波器的有理函数，必须实现有理代换，通常采用zmapping.m函数实现。32. 用zmapping函数实现高通滤波器实例zmappinghp.m33. 一个含有高频噪声的信号保存在数据文件noisbloc中，试用脉冲响应不变法设计低通滤波器对其进行除噪。lpchuzao.m34. 无限长单位冲击响应（IIR）数字滤波器的优点是可以利用模拟滤波器设计的结果，而模拟滤波器的设计有大量图表可查，方便简单。但是它也有明显的缺点，就是相位的非线性，在图像处理以及数据传输等要求信道具有线性相位特性的场合，IIR滤波器就不太适用了。有限单位冲击响应（FIR）数字滤波器则可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。此外，FIR滤波器的单位抽样相应是有限长的，因而FIR滤波器一定是稳定的。再有，只要经过一定的延时，任何非因果有限长序列都能变成因果的有限长序列，因而总能用因果系统来实现。最后，FIR滤波器由于单位冲击响应是有限长的，因而可以用快速傅里叶变换（FFT）算法来过滤信号，从而可大大提高运算效率。但是，要取得很好的衰减特性，FIR滤波器H(z)的阶次比IIR滤波器的要高。IIR滤波器设计中的各种变换法对FIR滤波器设计是不适用的，这是因为那里是利用有理分式的系统函数，而FIR滤波器的系统函数只是的多项式。因为最感兴趣的是具有线性相位的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是12的向量，则用来设计带通滤波器；若Wn是1L的向量，则可用来设计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相对应的所希望的幅频响应。不指定窗函数的类型时，将自动选择汉明窗。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的长度为奇数，即N应为偶数。38. remexor