《2023年数字信号处理上机指导.docx》由会员分享，可在线阅读，更多相关《2023年数字信号处理上机指导.docx（23页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、湖南科技学院电子工程系 数字信号处理试验试验一 离散时间信号的MATLAB实现1. 正弦序列离散正弦序列的MATLAB表示与连续信号类似，只不过是用stem函数而不是用plot函数6sin p k12来画出序列的波形。下面就是正弦序列的MATLAB源程序。程序运行结果如图1.1所示。%正弦序列实现程序k=0:39;fk=sin(pi/6*k); stem(k,fk)2. 指数序列离散指数序列的一般形式为ca图1.1正弦序列波形k ，可用MATLAB中的数组幂运算即点幂运算c* a. k 来实现。下面为用MATLAB编写绘制离散时间实指数序列波形的函数。function dszsu(c,a,k1
2、,k2)%c：指数序列的幅度%a：指数序列的底数%k1：绘制序列的起始序号%k2：绘制序列的终止序号k=k1:k2;x=c*(a.k); stem(k,x,”filled”) hold on plot(k1,k2,0,0) hold off5353利用上述函数，实现实指数波形MATLAB程序如下其中a 值分别为 ， ，-444 ，- 4 。%离散时间实指数序列实现程序subplot 221; dszsu(1,5/4,0,20);xlabel(”k”);title(”f1k”); subplot 222 dszsu(1,3/4,0,20);xlabel(”k”);title(”f2k”); su
3、bplot 223; dszsu(1,-5/4,0,20);xlabel(”k”);title(”f3k”); subplot 224; dszsu(1,-3/4,0,20);xlabel(”k”);title(”f4k”);分析程序运行结果,对于离散时间实指数序列ca k ，当a 确实定值大于1时，序列为随时间发散的序列，当a 确实定值小于1时，序列为随时间收敛的序列。同时可见，当a 的值小于零时，其波形在增长或衰减的同时，还交替地转变序列值的符号。对于离散时间虚指数序列，可用通过调用以下绘制虚指数序列时域波形的MATLAB函数。function=dxzsu(n1,n2,w)%n1:绘制波形
4、的虚指数序列的起始时间序号%n2：绘制波形的虚指数序列的终止时间序号%w：虚指数序列的角频率k=n1:n2;f=exp(i*w*k); Xr=real(f) Xi=imag(f) Xa=abs(f) Xn=angle(f)subplot(2,2,1), stem(k,Xr,”filled”),title(”实部”);subplot(2,2,3), stem(k,Xi,”filled”),title(”虚部”);subplot(2,2,2), stem(k,Xa,”filled”),title(”模”);subplot(2,2,4), stem(k,Xn,”filled”),title(”相角”
5、);利用上述函数，实现虚指数波形MATLAB程序如下其中虚指数分别为e j%离散时间虚指数实现程序figure(1); dxzsu(0,20,pi/4); figure(2); dxzsu(0,20,2);kp4 ，e j 2k 程序运行结果如图1.21a、b所示。由图可见，只有当虚指数序列的角频率满足2pw为有理数时，信号的实部和虚部和相角都为周期序列，否则为非周期序列。a e jkp4 波形(b) e j 2k 波形图1.2 虚指数序列波形对于复指数序列，其一般形式为f k = r k e jwk可以通过调用下面绘制复指数序列时域波形的MATLAB函数。function dfzsu(n1,
6、n2,r,w)%n1:绘制波形的虚指数序列的起始时间序号%n2：绘制波形的虚指数序列的终止时间序号%w：虚指数序列的角频率%r: 指数序列的底数k=n1:n2; f=(r*exp(i*w).k; Xr=real(f); Xi=imag(f); Xa=abs(f); Xn=angle(f);subplot(2,2,1), stem(k,Xr,”filled”),title(”实部”);subplot(2,2,3), stem(k,Xi,”filled”),title(”虚部”);subplot(2,2,2), stem(k,Xa,”filled”),title(”模”);subplot(2,2,
7、4), stem(k,Xn,”filled”),title(”相角”);利用上述函数，实现复指数序列波形MATLAB程序如下。%复指数序列实现程序(r1) figure(1); dfzsu(0,20,1.2,pi/4);%复指数序列实现程序(0r1时，复指数序列的实部和虚局部别为幅度按指数增长的正弦序列;当0r1时，复指数序列的实部和虚局部别为幅度按指数衰减的正弦序列;当r=1时，复指数序列的实部和虚局部别为等幅正弦序列。3. 单位抽样序列可以通过借助MATLAB中的零矩阵函数zeros表示。全零矩阵zeros(1,N)产生一个由N个零组成的列向量，对于有限区间的d k 可以通过以下MATLA
8、B程序表示% 单位抽样序列实现程序k=-30:30;delta=zeros(1,30),1,zeros(1,30); stem(k,delta)4. 单位阶跃序列可以通过借助MATLAB中的单位矩阵函数ones表示。单位矩阵ones(1,N)产生一个由N个1 组成的列向量，对于有限区间的uk 可以通过以下MATLAB程序表示% 单位阶跃序列实现程序k=-30:30;uk=zeros(1,30),ones(1,31); stem(k,uk)程序运行结果如图1.24所示。例 1.7编写程序来产生以下根本脉冲序列。（1） 单位抽样序列，起点n =0，终点n =10，在n=3 处有一单位脉冲。sf0（
9、2） 单位阶跃序列，起点n =0，终点n =10，在n =3 前为 0，在n =3 后为 1。s3复指数序列， =-0.2，解：程序清单如下。clear,n0=3; ns=0;nf=10;f0s。=0.50n1=ns:nf; x1=zeros(1,n0-ns), 1, zeros(1,nf-n0);%单位抽样序列的产生n2=ns:nf; x2=zeros(1,n0-ns), ones(1,nf-n0+1);%单位阶跃序列的产生n3=ns:nf; x3=exp(-0.2+0.5j)*n3);%复指数序列的产生subplot(2,2,1),stem(n1,x1);title(”单位脉冲序列(n-3
10、)”); %画图表示单位抽样序列subplot(2,2,2),stem(n2,x2);title(”单位阶跃序列 u(n-3)”); %画图表示单位阶跃序列subplot(2,2,3),stem(n3, real(x3); line(0, 10,0, 0)title(”复指数序列”);ylabel(”实部”);%画图表示复指数序列的实部subplot(2,2,4),stem(n3, imag(x3); line(0, 10,0, 0)title(”复指数序列”);ylabel(”虚部”);试验二 利用 MATLAB 实现频谱分析21-1例 1有限长序列 x(n)的长度 N = 4 ，且 x(n
11、)= 3n = 0n = 1n = 2n = 3，用 FFT 求 X (k )，再用 IFFT 求 x(n)。解：利用快速傅里叶变换函数求解的MATLAB 实现程序清单如下clearxn=1,2,-1,3;X=fft(xn) x=ifft(X)程序运行结果如下X =5.00002.0000 + 1.0000i-5.00002.0000 - 1.0000ix =12-13例 2设 x(n)是由两个正弦信号及白噪声的叠加，试用FFT 文件对其作频谱分析。图 2.1 运行结果解：程序清单如下% 产生两个正弦加白噪声；N=256; f1=.1;f2=.2;fs=1; a1=5;a2=3;w=2*pi/
12、fs;x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+randn(1,N);% 应用FFT 求频谱；subplot(2,2,1);plot(x(1:N/4);title(”原始信号”);f=-0.5:1/N:0.5-1/N;X=fft(x); y=ifft(X); subplot(2,2,2);plot(f,fftshift(abs(X); title(”频域信号”);subplot(2,2,3);plot(real(x(1:N/4); title(”时域信号”);运行结果如图 4.14 所示，该程序同时完成了傅里叶变换与傅里叶反变换。例 3 设 x(n)为
13、长度 N = 6 的矩形序列，用 MATLAB程序分析FFT 取不同长度时 x(n)的频谱变化。解： N = 8,32,64 时 x(n)的 FFT MATLAB实现程序如下x=1,1,1,1,1,1;N=8;y1=fft(x,N);n=0:N-1;subplot(3,1,1);stem(n,abs(y1),”.k”);axis(0,9,0,6); N=32;y2=fft(x,N);n=0:N-1;subplot(3,1,2);stem(n,abs(y2),”.k”);axis(0,40,0,6); N=64;y3=fft(x,N);n=0:N-1;subplot(3,1,3);stem(n,
14、abs(y3),”.k”);axis(0,80,0,6);运行结果 x(n)的频谱如图 2.2 所示图 2.2 不同长度时 x(n)的频谱图例 3 x (t) = e-1000|t| ，求其傅立叶变换 Xa( jw)。a解：严格说，在MATLAB 中不使用symbolic 工具箱是不能分析模拟信号的，但是当采样时间间隔充分小的时候，可产生平滑的图形，当时间足够长，可显示出全部的模型，也就h是可以近似的分析。此例中， x (t)为 f =2023Hz 的带限信号，因此取1aT = 5 10- 5 2 f，因此2ssh采样信号的频谱发生混叠造成取样信号不能无失真的恢复出原始信号。图 2.6 从 x
15、1(n)和 x2(n)分别重构模拟信号试验三 IIR 数字滤波器的设计 MATLAB 实现目前，设计 IIR 数字滤波器的通用方法是先设计相应的低通滤波器，然后再通过双线性变换法和频率变换得到所需要的数字滤波器。模拟滤波器从功能上分有低通、高通、带通及带阻四种，从类型上分有巴特沃兹Butterworth滤波器、切比雪夫ChebyshevI 型滤波器、切比雪夫II 型滤波器、椭圆Elliptic滤波器以及贝塞尔Bessel滤波器等。1 应用MATLAB 设计IIR 数字滤波器相关文件下面给出与IIR 数字滤波器设计有关的MATLAB文件。1. buttord.m用来确定数字低通或模拟低通滤波器的
16、阶次，其调用格式分别是1N,Wn=buttord(Wp,Ws,Rp,Rs)2N,Wn=buttord(Wp,Ws,Rp,Rs,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，因此，
17、它们实际上是频率。2. buttap.m用来设计模拟低通原型滤波器G( p) ，其调用格式是z,p,k=buttap(N)N 是欲设计的低通原型滤波器的阶次，z、p 和k 分别是设计出的G( p) 的极点、零点及增益。3. lp2lp.m4. lp2hp.m5. lp2bp.m6. lp2bs.m以上 4 个文件用来将模拟低通原型滤波器G( p) 分别转换为低通、高通、带通、及带阻滤波器。其调用格式分别是1B,A=lp2lp(b,a,Wo)或B,A=lp2hp(b,a,Wo)2B,A=lp2bp(b,a,Wo,Bw)或B,A=lp2bs(b,a,Wo,Bw)式中 b，a 分别是模拟低通原型滤波
18、器G( p) 有分子、分母多项式的系数向量，B，A 分别是转换后的 H (s)有分子、分母多项式的系数向量；在格式1中，Wo 是低通或高通滤波器的截止频率；在格式2中Wo 是带通或带阻滤波器的中心频率，Bw 是其带宽。7. bilinear.m实现双线性变换，即由模拟滤波器H(s)得到数字滤波器 H(z)。其调用格式是( )Bz,Az=bilinear(B,A,Fs)( )式中B、A 分别是 H s 的分子、分母多项式的系数向量；Bz、Az 分别是 H z 的分子、分母多项式的系数向量，Fs 是抽样频率。8. butter.m用来直接设计巴特沃兹数字滤波器，实际上它把buttord.m，but
19、tap.m，lp2lp.m 及bilinear.m等文件都包含进去，从而使设计过程更简捷。其调用格式是1B,A=butter(N,Wn)2B,A=butter(N,Wn,high)3B,A=butter(N,Wn,stop)4B,A=butter(N,Wn,s)格式13用来设计数字滤波器，B、A 分别是 H (z)的分子、分母多项式的系数向量，Wn 是通带截止频率，范围在 01 之间，1 对应抽样频率的一半。假设Wn 是标量，则格式1用来设计低通数字滤波器。假设 Wn 是1 2 的向量，则格式1用来设计数字带通滤波器；格式2用来设计数字高通滤波器；格式3用来设计数字带阻滤波器。格式4用来设计模
20、拟滤波器。9. cheb1ord.m求切比雪夫I 型滤波器的阶次。10. cheb1ap.m用来设计原型切比雪夫I 型模拟滤波器。11. cheby1.m直接设计切比雪夫I 型滤波器。以上 3 个文件的调用格式与对应的巴特沃兹滤波器的文件类似。12. cheb2ord.m13. ellipord.m14. cheb2ap.m15. ellipap.m16. besselap.m17. cheby2.m18. ellip.m19. besself.m以上分别为切比雪夫II 型滤波器、椭圆及贝塞尔滤波器设计函数。其格式与切比雪夫I型滤波器和巴特沃兹滤波器类同。此外，与本章内容有关的MATLAB文件
21、还有：20. impinvar.m用冲激响应不变法实现w 到W 及 s 到 z 的转换。21maxflat.m设计广义巴特沃兹低通滤波器。22yulewalk.m利用最小平方方法设计Yule-Walker滤波器。例 1 利用巴特沃斯模拟滤波器设计数字高通滤波器，要求通带截止频率为0.6 ，通带内衰减不大于 1dB，阻带起始频率为 0.4 ，阻带内衰减不小于 15dB,，采样周期Ts=1s。%using butterworth to design highpass DF clear all;Wp=0.6*pi;Ws=0.4*pi;Ap=1;As=15;N,wn=buttord(Wp/pi,Ws/
22、pi,Ap,As);%计算巴特沃斯滤波器阶次和截止频率b,a=butter(N,wn,”high”);%频率变换法设计巴特沃斯高通滤波器、C,B,A=sdir2cas(b,a)%把滤波器的系统函数转换为成级联型db,mag,pha,w=freqz_m2(b,a);%数字滤波器的频率响应subplot(211); plot(w/pi,mag);grid; axis(0,1,0,1.2);title(”数字滤波器的幅频响应”); subplot(212); plot(w/pi,db);grid; axis(0,1,10,-200);title(”数字滤波器的幅频响应dB”);function C,
23、B,A=sdir2cas(b,a);%把滤波器的系统函数转换为成级联型%a-数字滤波器系统函数分母系数%b-数字滤波器系统函数分子系数%C-数字滤波器阶联构造的常数因子%B-数字滤波器阶联构造的分母系数%A-数字滤波器阶联构造的分子系数Na=length(a)-1;Nb=length(b)-1; b0=b(1);b=b/b0;a0=a(1);a=a/a0;C=b0/a0;p=cplxpair(roots(a);K=floor(Na/2); if K*2=NaA=zeros(K,3); for n=1:2:NaArow=p(n:1:n+1,:);Arow=poly(Arow); A(fix(n+
24、1)/2),:)=real(Arow);end elseif Na=1A=0 real(poly(p);elseA=zeros(K+1,3); for n=1:2:2*KArow=p(n:1:n+1,:);Arow=poly(Arow); A(fix(n+1)/2),:)=real(Arow);endendA(K+1,:)=0 real(poly(p(Na);z=cplxpair(roots(b);K=floor(Nb/2); if Nb=0B=0 0 poly(z); elseif K*2=NbB=zeros(K,3); for n=1:2:NbBrow=z(n:1:n+1,:);Brow=
25、poly(Brow); B(fix(n+1)/2),:)=real(Brow);end elseif Nb=1B=0 real(poly(z);elseB=zeros(K+1,3); for n=1:2:2*KendBrow=z(n:1:n+1,:);Brow=poly(Brow); B(fix(n+1)/2),:)=real(Brow);B(K+1,:)=0 real(poly(z(Nb);end例 2 : fp=0.3kHz，Rp=1dB，fs=0.2kHz,Rs=20dB，利用双线性变换法设计一个 Chebyshev I 型数字高通滤波，分析作畸变处理的效果。clear all;clos
26、e all; Rp=1;Rs=20;T=0.001;fp=300;fs=200;wp=2*pi*fp;ws=2*pi*fs;%预畸变wp1=(2/T)*tan(wp*T/2);ws1=(2/T)*tan(ws*T/2);%设计模拟滤波器,无预畸变处理n,wn=cheb1ord(wp,ws,Rp,Rs,”s”);b,a=cheby1(n,Rp,wn,”high”,”s”);%设计模拟滤波器,有预畸变处理n1,wn1=cheb1ord(wp1,ws1,Rp,Rs,”s”);b1,a1=cheby1(n1,Rp,wn1,”high”,”s”);%双线性变换bz,az=bilinear(b,a,1/T
27、);db,w=freqz(bz,az); plot(w/pi,20*log10(abs(db),”LineWidth”,4,”color”,”k”); axis(0,1,-30,2);grid on; hold on;bz1,az1=bilinear(b1,a1,1/T); db1,w1=freqz(bz1,az1);plot(w1/pi,20*log10(abs(db1),”LineWidth”,2,”color”,”r”); axis(0,1,-30,2);hold off;(取T=0.001s)例 3.利用脉冲响应不变法设计一个数字巴特沃思低通滤波器，通带截止频率fp=0.1kHz，通带
28、最大衰减 Rp=1dB，阻带截止频率 fs=0.3kHz，阻带最小衰减Rs=10dB。争论不同采样频率对所设计数字滤波器频率响应的影响。设采样频率 fs 分别取 1kHz，2kHz，4kHz。clear all;close all; wp=2*pi*100;ws=2*pi*300;Rp=1;Rs=10; n,wn=buttord(wp,ws,Rp,Rs,”s”);b,a=butter(n,wn,”s”);%求模拟滤波器的频率响应db,w=freqs(b,a,500*2*pi);plot(w/(2*pi),20*log10(abs(db),”LineWidth”,2,”color”,”k”);
29、grid;axis(0,500,-20,1);hold on;%脉冲响应不变法fs1=1000;bz,az=impinvar(b,a,fs1);%求数字滤波器的频率响应db,w=freqz(bz,az); plot(0.5*fs1*w/pi,20*log10(abs(db),”r”); axis(0,500,-20,1);hold on;%脉冲响应不变法fs2=2023;bz,az=impinvar(b,a,fs2);%求数字滤波器的频率响应db,w=freqz(bz,az); plot(0.5*fs2*w/pi,20*log10(abs(db),”b”); axis(0,500,-20,1)
30、;hold on; fs3=4000;bz,az=impinvar(b,a,fs3);%求数字滤波器的频率响应db,w=freqz(bz,az); plot(0.5*fs3*w/pi,20*log10(abs(db),”g”); axis(0,500,-20,1);title(”冲激响应不变法”); hold on;hold off;试验四 应用 MATLAB 设计 FIR 数字滤波器与本章内容有关的 MATLAB 文件与 FIR 数字滤波器设计有关的MATLAB文件主要分为两类。一类是用于产生各种窗函数，另一类是用于设计FIR 数字滤波器。1窗函数用于产生窗函数的MATLAB文件有如下 8
31、个：（1） bartlett.m(三角窗)（2） blackman.m(布莱克曼窗)（3） boxcar.m(矩形窗)（4） hamming.m(海明窗)（5） hanning.m(汉宁窗)（6） triang.m(三角窗)（7） chebwin.m(切比雪夫窗)（8） kaiser.m(凯泽窗) 2FIR 数字滤波器的文件以下 MATLAB文件是用于FIR 数字滤波器的文件：（1） fir1.m本文件承受窗函数法设计FIR 数字滤波器，其调用格式是1b=fir1N,Wn2b=fir1N,Wn,high 3b=fir1N,Wn,stop式中N 为滤波器的阶次，因此滤波器的长度为 N+1；(W)
32、n 是通带截止频率，其值在 01 之间，1 对应抽样频率的一半；b 是设计好的滤波器系数h n 。对于格式1，假设 Wn 是一标量，则可用来设计低通滤波器；假设 Wn 是1 2 的向量，则用来设计带通滤波器；假设 Wn 是1 L 的向量，则可用来设计L 带滤波器，此时，格式将变为：b=fir1N,Wn,DC-1或b=fir1N,Wn,DC-0 其中，前者保证第一个带为通带，后者保证第一个带为阻带。格式2用来设计高通滤波器；格式3用来设计带阻滤波器。值得留意是：在上述全部格式中，假设不指定窗函数的类型，则fir1 自动选择汉明窗。（2） fir2.m本文件承受窗函数法设计具有任意幅频特性的FIR
33、 数字滤波器。其调用格式是b=fir1N,F,M其中 F 是频率向量，其值在 01 之间，M 是与 F 相对应的所期望的幅频响应。不指定窗函数的类型，则自动选择汉明窗。（3） remez.m本文件用来设计承受切比雪夫最正确全都靠近 FIR 数字滤波器。同时，还可以用来设计希尔伯特变换器和差分器。其调用格式是1b=remezN,F,A 2b=remezN,F,A,W 3b=remezN,F,A,W,hilbert 4b=remezN,F,A,W,differentiator其中，N 是给定的滤波器的阶次；b 是设计的滤波器的系数，其长度为 N+1；F 是频率向量， 其值在 01 之间；A 是对应
34、F 的各频段上的抱负幅频响应；W 是各频段上的加权向量。值得留意是：假设 b 的长度为偶数，设计高通和带阻滤波器时有可能消灭错误，因此最好保证 b 的长度为奇数，即N 应为偶数。（4） remexord.m本文件承受切比雪夫全都靠近设计 FIR 数字滤波器时所需要的滤波器阶次。其调用格式是N,Fo,Ao,W=remexordF,A,DEV,Fs式中，F、A 的含义同文件3，是通带和阻带上的偏差；该文件输出的是符合要求的滤波器阶次 N、频率向量 Fo、幅度向量 Ao 和加权向量 W。假设设计者事先不能确定自己要设计的滤波器的阶次，那么，调用remexord 后，就可利用这一族参数再调用remez
35、，即b=remezN,Fo,Ao,W，从而设计出所需要的滤波器。因此，通常remez 和 remexord 结合使用。 值得说明是：remexord 给出的阶次N 有可能偏低，这时适当增加N 即可；另外，假设N为奇数，就可令其加 1，使其变为偶数，这样b 的长度为奇数。（5） sgolay.m本文件用来设计Savitzky-Golay 平滑滤波器。其调用格式是b=sgolayk,f式中 k 是多项式的阶次，f 是拟合的双边点数。要求k f ，且 f 为奇数。（6） firls.m本文件用最小平方法设计线性相位FIR 数字滤波器。可设计任意给定的抱负幅频特性。（7） fircls.m用带约束的最
36、小平方法设计线性相位 FIR 数字滤波器。可设计任意给定的抱负幅频特性。（8） fircls1.m用带约束最小平方法设计线性相位 FIR 低通和高通滤波器。可设计任意给定的抱负幅频特性。（9） firrcos.m用来设计低通线性相位FIR 数字滤波器，其过渡带为余弦函数外形。例 1、选择适宜的窗函数设计一个FIR 数字低通滤波器，要求：通带截止频率为wp0.2p，Rp0.25 dB；阻带截止频率为ws0.3p，As50 dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。% 利用海窗设计数字低通滤波器clear all; Wp=0.2*pi; Ws=0.4*pi;tr_w
37、idth=Ws-Wp;%过渡带宽度N=ceil(6.6*pi/tr_width)+1%滤波器长度n=0:1:N-1;Wc=(Ws+Wp)/2;%抱负低通滤波器的截止频率hd=ideal_lp(Wc,N);%抱负低通滤波器的单位抽样响应w_ham=(hamming(N)”;%海明窗h=hd.*w_ham;db,mag,pha,w=freqz_m2(h,1);%计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1)%实际通带纹波As=-round(max(db(Ws/delta_w+1:1:501) %实际阻带纹波figure(1);grid ;subplot(2,2,1);stem(n,hd);grid ;subplot(2,2,2); stem(n,w_ham);subplot(2,2,3);stem(n,h);subplot(2,2,4); plot(w/pi,db); ax