数字信号处理实验与课程设计教程

1、数字信号处理实验与课程设计教程戴虹编工学部计算机与信息工程学院 2015年12月内容简介：数字信号处理是高校通信工程、电子信息工程、自动化、测控等专业 的一门非常重要的专业基础课程，是我校的校重点课程、校精品课程与上海市教 委重点课程。本书是数字信号处理课程的实验及课程设计教材，基于matlab编 程，配合理论授课内容与应用型人才培养的教学目标而编写。全书分为上下两篇， 上篇为数字信号实验指导，共有5个实验：信号、系统及系统响应、离散信号与 系统的复频域分析、dft和fft频谱分析、iir数字滤波器设计、fir数字滤波 器设计。下篇是数字信号处理课程设计教程，共3个课程设计题：双音多频(dtm

2、f)的产生与检测、音乐信号中啸叫噪声的消除、回音消除。力求使学生 对数字信号处理课程产生浓厚的兴趣，并掌握数字信号处理的mat i ab仿真实践 方法。课程设计以小组方式进行，要求学生任选一题进行设计并撰写课程设计报 告。本书适合作为普通高校电子信息工程、通信工程、测控、自动化等专业本科 数字信号处理实验与课程设计教材，也可作为相关领域工程技术人员的参考书 籍。目录上篇数字信号处理实验指导 实验一信号、系统及系统响应 实验二离散信号与系统的复频域分析实验三dft和fft频谱分析实验四iir数字滤波器设计实验五fir数字滤波器设计下篇数字信号处理课程设计教程课题1 双音多频cdtmf）信号的产生

3、与检测 课题2音乐信号中的啸叫噪声消除 课题3 回音消除附录a matlab程序设计入门附录b实验报告撰写格式参考文献上篇数字信号处理实验指导 实验一信号、系统及系统响应一、实验目的1. 掌握典型序列的产生方法。2. 掌握dft的实现方法，利用dft对信号进行频域分析。3. 熟悉连续信号经采样前后频谱的变化，加深对时域采样定理的理解。4. 分别利用卷积和dft分析信号及系统的时域和频域特性，验证时域卷积定理。二、实验环境1. windows2000 操作系统2. matlab6.0三、实验原理1. 信号采样对连续信号xa(t)=ae"atsiii(qut)ii(t)进行采样，采样周期

4、为t,采样点0<n<50,得采样序列 xa(n)=o2. 离散傅里叶变换(dft)设序列为x(n)，长度为n，则n-x(ehk)=dft|x(n)|= e x(n) e_j°k，1,/i=0其中(|<=0，1，2，.，1-1)，通常m>n，以便观察频谱的细节。|x(ejwk)|-x(n)的幅频谱。 m4. 连续信号采样前后频谱的变化zh=-ooa即采样信号的频谱xa(jq)是原连续信号xa(t)的频谱xa(jq)沿频率轴，以周期重复出现，幅度为原来的倍。5. 采样定理由采样信号无失真地恢复原连续信号的条件，即采样定理为：6. 时域卷积定理设离散线性时不变系统输

5、入信号为x(ii)，单位脉冲响应为h(n)，则输出信号y(n)=；由时域卷积定理，在频域中,y(ei<°)=fty(n)=。四、实验内容1. 分析采样序列特性(1)程序输入产生采样序列 xa(n)= ae'an 1 sin( q «nt)u(n)(o彡n<50)，其中 a=44.128,a=50 72 71, q o=50 a/i 兀 采样频率 (可变)，t=l/fs。(要求写程序注释)%程序 shiyanll.mdear%clc%a=444.128;a=§0*sqrt(2)*pi;%w0=50*sqrt(2)*pi;fs=inputf输入采样

