数字信号处理实验指导2student-201011new.doc

1、1实验项目表序号实验项目名称 内容提要实验要求实验类型实验时数每组人数实验室名称备注1利用 MATLAB语言对信号进行处理分析熟悉 MATLAB常用函数，并实现对数字信号的谱分析熟练 本科生 12信号、系统及系统响应连续信号经理想采样前后的频谱变化关系；时域离散系统时域特性；利用卷积方法观察分析系统的时域特性。熟练 本科生 13用 FFT作谱分析用 FFT对连续信号和时城离散信号进行谱分析掌握 本科生 14设计 IIR数字滤波器1 设计模拟低通滤波器。2 进行滤波器变换，得到数字低通滤波器。3 用 MATLAB语言实现。掌握 本科生 15设计 FIR滤波器利用窗函数设计法实现FIR滤波器设计，

2、并用MATLAB语言实现。掌握 本科生 167电子信息实验室目 录2实验一 离散时间的信号和系统 .2实验二 离散时间傅立叶变换 .4实验三 离散傅立叶变换（DFT） 7实验四 数字滤波器结构 10实验五 IIR 数字滤波器的设计 13实验六 FIR 滤波器的设计 窗函数法和频率抽样设计法 18实验七 语音处理系统仿真（综合滤波器设计） 22附录 自编函数 233实验一 离散时间的信号和系统一、实验目的1、复习离散时间的信号和系统，复习离散时间重要类型的信号和它们的运算的实现。2、复习离散时间信号理论中一些重要的结果，它们在数字信号处理中很有用。二、实验原理1、典型序列单位采样序列；单位阶跃序

3、列；实数指数序列；复数指数序列；正余弦序列；随机序列：MATLAB 可用 rand(1,N)和 randn(1,N)来生成；周期序列。2、序列的运算信号加；信号乘；改变比例 ；移位；折叠：fliplr(x)；取样和：sum(x(n1:n2)取样积：prod(x(n1:n2)；信号能量：sum(abs(x)2)；信号功率：sum(abs(x)2)/length(x)3、一些有用的结果单位采样合成：奇偶合成：几何级数：序列相关：kknxn)()()()()(oe1|,10aan对yx lnylr称 为 移 位),()(,4卷积运算：差分方程：在 Matlab 中：三、实验内容1.典型序列的实现单位

4、阶跃序列；实数指数序列；复数指数序列；正余弦序列；随机序列用 rand(1,N)和 randn(1,N)来生成；2.序列的运算给定序列 x1=1 2 3 4 5 6 7 8 9, ns1=-4; x2=9 8 7 6 5 4 3 2 1, ns2=4 求：1) x1+x2; 2) y3=x1x2；3) y1=0.5x1+0.8x2；4) y2=0.3x1(n)(n-6)+0.8(n-5)x2(n);5) x1 和 x2 的反折序列;6) x1(n)和 x2(n)的功率；7) y3=x1*x2 (线性卷积)；四、本实验用到的一些 MATLAB 函数Stem(x,y),Plot(x,y): x 轴