6、频率fs=f);t=l/fs;n=50;n=o:n-l;x a=a *exp(-a*n*t).*sin(wo*n*t); %subplot(221);stem(n，xa，.);grid; %m=100;xa,wkl=dft(xa,m);%f=wk*fs/(2 * pi);%subplot(222);plot(f,abs(xa);grid; %dft子函数:dft.mfunction |x,wk|=dft(x,m)n=length(x);%n=0:n-l;for k=0:m-lwk(k+l)=2*pi/m*k;x(k+l)=sum(x.*exp(-j*wk(k+l)*n); %end(2) 实验

7、及结果分析a. 取fs=1000(hz),绘出xa及ixal的波形。b. 取 fs=300(hz)，绘出 xa及|xa(ejwk)|的波形。c. 取 fs=200(hz)，绘出 xa及|xa(ejwk)|的波形。d. a，b，c中，哪几种情况出现了频谱混叠现象？ ；出现频谱混叠的原因是2. 时域离散信号和系统响应分析(1) hb(n)= 5 (n)+2.5 s (n-1 )+2.5 6 (n-2)+ 6 (n-3)程序语句为 hb=l，2.5,2.5，l;(2) 卷积语句：y=conv(x，h)其中x输入序列x(n)； h一单位脉冲响应h(n); y输出序列y(n)。3. 卷积定理验证(1)

8、编程实现 y(n)=xa(n)*hh(n),其中xa(n)= ae_antsin( q ont)u(n)(on<5g), a=l,a=0.4, q o=2.o734，t=lhb(n)= 6 (n)+2.5 5 (n-l)+2.5 5 (n-2)+ 6 (n-3)及 y(ejwk)=dfty(n)(m=100),(利用 dft.m);绘出波形。(2) 编程实现 xaykfdftixaoormmoo)及 hb(ejwk)=dfthb(n)j(m=100)；计算 y(ejwk)= xa(dwk) hb(ejuk)；绘出|y(ej"k)|波形。问：和(2)中的|y(ek)|波形一致吗？

9、 ；为什么？ 4. 正弦序列的产生设正弦序列x(n)=sin(wn)，取样频率fs=64hz，信号频率f=5hz, 样点n=0n-l，(n=64)，co=2n f/fs，编程产生x(n)并绘图。程序：%shiyanl41einclear clcn=o:l:n-l; w=2*pi*f/fs; x=sin(w*n); stem(x，.) 实验要求：n=64;fs=64;f=5; %1) 在后的空格内填入注释2)运行上述程序，绘出结果波形。3)设双音多频信号"2”为 x(n)=sin(<o in)+sin(a)2n)，其中 fl=697hz,f2=l336hz, 取样频率fs=800

10、0hz，编程产生x(n)并绘出x(n)的波形。(n=0799)5. 拓展题数字信号处理算法之时间平均时间平均可用来消除周期信号中的随机噪声。对m个周期的振动信号x(n)=s+u(n)其中s(n)-故障信号，为余弦序列； u-随机噪声|做时间平均，即把x(n)按周期n分段，将每个周期的对应点相加再做平均, 计算时间平均前后的信噪比。%时间平均程序programll. mclear% 定义n=10载入sip% 定义i从1至ijmclcn=10;m=input(m=); load sip x=sip;s=zeros(l，n); for i=l:ms=s+x(l+n*(i-l):i*n);end s=

11、s/m;k=0:n-ll;subplot(321);plot(k,x(l:n);grid;subplot(322);plot(k,s);grid;i=0:l:n-l;s0=cos(2*pi*i/n);ps=sum(so.a2)/n;pu=l;snro=l 0*logl 0(ps/pu) %原信噪比py=sum(s-so).a 2)/n;snr=10*logl0(ps/py)°/。时间平均后的信噪比%数据文件sip的产生程序shu.m clearclcn=10;m=1000;i=o:l:n-l;s=cos(2*pi*i/n);x=zeros(l，n*m);u=randn(size(l:

12、n*m);%生成 1 至j n*m 的矩阵for j=l:m% 定义 i 从 1 至lj mx(l+n*(j-l):j*n)=s+u(l+n*(j-l):j*n);end文件sipsave sip x -ascii% 产生数实验要求：1）在空格中填入注释。2）先运行shu.m,产生数据文件sip,再运行programll.m，分别绘出当m=10，100,500，10时， 输入/输出信号的波形。1o实验二离散信号与系统的复频域分析一、实验目的1. 掌握求解离散时间系统差分方程的两种方法：迭代法和filter函数法。2. 利用z变换对系统进行复频域分析。3. 掌握系统零极点的绘制方式及z反变换的求

13、解方法。4. 熟悉z变换的应用。二、实验环境1 .windows xp操作系统2. matlab2007a 软件三、实验原理3. 系统的输出y(n)与输入x(n)及单位脉冲响应h(n)的关系:4. 一个n阶常系数线性差分方程表达式：5.系统函数h与x，y的关系为：h(z)=。4. h(z)零点定义为：； h(z)极点定义为：5. 本实验涉及的matlab函数1) 求解差分方程：y=filter(b，a，x)其中：输入信号；y差分方程的解。b，a-输入/输出信号前的系数。2) 画零极点图:zplane(b，a，n)其中：b，ah(z)分子/分母系数。3) 求z反变换 r，p，k=residuez

14、(b，a)b-h(z)的分子多项式系数；a-h(z)的分母多项式系数；输出：r，p，k的含义如下：“residuez”可将h(z)分解为简单有理分式之和：h(z)l-/?(l)z1 - p(n)z+ (1) + a(2)z +. < 1 >式1中0(1)中(2).|(11)的集合为列向量0;1*(1)，1*(2).1*(11)的集合为列向量n k(l)，k(2).的集合为行向量k;己知r,p,k即可手工写出h(n)的表达式：wn)=r(l)|p(l)n + r(2)p(2)|n+".ilii)|p(ii)ru(ii)+k(l)wii)+k(2)6(n-l)+四、实验内容1

15、. 己知：系统的差分方程为：y(n)-y(n-l)+0.9v(n-2)=x(n)分别利用迭代法和filter函数法求此系统的单位脉冲响应h，并绘图。|1迭代法:%程序：programll.mclear clcn=100;x=|l zcros(l,n-l)|; %产生 x(n)=8(n)y(l)=0;%y(-2)=0y(2)=0;%y(-l)=0y(3)=i;%y(0)=lfor n=4:n%迭代法求解该差分方程，得y(n)y(n)=x(n)+y(n-l)-0.9*y(n-2);endk=-2:n-3;%时间轴范围kstem(k，y，.)；画 y(n)=h(n)titlef系统的单位脉冲响应h(

16、n)1); %写标题 h(n)波形：2filter 函数法：%程序：program 12.mclearclcn=100;x=l zeros(l，n-l); °/o产生 x(n)=6(n)b=l;%差分方程系数b，aa=l,-1,0.9；h=filter(b,a,x);%用filter函数求解系统的单位脉冲响应h(n)k=-2:n-3;stem(k，h，j);title(系统的单位脉冲响应h(n)；h(n)波形：2)修改输入信号为x(n)=u(n),编程实现：求系统的阶跃响应yl(n)。 u(n)由：u=ones(l，n)产 生i。%程序：program 13.ni3)利用由1)产生的

17、h(n)，利用卷积法产生系统的阶跃响应：y2(n)=n(n)*h(n)。%程序：program 14.my2(n)波形：问：y 1 (n)和y2(n)是否相同？ 2. 求解差分方程(迭代法及fdter函数法)，实现一个数字振荡器(正弦波信号): h(n)=sin(okn)，其中 ok=2nfk/fs，频率 fk=1000hz，取样频率 fs=10khz。|提示：h(n)对应的差分方程为：h(n)-ah(n-l)+h(n-2)=b6(n-l)<1>a=2cos(ok)=1.618,b=sin(ok)=0.588画出h并用“zplane”画该系统的零极点图。%迭代法实现程序progra

18、m21.mh(n)波形：o/ofilter函数法实现程序 program22.m问：极点位置处于单位圆(内/上/外)？3.求下列序列的z变换并用zplane”画出零极点分布图。 l)x(n)=l，2，l，3，该序列的 z 变换为：%程序 programsl.mclearclcb=l 2 1 3h°/«h分子系数 a=|1 0 0 0;%h(z)分母系数 zplane(b，a);%画零极点分布图 结果图形：2) x(n)=0.2nu(n)-u(n-5)|，贝！):x(z)=%程序 program32.m结果图形：4. z反变换473 -2172 +1s7用”residuez”

19、函数求：h(z) =；的z反变换。(z -5z + 6)(z-l)%程序 program4.mclear clcb=4 -21 18 0|;%h分子系数 a=l -6 11 -6;%h分母系数 r，p，k=residuez(b，a)°/<z 反变换中的系数 运行该程序，得：p=k=贝lj h(n)=实验三dft和fft频谱分析一、实验目的1. 掌握dft频谱分析的原理与编程方法。2. 理解fft算法的编程思想。2. 熟练掌握利用fft对信号作频谱分析，包括正确地进行参数选择、画频谱 及读频谱图。3. 利用fft频谱分析进行快速卷积和太阳黑子周期性检测。二、实验环境1. wind