5、和 y 轴均为线性刻度（Linear scale） Xlable (x), Ylable(y)Title(x)axis(xmin,xmax,ymin,ymax):调整图轴的范围 Subplot(x,y,z):同时画出数个小图形在同一视窗中。Real(x):求复数 x 的实数部分；),(yxconvMmNkknaby01)()( ),(xfilter5Imag(x):求复数 x 的虚数部分；Abs(x):求复数 x 的模；Angle(x):求复数 x 的相位；Conv(x,y):求 x 和 y 的卷积，注意下标是从 1 开始；Fliplr(x):信号反折；Fliter(b,a,x):差分方程的实

6、现。五、MATLAB 编程和调试技巧因为 MATLAB 语言是一种解释性语言，所以有时 MATLAB 程序的执行速度不是很理想。因此尽量避免使用循环，用向量化的运算来代替循环操作。 注意;的使用。如果有的同学对 MATLAB 的使用不熟悉，请在老师处考取语言学习 PPT。六、实验报告的要求实验目的；实验原理；实验内容：要求有程序，有图形（坐标要标识清楚） 。实验结论。实验二 离散时间傅立叶变换一、实验目的1.复习离散时间傅立叶正反变换2.复习 DTFT 的两个重要特性63.复习 DTFT 的其它特性4.离散 LTI 系统的频率响应5.采样及重构信号二、实验原理1、信号的离散时间傅立叶变换（DT

7、FT）2、DTFT 的两个重要特性周期性：离散时间傅立叶变换是 w 的周期函数，其周期为 2。对称性：对于实值的 X(n)，是共扼对称的。即实部为偶对称，虚部为奇对称。3、DTFT 的其他特性线性；时移、共扼、折叠、卷积、乘法、能量等；4、LTI 系统的频率响应5、模拟信号的采样与重构采样定理重构: 步骤如下（a）先把样本集转换成一个加权脉冲串列（b）再将此脉冲串列通过一个带宽为 F的低通滤波器进行滤波。deXnxnxeXjwnjnjnjw)(21)( )()()( 2121 nxbFanbxaFnjwnjweheH)()(na Tstcxtx)(si)()(7以上两个步骤可用插值公式来描述：

8、三、实验内容1.求信号的离散时间傅立叶变换并分析其周期性和对称性；给定正弦信号 x(t)=2*cos(2*pi*10*t),fs=100HZ,求其 DTFT。2.用以下两个有限长序列来验证 DTFT的线性、卷积和共轭特性；x1(n)=1 2 3 4 5 6 7 8 9 10 11 12;x2(n)=R10(n)3. 求 LTI系统的频率响应给定系统 H（Z）=B（Z）/A（Z），A=0.98777 -0.31183 0.0256B=0.98997 0.989 0.98997，求系 统的幅频响应和相频响应。 （要求使用 filter(B,A,(n)求解。4. 采样和频谱混叠给定信号 x(t)=1

9、00*exp(-100*t)*cos(2*pi*500*t),求该信号的频谱；当采样频率分别为 fs1=2000HZ,fs2=1000HZ;fs3=500HZ;fs4=200HZ，时输出序列的 DTFT。四、本实验用到的函数有限长序列的 DTFT：dtft(x,n,w)（该函数需自己编写）算法： njwnjwexeX1)()( )(2)1(*)( )()(21)2()2(1)1()(2)2(1 nxxeeeeXnkjwnjjwnkjwnjjnkjnjjkjwjj 8用 Matlab语言来实现：function y=dtft(x,n,w)y=x*(exp(-j).(n*w); 五、实验报告的要求

10、实验目的；实验原理；实验内容：要求有程序，有图形（坐标要标识清楚） 。实验结论。实验三 离散傅立叶变换（DFT）一、实验目的1.掌握离散傅立叶级数(DFS)2.掌握 DFT变换3.掌握 DFT性质4.掌握利用 DFT计算线性卷积5.掌握快速傅立叶变换（FFT）二、实验原理1、离散傅立叶级数DFS 10102)()()()(NnnkNNnkj WxexxDFSkX9IDFS2、离散傅立叶变换（DFT）X(n)为有限长序列：DFT：IDFT：3、DFT 的性质线性性循环折叠性共轭性时序列的对称性序列的圆周移位频域中的圆周移位时域循环卷积频域循环卷积帕赛瓦尔定理4、用 DFT计算线性卷积10102)

11、()()()( NknkNNknj WXeXIDFSnx 10102)()()()( NnnkNnkj WxexDFTkX10102)()()()( NknkNknjXeXInx 101）设 x1(n)为 N1点的序列，x2(n)为 N2点的序列，则线性卷积的长度为 N1 N21，而当其两者的圆周卷积的长度取为N1N21 时，圆周卷积与线性卷积相等。2）实际中，当 x1(n)的长度远远大于和 x2(n)的长度时，往往采用分块卷积的方法，即重叠保留法和重叠相加法。5、快速傅立叶变换（FFT）在 MATLAB中提供 fft函数来计算 x的 DFT。Xfft(x,N);三、实验内容1.离散傅立叶级数

12、给定有限长序列1 2 3 4，延拓为周期 N=6 的周期序列，并求其DFS。2.求以下序列的 N=16,32,64点的快速傅立叶变换x(n)=exp(-i*pi/8*n);x2(n)=cos(pi/8*n);x3(n)=sin(pi/8*n)3.利用 DFT计算线性卷积给定序列 x1(n)=0.9n,n=0:11;h(n)=R9(n) 求 x1(n)*h(n)；x1(n)与h(n)的 10 点 圆周卷积。四、本实验用到的一些函数1.离散傅立叶级数：Xk=dfs(xn,N) 算法：xn=idfs(Xk,N) 2.离散傅立叶变换：Xk=dft(xn,N)10)()(NnnkWxkX )1()0(*

13、11)()0( 2)1(1NxWNXNN 11xn=idft(Xk,N)3.序列的循环移位：cirshftt(x,n,N)4.序列的循环卷积：circonvt(x1,x2,N)5.分块卷积：ovrlpsav(x,h,N)6.快速傅立叶变换：fft(x,N)五、实验报告的要求实验目的；实验原理；实验内容：要求有程序，有图形（坐标要标识清楚） 。实验结论。实验四 数字滤波器结构一、实验目的1.掌握 IIR滤波器的三种结构（直接形式、级联形式、并联形式）及其它们的互相转换。2.掌握线性相位 FIR滤波器的四种结构（横截形、级联形、线性相位形、频率抽样形）及其它们互相转换。二、实验原理1、IIR 滤波

14、器IIR滤波器可以写成：其差分方程可表示为：NMzazbbzABH10)(Mmmnynxby01)()()(12直接形式：可利用 y=filter(B,A,x)直接实现。级联形式：并联形式：2、FIR 滤波器FIR滤波器可表示为：其差分方程为：直接形式：可用函数 filter来实现。级联形式：线性相位形式：线性相位的条件：可表示为：h(n)=h(M-1-n) 偶对称h(n)=-h(M-1-n) 奇对称频率取样形式：冲激响应 h(n)的 M点 DFT为 H(k)则有：H(z)=Zh(n)=ZIDFT(H(k)利用内插公式可得：KkkkzABbz12,1,0)(HKk NMkkkk zCzz1 0

15、2,1,0)( )1(10)(HMzbzbz)110nxnxny21 21010 )()()Ni iiiNnn zZhzH（ weHj)( 10 2,1,0101)()() Mk iiMkk zAbBCzWHzH13三、实验内容IIR滤波器给定 IIR 滤波器 ，求其典范型、)8.09.1)(5.0( 61482) 232ZZH级联型、并联型结构。FIR滤波器给定 FIR 滤波器 ，求其 级联型、)3531(*2.0)( 42ZZZH频率抽样型结构。四、本实验用到的一些函数求多项式的根：b=roots(a)构造指定根的多项式：a=poly(b)部分分式展开：r1,p1,k=residuez(b

16、,a)b,a=residuez(r1,p1,k)对复数进行重新排序：p1=cplxpair(p)变直接形式为级联形式：b0,B,A=dir2cas(b,a)滤波器的级联实现：y=casfiltr(b0,B,A,x)变级联形式为直接形式：b,a=cas2dir(b0,B,A)变直接形式为并联形式：C,B,A=dir2par(b,a)NMzazbbzABH10)(KkkkzABbz12,1,0)(14滤波器的并联实现：y=parfiltr(b0,B,A,x)变并联形式为直接形式：b,a=par2dir(C,B,A)变 h(n)值形式为频率取样形式：C,B,A=dir2fs(h)五、实验报告的要求实

17、验目的；实验原理；实验内容：要求有程序，有图形（坐标要标识清楚） 。实验结论。实验五 IIR 数字滤波器的设计一、模拟滤波器的设计设滤波器的通带截止频率为 wp，阻带截止频率为 ws，通带衰减为 Rp，阻带衰减为 As，通带波纹为 Rp。1、根据给定指标得出，利用函数 buttord,cheb1ord,cheb2ord可以计算出低通模拟原型滤波器的阶数和截止频率。n,wn=buttord(wp,ws,Rp,As,s)NMzazbbzABH10)(Kk Mkkkk zCzz1 02,1,0)(10 2,1,0)()Mk iizAbBCzH15n,wn=cheb1ord(wp,ws,Rp,As,s

18、)2、求出低通模拟原型滤波器后，再利用函数z,p,k=Buttap(n)；z,p,k=Cheb1ap(n,Rp) 求出零极点型系统函数，然后可用函数zp2tf-bap,aap=zp2tf(z,p,k)将其转化为 b/a型3、将模拟低通原型滤波器经频率变换为所要求的模拟滤波器（低通，高通，带通，带阻）可用函数 lp2lp，lp2hp，lp2bs， lp2bp。b,a=lp2lp(bap,aap,wo)b,a=lp2bs(bap,aap,wo)4、求出滤波器的幅频，相频及冲激响应。db,mag,pha,w=freqs_m(b,a,Wmax)(见附录）ha,x,t=impulse(b,a)(matl

19、ab自带）直接画出幅频特性和相频特性：freqs(b,a) (matlab 自带）例：wp=0.2*pi;Rp=1;ws=0.3*pi;As=15;n,wn=buttord(wp,ws,Rp,As,s)z,p,k=buttap(n);bap,aap=zp2tf(z,p,k)b,a=lp2lp(bap,aap,wn);db,mag,pha,w=freqs_m(b,a,0.5*pi);ha,x,t=impulse(b,a);figure(1)16subplot(2,2,1),plot(w/pi,mag);title(幅度); subplot(2,2,2),plot(w/pi,db);title(幅

20、度 in db);subplot(2,2,3),plot(w/pi,pha/pi);title(相位); subplot(2,2,4),plot(t,ha),title(冲激响应);figure(2)freqs(b,a);二、数字滤波器的设计设滤波器的通带截止频率为 Wp，止带截止频率为 Ws，通带衰减为 Rp，止带衰减为 As，通带波纹为 Rp。设计方法一：1、根据给定指标得出，低通模拟原型滤波器的阶数和截止频率。可利用函数 buttord, cheb1ord,cheb2ord完成。n,wn=buttord(Wp,Ws,Rp,As)n,wn=cheb1ord(Wp,Ws,Rp,As)2、直接

21、求出数字滤波器，利用函数 butter,cheby1,cheby2。b,a=butter(n,Wn， ftype)b,a=cheby1(n, Rp,Wn,ftype)ftype指的是：low,bandpass,high,stop3、求出滤波器的幅频，相频及群延时。db,mag,pha,grd,w=freqz_m(b,a);(见附录)直接画出幅频特性和相频特性：freqz(b,a) (matlab 自带）例：17wp=0.2;Rp=1;ws=0.3;As=15;n,wn=buttord(wp,ws,Rp,As)b,a=butter(n,wn);db,mag,pha,grd,w=freqz_m(b

22、,a);figure(3)subplot(2,2,1),plot(w/pi,mag);title(幅度);subplot(2,2,2),plot(w/pi,db);title(幅度 in db); subplot(2,2,3),plot(w/pi,pha/pi);title(相位);subplot(2,2,4),plot(w/pi,grd),title(群延时);figure(4);freqz(b,a)设计方法二：1冲激响应不变法使用冲激响应法来实现由模拟滤波器滤波器转变为数字滤波器。利用所编写的函数实现：b,a=imp_invr(c,d,T)b,a=impinvar(cs,ds,T)( Ma

23、tlab自带)2、双线性变换法利用双线性法来实现模拟滤波器到数字滤波器的转换。 12zTsTskkezAs118b,a=bilinear(cs,ds,T)(Matlab自带）三、设计要求1、设计低通模拟滤波器，具体要求如下：3db带宽 10KHZ，阻带起始频率 12KHZ，阻带应达到的最小衰减-30db 方法一：采用 Matlab自带函数来实现。方法二：采用自编函数实现。c,d=afd_butt(wp,ws,Rp,As)c,d=afd_cheb1(wp,ws,Rp,As)2、设计数字滤波器（1）带通滤波器：通带范围：100HZ-200HZ ； 阻带范围：截止频率 60 HZ；240HZ ，通带

24、允许最大衰减 3dB, 阻带应达到的最小衰减-30dB（2）低通滤波器 wp=pi/3, ws=pi/2通带允许最大衰减 3dB, 阻带应达到的最小衰减-30dB。（3）阻带范围：100HZ-200HZ 通带截止频率 60HZ和240HZ。通带允许最大衰减 3dB, 阻带应达到的最小衰减-30dB四、实验报告的要求实验目的；实验原理；实验内容：要求有程序，有图形（坐标要标识清楚） 。实验结论。1920实验六 FIR 滤波器的设计窗函数法和频率抽样设计法求四种线性相位 FIR滤波器的振幅响应函数（见附录）hr,w,a,L=hr_type1(h);h偶对称，N 为奇数，h(n)=h(N-1-n)h

25、r,w,a,L=hr_type2(h); h偶对称，N 为偶数，h(n)=h(N-1-n)hr,w,a,L=hr_type3(h); h奇对称，N 为奇数,h(n)=-h(N-1-n)hr,w,a,L=hr_type4(h); h奇对称，N 为偶数,h(n)=-h(N-1-n)一、窗函数设计法原理窗函数法是从时域出发，把理想的 hd(n)用一定形状的窗函数截取成有限长的 h(n)，以此来近似理想的 hd(n),这样得到的频率响应 H(ejw)逼进于所要求的理想频率响应。(1)先给定所要求设计的理想滤波器的频率响应 Hd(ejw).(2)由于 hd(n)是无限长的，而我们必须用有限长的 h(n)

26、来逼近hd(n)，最有效的方法就是截断 hd(n)。h(n)=w(n)* hd(n)dweHnhnheHjnjwddnjjw)(21)(21二、窗函数设计法设计方法一1、根据实际阻带衰减指标，来确定所使用的窗函数。Matlab 提供了几个函数来实现这些窗函数。W=boxcar(N); 矩形窗 21dBW=triang(N); 三角窗 25dBW=hanning(N); 汉宁窗 44dBW=hamming(N); 海明窗 53dBW=blackman(N); 布莱克曼窗 74dBW=kaiser(N,beta); 凯泽窗(beta=0.7865) 80dB2、根据过渡带来计算出 N值。例： 通带

27、截止频率 wp, 阻带截止频率 ws，已知为汉宁窗：N=3.1*2*pi/(wp-ws);）（ 凯 泽 窗布 拉 克 曼 窗海 明 窗汉 宁 窗三 角 形 窗矩 形 窗 阻 带 最 小 衰 减 （过 渡 带 宽旁 瓣 峰 值 幅 度 （窗 函 数 865.7 80/2*574/. 53/2*3.41 4/1.3 25/2*.25 1/9.013 ) NNdBwdB223、求出理想的 hd(n)。可使用 ideal_lp 来实现理想低通滤波器的冲激响应。例：由 ideal_lp来实现理想带阻滤波器的冲激响应。hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(w

28、c2,N);3、求得所设计的 FIR滤波器的单位抽样响应:h=hd.*w例 1ws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;wd=min(wp1-ws1),(ws2-wp2);wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;w_bla=(blackman(N+1);N=ceil(11*pi/wd);hd=ideal_lp(wc2,N+1)-ideal_lp(wc1,N+1);h=hd.*w_bla;设计方法二1、根据实际阻带衰减指标，来确定所使用的窗函数。2、根据过渡带来计算出 N值。3、利用 Matlab所提供的函数 fir1，来

29、实现 fir滤波器。h=fir1(N,wn,ftype,windows(N+1),对于高通滤波器和带阻滤波器，N 必须为偶数， ftype指的、是：low,bandpass,high,stop例 2设计理想带通滤波器：23ws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;wd=min(wp1-ws1),(ws2-wp2);wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;N=ceil(11*pi/wd);wn=wc1,wc2/pi;b=fir1(N,wn,blackman(N+1);freqz(b,1)三、频率抽样设计法频率抽样法从频域出

30、发把给定的理想频率响应 Hd(ejw)加以等 间隔抽样，然后以此 Hd(k)作为实际 FIR滤波器的频率特性的抽样值 H(k)，知道 H(k)后，由 DFT定义，可用频域的这 N个抽样值 H(k)求唯一确定的有限长序列 h(n)，利用这 N个频域抽样值 H(k)同样可得FIR滤波器的系统函数 H（Z）及频率响应 。1、根据给定的 N值和过渡点来求得 Hd(k) 。例 3N=40;T1=0.5925;T2=0.1099;ws1=0.2*pi;wp1=0.35*pi;ws2=0.8*pi;wp2=0.65*pi;alpha=(N-1)/2;L=0:N-1;wl=(2*pi/N)*L;hrs=zer

31、os(1,5),T2,T1,ones(1,7),T1,T2,zeros(1,9),T2,T1,ones(1,7),T1,T2,zeros(1,4);k1=0:floor(N-1)/2);24k2=floor(N-1)/2)+1:N-1;angh=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2);H=hrs.*exp(j*angh);2、利用 DFT反变换求的所需的 h(n)。h=real(ifft(H,N);四、设计任务完成例 1-例 3程序的调试，并分析原理。 五、实验报告的要求实验目的；实验原理；实验内容：要求有程序，有图形（坐标要标识清楚） 。实验结论。实

32、验七 语音处理系统仿真（综合滤波器设计） 一、实验目的：通过本实验,使学生加深对数字信号处理系统结构的理解，通过不同采样频率对输出频谱和重构信号的影响，加深学生对采样定理的认识。通过不同类型滤波器的设计得到不同的输出频谱变化和语音处理效果，使学生进一步了解 FIR 和 IIR 滤波器的特性，通过本实验， 使学生掌握基本的时域、频域信号处理方法，同时提高学生综合应用所学理论知识分析问题解决问题的能力。二、实验要求：251） 掌握采样定理、FFT、滤波器设计原理及数字信号处理系统的结构。2） 熟悉常用数字信号处理中的 MATLAB 函数。3） 设计系统框图，对各个模块的输出输入的时域和频域图形进行

33、记录。4） 进行数据分析，完成实验报告。三、实验内容：（1） 信号预处理与采样：通过 MATLAB 函数读取一段语音信号，通过不同模拟滤波器的设计完成抗混叠滤波，通过设置不同的采样频率实现采样。（2） 均衡处理（综合滤波器设计）通过设计 FIR、IIR 低通、带通、高通滤波器实现信号在低、中、高频段的均衡。四、实验方案（任务）：任务 1：用声音编辑工具软件录制 Wav文件，获取处理语音信号。 任务 2：抗混叠滤波器设计（IIR， fc=10KHZ） 。任务 3：信号的采样（fs=5KHZ,10KHZ,20KHZ,40KHZ） ，并观测频谱混叠现象。任务 4：设计均衡滤波器（低通-5KHZ，高通

34、-5KHZ，带通-1KHZ8KHZ）五、实验仪器设备 ：MATLAB 常用函数附录 自编函数26function b,a = afd_butt(Wp,Ws,Rp,As);% Analog Lowpass Filter Design: Butterworth% -% b,a = afd_butt(Wp,Ws,Rp,As);% b = Numerator coefficients of Ha(s)% a = Denominator coefficients of Ha(s)% Wp = Passband edge frequency in rad/sec; Wp 0% Ws = Stopband

35、edge frequency in rad/sec; Ws Wp 0% Rp = Passband ripple in +dB; (Rp 0)% As = Stopband attenuation in +dB; (As 0)%if Wp 0% Ws = Stopband edge frequency in rad/sec; Ws Wp 0% Rp = Passband ripple in +dB; (Rp 0)% As = Stopband attenuation in +dB; (As 0)%if Wp 0% Ws = Stopband edge frequency in rad/sec;

36、 Ws Wp 0% Rp = Passband ripple in +dB; (Rp 0)% As = Stopband attenuation in +dB; (As 0)%if Wp 0% Ws = Stopband edge frequency in rad/sec; Ws Wp 0% Rp = Passband ripple in +dB; (Rp 0)% As = Stopband attenuation in +dB; (As 0)%if Wp = 0error(Passband edge must be larger than 0)endif Ws = Wperror(Stopb

37、and edge must be larger than Passband edge)endif (Rp = 0) | (As 0)error(PB ripple and/or SB attenuation ust be larger than 0)endep = sqrt(10(Rp/10)-1);A = 10(As/20);OmegaC = Wp;k = Wp/Ws;k1 = ep/sqrt(A*A-1);capk = ellipke(k.2 1-k.2); % Version 4.0 codecapk1 = ellipke(k1 .2) 1-(k1 .2); % Version 4.0

38、codeN = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1);fprintf(n* Elliptic Filter Order = %2.0f n,N)b,a=u_elipap(N,Rp,As,OmegaC);function w_black = Blackman(M);% M-point Blackman window% -% w_black = Blackman(M);%M1 = M-1;m = 0:1:M1;w_black = abs(0.42 - 0.5*cos(2*pi*m/(M1) + 0.08*cos(4*pi*m/(M1);function b

39、,a = cas2dir(b0,B,A);% CASCADE-to-DIRECT form conversion% -29% b,a = cas2dir(b0,B,A)% b = numerator polynomial coefficients of DIRECT form% a = denominator polynomial coefficients of DIRECT form% b0 = gain coefficient% B = K by 3 matrix of real coefficients containing bks% A = K by 3 matrix of real

40、coefficients containing aks%K,L = size(B);b = 1;a = 1;for i=1:1:Kb=conv(b,B(i,:);a=conv(a,A(i,:);endb = b*b0;function y = casfiltr(b0,B,A,x);% CASCADE form realization of IIR and FIR filters% -% y = casfiltr(b0,B,A,x);% y = output sequence% b0 = gain coefficient of CASCADE form% B = K by 3 matrix of

41、 real coefficients containing bks% A = K by 3 matrix of real coefficients containing aks% x = input sequence%K,L = size(B);N = length(x);w = zeros(K+1,N);w(1,:) = x;for i = 1:1:Kw(i+1,:) = filter(B(i,:),A(i,:),w(i,:);endy = b0*w(K+1,:);function b,a = cheb1hpf(wp,ws,Rp,As)% IIR Highpass filter design

42、 using Chebyshev-1 prototype% b,a = cheb1hpf(wp,ws,Rp,As)% b = Numerator polynomial of the highpass filter% a = Denominator polynomial of the highpass filter% wp = Passband frequency in radians% ws = Stopband frequency in radians% Rp = Passband ripple in dB30% As = Stopband attenuation in dB% Determ

43、ine the digital lowpass cutoff frequecies:wplp = 0.2*pi;alpha = -(cos(wplp+wp)/2)/(cos(wplp-wp)/2);wslp = angle(-(exp(-j*ws)+alpha)/(1+alpha*exp(-j*ws);% Compute Analog lowpass Prototype Specifications:T = 1; Fs = 1/T;OmegaP = (2/T)*tan(wplp/2);OmegaS = (2/T)*tan(wslp/2);% Design Analog Chebyshev Pr

44、ototype Lowpass Filter:cs,ds = afd_chb1(OmegaP,OmegaS,Rp,As);% Perform Bilinear transformation to obtain digital lowpassblp,alp = bilinear(cs,ds,Fs);% Transform digital lowpass into highpass filterNz = -alpha,1; Dz = 1,alpha;b,a = zmapping(blp,alp,Nz,Dz);function xec, xoc = circevod(x)% signal decom

45、position into circular-even and circular-odd parts% -% xec, xoc = circecod(x)%if any(imag(x) = 0)error(x is not a real sequence)endN = length(x); n = 0:(N-1);xec = 0.5*(x + x(mod(-n,N)+1);xoc = 0.5*(x - x(mod(-n,N)+1);function y = circonvt(x1,x2,N)% N-point circular convolution between x1 and x2: (time-domain)% -% y = circonvt(x1,x2,N)% y = output sequence containing the circular convolution% x1 = input sequence of length N1 = N% x2 = input sequence of length N2 = N% N = size of circular buffer% Method: y(n) = sum (x1(m)*x2(n-m) mod N)