20、ows xp以上操作系统2. 安装 matlab2007a 软件三、实验原理1.离散傅里叶变换(dft)设序列为x(n)，长度为n，则n-xcevdftwn)卜x(n) e如k n，/r=0q jr其中ok= a(k=0，l，2，.，m-l)，通常m>n，以便观察频谱的细节。|x(ejwk)卜-x(n)的幅频谱。 m2. 谱分析参数选择1) 设信号x(t)最髙频率为fc，对其进行取样得x(n)，根据取样定理，取样频率fs必须满足： o2) 设谱分辨率为f，则最小记录时间tpmin=；取样点数n>为使用快速傅里叶变换(fft)进行谱分析，n还须满足：n=。3. 用fft计算信号x(n

21、)的频谱。|设x(n)为实信号快速傅里叶变换(fft)是dft的一种快速算法，其使得dft的运算速度大为加快。1)对信号x(n)作n点fft，得频谱x(k)(k=on-1)x(k)=xr(k)+jx,(k) (k=0n/2-1), xr(k) x(k)的实部；xkk) x(k)的虚部。matlab语句：y=fft(x，n)其中：x-x(n);yx(k)2)幅频谱:|x(k)|= ，由于x为实信号，因此|x(k)| 对称,matlab语句：abs(y) iii)功率谱:psd(k)=|x(k)|2/n=x(k)x*(k)/n matlab语句：psd=y.*coni(y)/n 其中：conj(y

22、)-x'x(k)的共轭14. 读频谱图频谱图中任意频率点k对应实际频率为:fk=。5. 用fft实现线性卷积运算用fft实现y(n)=x(n)*h的步骤为：1) 设x(n)及h的长度分别为和乂。为使循环卷积等于线性卷积，用补0的方法使x(n)，h长度均为n,则n须满足伦；为用fft计算dft，则n还须满mn=o2) 用 fft计算x(k)，h(k) (n点)。3) y(k)=； v(n)=o四、实验内容1.根据公式设计dft原理程序，并计算：x(n)=|l，l，l，ll的4，16,64点dft并绘图。 %dft/idft程序 dft.m clc%输入序列x(n)=l 1 1 1 %x(

23、n)的长度m %变换区间n %补0,使xn长度为n%旋转因子wn%作(11)的 dft=xkclearxn=input(x(n)= )； m=length(xn);n=input(变换区间=); xn=xn zeros(l，n-m)|; n=0:n-l;k=0:n-l; nk=n，*k;wn=exp(-j*2*pi/n); wnk=wn八 nk; xk=xn*wnk subplot(211);stem(k，abs(xk)，'.);grid on; %显示xk的幅频谱(离散曲线)subplot(212);plot(k，abs(xk);grid on;%显示xk的幅频谱(连续曲线)运行结果

24、：问：由此得出怎样的结论？2.理解dit-fft算法原理程序，并用它计算x0<)=fft|r4(n)i，分别取n=4,8，16 和64,绘出幅频谱|x(k)|。%程序dit.m clearx=i n p u t( x= ) ；%n=input(fn= f);%x(length(x)+l:n)=zeros(l,n-iength(x); %补0 x(l:n) l=iog2(n);xl=zeros(l，n);forjl=l:n% 倒序xl(jl)=x(bin2dec(fliplr(dec2bin(jl-l,l)+l);end %fft(dit)%m=2;while(m<=n)w=exp(

25、-2*j*pi/m);v=l;for k=0:l:m/2-l for i=0:m:n-lp=k+i;q=p+m/2;a=xl(p+1);b=xl(q+l)*v;xl(p+l)=a+b;xl(q+l)=a-b;endv=v*w；endm=2*m;end%旋转因子w%k为每级蝶形运算旋转因子的个数 %1为各群的首序号%本级蝶形运算，xl最终存放x(k)%旋转因子w的变化%第11级subplot(211);stem(x，.);grid on; title(x(n);subplot(212);stem(abs(xl)/et);grid on;%title(,|x(k)|,);%3.fft谱分析设信号为

26、x(t)=sin(27rf1t)+sin(2nf2t)+随机噪声，fi=50hz, f2=120hz，以取样频率fs=lkhz对x(t)进行取样，样本长度tp=0.25s，得x(n)，对x(n)作256点fft，得频谱x(k)，画原信号 x(n)，幅频谱|x(k)|以及功率谱psd(k)，对信号进行谱分析。%程序 pufenxi.m clearfs=looo;t=0:l/fs:0.25;n=256;fl=50;f2=120;s=sin(2*pi*fl*t)+sin(2*pi*f2*t);x=s+randn(size(t);y=fft(x,n);psd=y.*conj(¥)/n; f=

27、fs/n*(0:n/2-l); subplot(3 ll);plot(x);%信号+噪声x(n)%clcsubplot(312);plot(f，abs(y(l:n/2);subplot(313);plot(f,psd(l:n/2);%1) 画出图形窗口显示的图形，并注名每个图形的含义。2) 回答下列问题：i) 观察幅频谱图，可以发现，信号x(ii)含有的两个频率分量分别是riz 和iizoii) 在该程序中的“f=fs/n*(0:n/2-l);”下添加“k=0:n/2-l;”,“plot(f，abs(y(l:n/2);，改为“plot(k，abs(y(l:n/2);，重新运行该程序并观察幅频谱

28、图，图中两峰值对应的下标分别是和o它们的含义为o再将该程序中的n改为512,重新运行该程序并观察幅频谱图，这时图中两峰值对应的下标分别是和。结果是否和上面的相同？ 为什么？ .hz，改变f2=60hz，问：在幅频谱中，能否分辨f,和f2为什么？iii) 本例的频谱分辨率f是. 对应的频率分量？ 再改变f2=52hz，问：在幅频谱中，能否分辨ff2对应的频率分量？为什么？ ohz；再改变f2=600hz，在幅频谱中，f2对应的频率分量出现在 问：在ft=1000hz的情况下，能否正确检测f2对应的频率分量？为什么？ ohzo在该程序中改变fs,验为了正确检测f2对应的频率分量，则fs至少取多少h

29、z?证你的结论。iv)比较幅频谱和功率谱，可以发现功率谱具有的特性。4. fft实现任意两个序列的快速卷积。%程序 fftjuanji.mclear%序列xl(n)，x2的长度 %ceil向+oc方向取整 %clcxl=input(xl=);x2=input(x2=); n l=length(xl);n2=length(x2); e=ceil(log2(nl+n2-l);n=2 八 e;xl=xl,zcros(l,n-nl)|;x2=lx2,zeros(l,n-n2)j;xl=fft(xl,n);x2=fft(x2,n);y=x1.*x2;y=ifft(¥,n)结果分析：1) 回到m

30、atlab窗口，键入：xl=l 1 l,x2=l 2,回车。结果:y=2) 问：可用matlab中的什么函数验算上述卷积结果？ 5. 利用谱分析观察太阳黑子周期性。以100年中记录到的太阳黑子出现次数为信号x(n)，对x(ii)作功率谱，从中观察太阳黑子周期 性。%程序 taiyangheizi.m clearclcx=101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 .132 131 118 90 67 60 47 41 21 16 64 7 1434 45 43 48 .4228 10820 1 5 12 14 35 46 41 30

31、 24 16 7 4 2 8 17 36 50 62 67 71 48 28 8 13 57 122 138 103 86 63 37 24 .11 15 40 62 98 124 96 66 64 54 39 21 7 4 23 55 94 96 .77 59 44 47 30 16 7 37 74|;%100年中太阳黑子出现的次数subplot(211);plot(x)%画又(11)n=128; fs=l;%fs=lhz，n=128 点s=x-mean(x);%对*作零均值化处理(去除直流分量)y=;%对8做 n 点 fftpsd=;%做功率谱psdf=；%将频率定标为实际频率fsubpl

32、ot(212);%画功率谱(n/2点)1) 填写空格中的画图语句并绘出结果图形。2) 从s的功率谱观察到，其幅度最高处对应的横坐标f=hz,则太阳黑子每隔年出现一次最高峰。3) 在对s做fft时，为何可取fs=lhz，n=128点？实验四iir数字滤波器设计一、实验目的(1) 掌握用双线性变换法设计iir数字低通和髙通滤波器。(2) 设计低通滤波器对实际心电图信号进行滤波。(3) 设计低通滤波器对含有啸叫噪声的音乐信号进行消噪。(4) 设计iir数字低通和高通滤波器对某个dtmf(双音多频)信号进行频带分离。二、实验环境1. windows98以上操作系统2. 安装matlab6.0以上版本三

33、、实验原理1. 选频型数字滤波器的种类有、和滤波器。2. 从实现方法上，数字滤波器通常分为和滤波器。3. iir滤波器的设计目的是根据技术指标，找到； iir滤波的matlab 语句为 y=;4. 双线性变换法设计低通iir滤波器的步骤是qj(2)(3) 四、实验内容1. 人体心电图信号在测量过程中往往受到工业高频干扰，必须经过低通滤波处理后才能作 为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(n),其中存在高频干 扰。试以x(n)作为输入序列，滤除其中的干扰成分。x(n)= -4， -2， 0， -4， -6， -4， -2， -4， -6， -6, -4， -4， -6，

34、-6， -2， 6， 12， 8， 0, -16, -38, -60， -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0,0,0,0, -2, -4, 0, 0, 0, -2, -2, 0,0, -2, -2, -2, -2, 0。低通滤波器设计指标：o p=0.2 h rad,s=0.3 rad，ap=ldb，as=15dbofj(z) =0.0007378(1+ z'1)6(1 -1.268z_, + 0.705z-2 )(1 -1.0106z_, + 0.3583z-2 )(1 - 0.904z_1 +

35、 0.215z己设计出h(p300)3=n h相k=hk(z)a(l + 2z-1 + z-2) 1_5az_1 _gz2其中a=0.09036 ； bfl.2686,<=0.7051 ; b2=1.0106,c2=-0.3583； b3=0.9044,c3=-0.2155 *iir 滤波的 mat lab 语句：y=f liter (b, a, x)b, ahk(z)分子/分母系数；x输入信号x (n); y滤波结果y (n)。* 求频响特性的 mat lab 语句：h, w=freqz (b, a, n)b,ahk(z)分子/分母系数；n频率点数；h滤波器的频响特性h(w);w数字频

36、率轴w。用iir低通滤波器(原理设计)对实际的心电信号进行滤波(iirecg.m)%程序：(iirecg.m) clearclc-4,-2,0,-4,-6, -4, 2, -4, -6, 6,-4,-4,-6,-6, -2, 6,12, 8, 0,-16-38,-60,-84,-90, -66, -32, -4, -212,12,10,6, 6, 6, 4, 0, 0,0,o,-2,-4, 0, 0, 0, -2, -20,一 2,一 2,-2, -2, 0;%心电信a=0. 09036;%iir低通滤波器系数b1=1.2686;c1=-o.7051;b2=l. 0106;c2=-0.3583

37、;b3=0. 9044;c3=-0. 2155;b=a 2*a a;al=l -bi -cl;n=128;hl, w=freqz (b, al, n) ;%频响特性 hl (w)yl=f ilter (b, al, x) ;%yl (n)是 x (n)经 hl (z)滤波的结果a2=l -b2 -c2;h2, w =freqz (b, a2, n); %频响特性 h2 (w)y2=filter(b, a2, yl):%y2(n)是 yl (n)经 h2(z)滤波的结果a3=l -b3 -c3;h3, w=freqz (b, a3, n) ;%频响特性 h3 (w)h=h1. *h2.*h3;%

38、总的频响特性 h (w) =h1 (w) h2 (w) h3 (w)mag=abs (h);%滤波器的幅频特性db=20*logl0( (mag+eps) /max (mag) ;%幅频特性(db)y3=f ilter (b, a3, y2);%y3(n)是 y2(n)经 h3(z)滤波的结果x=fft (x, n) ;%对 x (n)做 n 点 fft,得其频谱 xwx=2*pi* (0:n/2-1) /n; %将坐标轴从频率点k转换为数字频率wx (wx=2*pi*k/n)%绘图%subplot(221);plot(x):grid on;titlec x(n)；subplot(222):p

39、lot(wx/pi, abs(x(l:n/2);grid on;title( |x(wx)|)；subplot (223);plot (w/pi, db) ;grid on;title( |h(w) | (db)；subplot (224);plot (y3) ;grid on; title( y3(n)；要求：绘出并分析结果图形。(2)用现有的matlab函数设计上述iir低通滤波器，对实际的心电信号进行滤波，同样绘出 题1的结果图形。matlab 函数：*根据技术指标,os ,<xp ,<zs计算巴特沃思滤波器的阶数n和3db截止频率语句： n,wc=buttord(wp/pi

40、,ws/pi,ap,as)注：数字频率以n为单位。*根据n, co。确定数字滤波器h (z)分子/分母多项式的系数b, a低通b, a =butter (n, wc)髙通b, a =butter (n, wc, high)%程序 iirecgbt. mclearclcx= -4, 2, 0, 4, 6, 一4, 2, 4, 6, -6,.一4，一4，一6，一6, 2, 6, 12, 8, 0, 16，.一38, -60, 84, -90, 一66, -32, 4, _2, _4，8,.12, 12, 10, 6, 6, 6, 4, 0, 0, 0,.0, 0, -2, -4, 0, 0，0,

41、-2，-2, 0,.0, -2, -2, -2, -2, 0;%心电信号x(n)含髙频噪声wp=0. 2*pi;ws=0. 3*pi;ap=l;as=15;%低通滤波器技术指标n=128;nl, wc=buttord(wp/pi, ws/pi, ap, as)%b, a =butter (nl, wc) %h, w=freqz (b, a, n):%y3=fliter (b, a, x);%mag=abs(h);%db=20*logl0( (mag+eps) /max (mag) ;%幅频特性(db)x=fft (x, n) ;%wx=2*pi*(0:n/2-l)/n; %绘图%subplot

42、 (221) ;plot (x):grid on; title(，x(n”)；subplot(222);plot(wx/pi, abs(x(l:n/2);grid on;title( |x(wx)|)；subplot(223);plot(w/pi, db);grid on;title(，|h(w)|(db)；subplot (224);plot (y3) ;grid on; title( y3(n)；1) 在空格中填写程序注释。2) 运行上述程序，得：滤波器阶数n=； 3db截止频率wc=；b=； a=3)修改as=50，得：n=； wc=；b=;2. 文件”yinyue.wav”是含有啸叫噪

43、声的音乐信号，首先对其进行谱分析，观察信号和噪声的 频带范围，再设计适当的iir数字低通滤波器对其进行消噪处理恢复原信号，将结果存入 音乐文件” yinyuexiaozao.wav”并用耳机监听消噪前后的音乐信号。% 程序:iirxiaozaoy.nlclear%做fft的点数%读入含噪音乐信号%播放yl %暂停10sclcnl=81920;yl=wavtead(yinyue.'vav);sound(yl)yl=yl*6;yl=fft(yl,nl);psdl=yl.*conj(yl)/nl;fs=8192;f=8192/nl*(0:nl/2-l); % 绘图 °/。°

44、;/。°/。pause(lo)%figure(l)subplot(211);plot(real(yl);grid on;subplot(212);plot(f,psdl(1 :nl/2);axis(0 5000 0 500);grid on;实验要求：(1) 运行上述程序，绘出结果波形并监听含噪音乐信号。(2) 从第二张子图可见音乐信号的频带范围在hz，啸叫噪声频率为hz。(3) 根据后的要求编写程序，将这些程序连在上述程序之后。%iir低通滤波器%fp=；%通带截止频率fp=1000hzfsl=;%阻带截止频率fsl=1200hzwp=2*pi*fp/fs;%化为数字频率wp和ws

45、ws=;rp=;as=; %rp=l db,as=50dbn,wn=；%根据技术指标计算巴特沃思滤波器的阶数n和3db截止频率o。|b,a| =；%根据n, co e确定数字滤波器h (z)分子/分母多项式的系数b, a%h，w=freqz(b，a，fs，whole)；频响特性 h(w)mag=；°/。幅频特性db=; %幅频特性(以db为单位)y=;%用所设计的滤波器对yl进行低通滤波消噪，得y % 绘图。/。°/。°/。°/。figure(2)subplot(211);plot(w*fs/(2*pi),db);axis(|0 fs/2 -400 10

46、0);grid on;subplot(212);plot(real(y);grid on;sound(real(y)%播放 y叭狄1¥1*加(代31(丫)，41113(30.外3'');将结果存入音乐文件：yinyuexiaozao.wav”(4) 绘出结果波形并监听己消噪的音乐信号。'3. dtmf(双音多频)信号是按键电话拨号音，其中”2”键是由下列低频音和髙频音所构成： x(n)=sin(2 h fvfsxnh sin(2 h f2/fsxn)，fi=697hz，f2=1336hz,取样频率 fs=8000hzo (n=0 1599)编写程序，完成以下功

47、能：(1) 产生x(n)并绘图，用”sound”语句播放x(n)，用耳机监听，并将x(n)存入音乐文 件 ”dtmf2.wav”。(2) 设计适当的低通和髙通滤波器对x(n)进行滤波，将低频音和高频音分离，得分离后的信号: x(n)=sn(2 h fj/fsxn), x2(n)=sin(2 冗 f2/fs x n),画出 xi(n)和 x2(n),并分别将 xn)和 x2(n)存入 音乐文件”dtmfdi.wav”和”dtmfgao.wav”，用耳机监听这两种声音。%程序：dtmf.m波形:实验五hr数字滤波器设计一、实验目的1. 掌握用窗函数法设计fir低通数字滤波器。2. 了解各种窗函数对

48、滤波特性的影响。3. 设计f1r低通滤波器对含噪音乐信号进行消噪。4. 用二维低通fir滤波器连接图像中文字的断裂部分。二、实验环境1. windows98以上操作系统2. 安装matlab6.0以上版本三、实验原理1. fir滤波器的设计目的是根据技术指标，设计。fir滤波的matlab语句为v=o2. 窗函数矩形窗的表达式为：(2)汉宁窗的表达式为：哈明窗的表达式为：四、实验内容1. 人体心电图信号在测量过程中往往受到工业髙频干扰，必须经过低通滤波处理后才能作 为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(ii),其中存在高频干 扰。试以x(n)作为输入序列，用窗函数法设计

49、一个fir低通滤波器滤除其中的干扰成分。 x(n)= -4, -2,0, -4, -6, -4, -2,-4, -6, -6, -4, -4,-6, -6, -2, 6,12, 8,0,-16, -38, -60,-84， -90, -66， -32,-4， -2， -4， 8, 12，12， 10， 6， 6,6， 4，0,0,0,0,0, -2，-4，0，0，0，-2，-2，0,0，-2，-2，-2，-2，0。低通滤波器设计指标：<*> p=0.2 n rad, s=0.3 h rad，ap=ldb，as=50db。选哈明窗，即 w(n)=0.54-0.46cos(2 n/(n

50、-l)rn(n)o*fir 滤波的 mat lab 语句：y=conv (h, x)h滤波器单位脉冲响应h(n); x 输入信号x(n); y 滤波结果y(n)。程序：%用fir低通滤波器对实际的心电信号进行滤波(firecg. m)clearclcwp=o. 2*pi;ws=o. 3*pi;%tr=wswp;%n=cei1(8*pi/tr)+l %n=o:l:n-l;wc=(ws+wp) /2;%m=n-(n-l)/2+eps;hd=sin(wc*m). /(pi*m);%w_ham=(0. 54-0. 46*cos(2*pi*n/(n-1);%h=hd. *w_ham;%h, w=freq

51、z (h, 1, 1000,，whole,);%mag=abs(h);db=20*logl0(mag+eps)/max(mag);%delta_w=2*pi/1000;ap=- (min (db (1:1: wp/delta_w+l)%技术指标验证 ap, asas=-round(max(db(ws/delta_w+l:1:501)%绘图%figure (1)subplot (221) ;stem(n, hd，. );grid on;title(，hd(n)；subplot (222):stem(n, w_ham, . ):grid on;title( 哈明窗 w_ham(n); subplo

52、t(223);stem(n, h, );grid on;titlec h(n); subplot (224) :plot (w(l :501)/pi, db(1:501):grid on;title( h(w); %-4, 8, 0,0,.4, 一4, _6, -6, 2, 6,12, 8, 0, -16,一38, 一60, 84, -90, -66, 一32, 4, -2,12, 12, 10, 6, 6，6, 4, 0, 0,0, 0, -2，4, 0, 0，0, 2, -2,y=conv (x, h);x=f f t (x, 128);y=f f t (y, 128); wx=2*pi*

53、(0:n/2-l)/n; %绘图%0, -2, -2, 2, 2, 0;% 含噪的 ecg 信号%figure (2)subplot(221);plot(x);axis(0 60 -100 20);grid on;title(x(n);subplot(222):stem(wx/pi, abs(x(l:n/2), . ):grid on;title(，|x(w)|)；subplot(223):plot(y);axis(n-l)/2 (n-l)/2+60 -100 20);grid on;titlec y(n)*) subplot(224) ;stem(wx/pi, abs(y(l:n/2), . );grid on;title(，|y(w) |)；%绘图%要求：1) 绘出结果图形，并观察图2的第1和第3张小图，可见，y(n)和x(n)相比，延迟了 点。2) 运行上述程序，得：滤波器阶数n=；实际的ap=；实际的as=;可见，在设计相同技术指标的滤波器