数字信号处理实验-滤波器设计

数字信号处理实验指导书（修订版1）PAGE1数字信号处理（本科）实验指导书（修订版1）戴虹编2010年6月目录实验一信号、系统及系统响应实验二FFT频谱分析实验三IIR数字滤波器设计实验四FIR数字滤波器设计实验五离散系统的时域及复频域分析综合实验双音多频（DTMF）通信设计的MATLAB仿真附录实验报告格式实验一信号、系统及系统响应一、实验目的掌握典型序列的产生方法。掌握DFT的实现方法，利用DFT对信号进行频域分析。熟悉连续信号经采样前后频谱的变化，加深对时域采样定理的理解。分别利用卷积和DFT分析信号及系统的时域和频域特性，验证时域卷积定理。二、实验环境Windows2000操作系统MATLAB6.0三、实验原理信号采样对连续信号xa(t)=Ae-atsin(Ω0t)u(t)进行采样，采样周期为T,采样点0≤n<50，得采样序列xa(n)=。离散傅里叶变换(DFT)设序列为x(n),长度为N,则X(ejωk)=DFT[x(n)]=x(n)e-jωkn,其中ωk=(k=0,1,2,…,M-1),通常M>N,以便观察频谱的细节。|X(ejωk)|x(n)的幅频谱。4.连续信号采样前后频谱的变化a(jΩ)=即采样信号的频谱a(jΩ)是原连续信号xa(t)的频谱Xa(jΩ)沿频率轴，以周期重复出现，幅度为原来的倍。采样定理由采样信号无失真地恢复原连续信号的条件，即采样定理为：。6．时域卷积定理设离散线性时不变系统输入信号为x(n),单位脉冲响应为h(n),则输出信号y(n)=；由时域卷积定理，在频域中，Y(ejω)=FT[y(n)]=。四、实验内容1.分析采样序列特性(1)程序输入产生采样序列xa(n)=Ae-anTsin(Ω0nT)u(n)(0≤n<50),其中A=44.128,a=50,Ω0=50采样频率fs（可变），T=1/fs。(要求写%程序注释)程序shiyan11.mclear%clc%A=444.128;a=50*sqrt(2)*pi;%w0=50*sqrt(2)*pi;fs=input('输入采样频率fs=');T=1/fs;N=50;n=0:N-1;xa=A*exp(-a*n*T).*sin(w0*n*T);%subplot(221);stem(n,xa,'.');grid;%M=100;[Xa,wk]=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-1;fork=0:M-1wk(k+1)=2*pi/M*k;X(k+1)=sum(x.*exp(-j*wk(k+1)*n));%end(2)实验及结果分析取fs=1000(Hz),绘出xa(n)及|Xa(ejωk)|的波形。b.取fs=300(Hz),绘出xa(n)及|Xa(ejωk)|的波形。c.取fs=200(Hz),绘出xa(n)及|Xa(ejωk)|的波形。d.a,b,c中，哪几种情况出现了频谱混叠现象？；出现频谱混叠的原因是。2.时域离散信号和系统响应分析(1)hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)程序语句为hb=[1,2.5,2.5,1];(2)卷积语句：y=conv(x,h)其中x输入序列x(n)；h单位脉冲响应h(n);y输出序列y(n)。3.卷积定理验证(1)编程实现y(n)=xa(n)*hb(n),其中xa(n)=Ae-anTsin(Ω0nT)u(n)(0≤n<50),A=1,a=0.4,Ω0=2.0734,T=1hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)及Y(ejωk)=DFT[y(n)](M=100)，(利用DFT.m)；绘出|Y(ejωk)|波形。(2)编程实现Xa(ejωk)=DFT[Xa(n)](M=100)及Hb(ejωk)=DFT[hb(n)](M=100)；计算Y(ejωk)=Xa(ejωk)Hb(ejωk)；绘出|Y(ejωk)|波形。问：(1)和(2)中的|Y(ejωk)|波形一致吗？；为什么？。4.正弦序列的产生设正弦序列x(n)=sin(ωn),取样频率fs=64Hz,信号频率f=5Hz,样点n=0~N-1,(N=64),ω=2πf/fs,编程产生x(n)并绘图。程序：%shiyan141.mclearclcN=64;fs=64;f=5;%n=0:1:N-1;%w=2*pi*f/fs;%x=sin(w*n);%stem(x,'.')%实验要求:在%后的空格内填入注释2)运行上述程序，绘出结果波形。3)设双音多频信号"2"为x(n)=sin(ω1n)+sin(ω2n),其中f1=697Hz,f2=1336Hz,取样频率fs=8000Hz,编程产生x(n)并绘出x(n)的波形。(n=0~799)5.拓展题数字信号处理算法之一时间平均时间平均可用来消除周期信号中的随机噪声。对m个周期的振动信号x(n)=s(n)+u(n)[其中s(n)--故障信号，为余弦序列；u(n)随机噪声]做时间平均,即把x(n)按周期n分段，将每个周期的对应点相加再做平均，计算时间平均前后的信噪比。%时间平均程序program11.mclearclcn=10;%m=input('m=');%loadsip%x=sip;s=zeros(1,n);fori=1:m%s=s+x(1+n*(i-1):i*n);ends=s/m;k=[0:n-1];subplot(321);plot(k,x(1:n));grid;subplot(322);plot(k,s);grid;i=0:1:n-1;s0=cos(2*pi*i/n);ps=sum(s0.^2)/n;pu=1;snr0=10*log10(ps/pu)%原信噪比py=sum((s-s0).^2)/n;snr=10*log10(ps/py)%时间平均后的信噪比%数据文件sip的产生程序shu.mclearclcn=10;m=1000;i=0:1:n-1;s=cos(2*pi*i/n);x=zeros(1,n*m);u=randn(size(1:n*m));%forj=1:m%x(1+n*(j-1):j*n)=s+u(1+n*(j-1):j*n);endsavesipx-ascii%实验要求：1)在空格中填入注释。2)先运行shu.m，产生数据文件sip,再运行program11.m,分别绘出当m=10,100,500,1000时，输入／输出信号的波形。实验二FFT频谱分析一、实验目的1.理解FFT算法的编程思想。2.熟练掌握利用FFT对信号作频谱分析。包括正确地进行参数选择、作频谱图以及读频谱图。3.了解FFT的应用。二、实验环境1.Windows98以上操作系统2.安装MATLAB6.0以上版本三、实验原理谱分析参数选择1)设信号x(t)最高频率为fc,对其进行取样得x(n),根据取样定理，取样频率fs必须满足：。2)设谱分辨率为F,则最小记录时间tpmin=；取样点数N≥;为使用快速傅里叶变换(FFT)进行谱分析，N还须满足：N=。2.用FFT计算信号x(n)的频谱。[设x(n)为实信号]1)对信号x(n)作N点FFT,得频谱X(k)(k=0~N-1)X(k)=XR(k)+jXI(k)(k=0~N/2-1)，XR(k)—X(k)的实部；XI(k)—X(k)的虚部。Matlab语句：Y=fft(x,N)其中：xx(n)YX(k)2)幅频谱:|X(k)|=，由于x(n)为实信号，因此|X(k)|对称，Matlab语句：abs(Y)iii)功率谱:PSD(k)=|X(k)|2/N=X(k)X*(k)/NMatlab语句：PSD=Y.*conj(Y)/N其中：conj(Y)--X*(k)[X(k)的共轭]3．读频谱图频谱图中任意频率点k对应实际频率为:fk=。三、用FFT实现线性卷积运算用FFT实现y(n)=x(n)*h(n)的步骤为：1)设x(n)及h(n)的长度分别为N1和N2。为使循环卷积等于线性卷积，用补0的方法使x(n),h(n)长度均为N，则N须满足N≥；为用FFT计算DFT,则N还须满足N=。2)用FFT计算X(k),H(k);(N点)3)Y(k)=。4)y(n)=。四、实验内容理解DIT-FFT算法原理程序，并用它计算X(k)=FFT[R4(n)],分别取N=4,8,16和64,绘出幅频谱|X(k)|。%程序DIT.mclearclcx=input('x=');%N=input('N=');%x(length(x)+1:N)=zeros(1,N-length(x));%补0x(1:N)l=log2(N);x1=zeros(1,N);forj1=1:N%倒序x1(j1)=x(bin2dec(fliplr(dec2bin(j1-1,l)))+1);end%%FFT(DIT)%%M=2;while(M<=N)W=exp(-2*j*pi/M);%旋转因子WV=1;fork=0:1:M/2-1%k为每级蝶形运算旋转因子的个数fori=0:M:N-1%i为各群的首序号p=k+i;q=p+M/2;A=x1(p+1);B=x1(q+1)*V;x1(p+1)=A+B;%本级蝶形运算，x1最终存放X(k)x1(q+1)=A-B;endV=V*W;%旋转因子W的变化endM=2*M;%第M级end%%%%%%subplot(211);stem(x,'.');gridon;%title('x(n)');%subplot(212);stem(abs(x1),'.');gridon;%title('|X(k)|');%2.FFT谱分析设信号为x(t)=sin(2πf1t)+sin(2πf2t)+随机噪声，f1=50Hz,f2=120Hz,以取样频率fs=1kHz对x(t)进行取样,样本长度tp=0.25s,得x(n),对x(n)作256点FFT,得频谱X(k),画原信号x(n),幅频谱|X(k)|以及功率谱PSD(k),对信号进行谱分析。%程序pufenxi.mclearclcfs=1000;t=0:1/fs:0.25;%N=256;%f1=50;f2=120;%s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%x=s+randn(size(t));%信号+噪声x(n)Y=fft(x,N);%PSD=Y.*conj(Y)/N;%f=fs/N*(0:N/2-1);%subplot(311);plot(x);%subplot(312);plot(f,abs(Y(1:N/2)));%subplot(313);plot(f,PSD(1:N/2));%1)画出图形窗口显示的图形，并注名每个图形的含义。2)回答下列问题：i)观察幅频谱图，可以发现，信号x(n)含有的两个频率分量分别是Hz和Hz。ii)在该程序中的“f=fs/N*(0:N/2-1);”下添加“k=0:N/2-1;”，“plot(f,abs(Y(1:N/2)));”改为“plot(k,abs(Y(1:N/2)));”重新运行该程序并观察幅频谱图，图中两峰值对应的下标分别是和。它们的含义为。再将该程序中的N改为512,重新运行该程序并观察幅频谱图，这时图中两峰值对应的下标分别是和。结果是否和上面的相同？为什么？。iii)本例的频谱分辨率F是Hz,改变f2=60Hz,问：在幅频谱中，能否分辨f1和f2对应的频率分量？。为什么？。再改变f2=52Hz,问：在幅频谱中，能否分辨f1和f2对应的频率分量？。为什么？。再改变f2=600Hz，在幅频谱中，f2对应的频率分量出现在Hz；问：在fs=1000Hz的情况下，能否正确检测f2对应的频率分量？。为什么？。为了正确检测f2对应的频率分量，则fs至少取多少Hz?Hz。在该程序中改变fs，验证你的结论。iv)比较幅频谱和功率谱，可以发现功率谱具有的特性。FFT实现任意两个序列的快速卷积。%程序fftjuanji.mclearclcx1=input('x1=');x2=input('x2=');%N1=length(x1);N2=length(x2);%序列x1(n),x2(n)的长度E=ceil(log2(N1+N2-1));%ceil向+∞方向取整N=2^E;%x1=[x1,zeros(1,N-N1)];%x2=[x2,zeros(1,N-N2)];X1=fft(x1,N);%X2=fft(x2,N);Y=X1.*X2;%y=ifft(Y,N)%结果分析：1)回到MATLAB窗口，键入：x1=[111],x2=[12],回车。结果:y=2)问：可用Matlab中的什么函数验算上述卷积结果？4.利用谱分析观察太阳黑子周期性。以100年中记录到的太阳黑子出现次数为信号x(n),对x(n)作功率谱，从中观察太阳黑子周期性。%程序taiyangheizi.mclearclcx=[101826635317209215412585683823102483...132131118906760474121166471434454348...4228108201512143546413024167428...17365062677148288135712213810386633724...11154062981249666645439217423559496...77594447301673774];%100年中太阳黑子出现的次数subplot(211);plot(x)%画x(n)N=128;fs=1;%fs=1Hz,N=128点s=x-mean(x);%对x作零均值化处理（去除直流分量）Y=;%对s做N点fftPSD=;%做功率谱PSDf=;%将频率定标为实际频率fsubplot(212);;%画功率谱（N/2点）1)填写空格中的画图语句并绘出结果图形。2)从s的功率谱观察到，其幅度最高处对应的横坐标f=Hz,则太阳黑子每隔年出现一次最高峰。3)在对s做FFT时，为何可取fs=1Hz,N=128点？实验三IIR数字滤波器设计一、实验目的(1)掌握用双线性变换法设计IIR数字低通和高通滤波器。(2)设计低通滤波器对实际心电图信号进行滤波。(3)设计低通滤波器对含有啸叫噪声的音乐信号进行消噪。(4)设计IIR数字低通和高通滤波器对某个DTMF(双音多频)信号进行频带分离。二、实验环境1.Windows98以上操作系统2.安装MATLAB6.0以上版本三、实验原理1．选频型数字滤波器的种类有、、和滤波器。2.从实现方法上，数字滤波器通常分为和滤波器。3．IIR滤波器的设计目的是根据技术指标，找到；IIR滤波的MATLAB语句为y=;4．双线性变换法设计低通IIR滤波器的步骤是(1)(2)(3)四、实验内容1.人体心电图信号在测量过程中往往受到工业高频干扰，必须经过低通滤波处理后才能作为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。试以x(n)作为输入序列，滤除其中的干扰成分。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πrad,ωs=0.3πrad,ap=1dB,as=15dB。已设计出H(z)(p300)其中A=0.09036；B1=1.2686,C1=-0.7051；B2=1.0106,C2=-0.3583；B3=0.9044,C3=-0.2155*IIR滤波的Matlab语句：y=filter(b,a,x)b,aHk(z)分子/分母系数；x输入信号x(n)；y滤波结果y(n)。**求频响特性的Matlab语句：[H,w]=freqz(b,a,N)b,aHk(z)分子/分母系数；N频率点数；H—滤波器的频响特性H(w);w数字频率轴w。(1)用IIR低通滤波器(原理设计)对实际的心电信号进行滤波(IIRECG.m)%程序：(IIRECG.m)clearclcx=[-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];%心电信号x(n)[含高频噪声]A=0.09036;%IIR低通滤波器系数B1=1.2686;C1=-0.7051;B2=1.0106;C2=-0.3583;B3=0.9044;C3=-0.2155;b=[A2*AA];a1=[1-B1-C1];N=128;[H1,w]=freqz(b,a1,N);%频响特性H1(w)y1=filter(b,a1,x);%y1(n)是x(n)经H1(z)滤波的结果a2=[1-B2-C2];[H2,w]=freqz(b,a2,N);%频响特性H2(w)y2=filter(b,a2,y1);%y2(n)是y1(n)经H2(z)滤波的结果a3=[1-B3-C3];[H3,w]=freqz(b,a3,N);%频响特性H3(w)H=H1.*H2.*H3;%总的频响特性H(w)=H1(w)H2(w)H3(w)mag=abs(H);%滤波器的幅频特性db=20*log10((mag+eps)/max(mag));%幅频特性(dB)y3=filter(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);gridon;title('x(n)');subplot(222);plot(wx/pi,abs(X(1:N/2)));gridon;title('|X(wx)|');subplot(223);plot(w/pi,db);gridon;title('|H(w)|(db)');subplot(224);plot(y3);gridon;title('y3(n)');要求：绘出并分析结果图形。(2)用现有的Matlab函数设计上述IIR低通滤波器,对实际的心电信号进行滤波，同样绘出题1的结果图形。Matlab函数:*根据技术指标ωp,ωs,αp,αs计算巴特沃思滤波器的阶数N和3dB截止频率ωc语句：[N,wc]=buttord(wp/pi,ws/pi,ap,as)注：数字频率以π为单位。**根据N,ωc确定数字滤波器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,-2,-2,0,...0,-2,-2,-2,-2,0];%心电信号x(n)[含高频噪声]wp=0.2*pi;ws=0.3*pi;ap=1;as=15;%低通滤波器技术指标N=128;[N1,wc]=buttord(wp/pi,ws/pi,ap,as)%[b,a]=butter(N1,wc)%[H,w]=freqz(b,a,N);%y3=filter(b,a,x);%mag=abs(H);%db=20*log10((mag+eps)/max(mag));%幅频特性(dB)X=fft(x,N);%wx=2*pi*(0:N/2-1)/N;%%%%绘图%%%subplot(221);plot(x);gridon;title('x(n)');subplot(222);plot(wx/pi,abs(X(1:N/2)));gridon;title('|X(wx)|');subplot(223);plot(w/pi,db);gridon;title('|H(w)|(db)');subplot(224);plot(y3);gridon;title('y3(n)');1)在空格中填写程序注释。2)运行上述程序，得：滤波器阶数N=；3dB截止频率wc=；b=；a=。3)修改as=50,得：N=；wc=；b=；a=。2.文件”yinyue.wav”是含有啸叫噪声的音乐信号，首先对其进行谱分析，观察信号和噪声的频带范围，再设计适当的IIR数字低通滤波器对其进行消噪处理恢复原信号，将结果存入音乐文件”yinyuexiaozao.wav”并用耳机监听消噪前后的音乐信号。%程序:IIRxiaozaoy.mclearclcN1=81920;%做fft的点数y1=wavread('yinyue.wav');%读入含噪音乐信号sound(y1)%播放y1pause(10)%暂停10sy1=y1*6;Y1=fft(y1,N1);%PSD1=Y1.*conj(Y1)/N1;%fs=8192;f=8192/N1*(0:N1/2-1);%%%%绘图%%%figure(1)subplot(211);plot(real(y1));gridon;subplot(212);plot(f,PSD1(1:N1/2));axis([050000500]);gridon;实验要求：(1)运行上述程序，绘出结果波形并监听含噪音乐信号。(2)从第二张子图可见音乐信号的频带范围在Hz,啸叫噪声频率为Hz。(3)根据%后的要求编写程序，将这些程序连在上述程序之后。%%IIR低通滤波器%%fp=;%通带截止频率fp=1000Hzfs1=;%阻带截止频率fs1=1200Hzwp=2*pi*fp/fs;%化为数字频率wp和wsws=;Rp=;As=;%Rp=1dB,As=50dB[N,wn]=;%根据技术指标计算巴特沃思滤波器的阶数N和3dB截止频率ωc[b,a]=;%根据N,ωc确定数字滤波器H(z)分子/分母多项式的系数b,a%%%%%%%%%%%%%%%%[H,w]=freqz(b,a,fs,'whole');%频响特性H(w)mag=;%幅频特性db=;%幅频特性(以dB为单位)y=;%用所设计的滤波器对y1进行低通滤波消噪,得y%%%%绘图%%%%figure(2)subplot(211);plot(w*fs/(2*pi),db);axis([0fs/2-400100]);gridon;subplot(212);plot(real(y));gridon;sound(real(y))%播放ywavwrite(real(y),'yinyuexiaozao.wav');%将结果存入音乐文件："yinyuexiaozao.wav"(4)绘出结果波形并监听已消噪的音乐信号。3.DTMF(双音多频)信号是按键电话拨号音，其中”2”键是由下列低频音和高频音所构成：x(n)=sin(2πf1/fs×n)+sin(2πf2/fs×n),f1=697Hz,f2(n=0~1599)编写程序，完成以下功能：(1)产生x(n)并绘图，用”sound”语句播放x(n),用耳机监听,并将x(n)存入音乐文件”DTMF2.wav”。(2)设计适当的低通和高通滤波器对x(n)进行滤波，将低频音和高频音分离,得分离后的信号:x1(n)=sin(2πf1/fs×n),x2(n)=sin(2πf2/fs×n),画出x1(n)和x2(n)，并分别将x1(n)和x2(n)存入音乐文件”DTMFdi.wav”和”DTMFgao.wav”，用耳机监听这两种声音。%程序：DTMF.m波形：实验四FIR数字滤波器设计一、实验目的1.掌握用窗函数法设计FIR低通数字滤波器。2.了解各种窗函数对滤波特性的影响。3.设计FIR低通滤波器对含噪音乐信号进行消噪。4.用二维低通FIR滤波器连接图像中文字的断裂部分。二、实验环境1.Windows98以上操作系统2.安装MATLAB6.0以上版本三、实验原理1．FIR滤波器的设计目的是根据技术指标，设计。FIR滤波的MATLAB语句为y=。2.窗函数(1)矩形窗的表达式为：(2)汉宁窗的表达式为：(3)哈明窗的表达式为：四、实验内容1.人体心电图信号在测量过程中往往受到工业高频干扰，必须经过低通滤波处理后才能作为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。试以x(n)作为输入序列，用窗函数法设计一个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πrad,ωs=0.3πrad,ap=1dB,as=50dB。选哈明窗，即w(n)=[0.54-0.46cos(2πn/(N-1))]RN(n)。*FIR滤波的Matlab语句：y=conv(h,x)h滤波器单位脉冲响应h(n)；x输入信号x(n)；y滤波结果y(n)。程序：%用FIR低通滤波器对实际的心电信号进行滤波(FIRECG.m)clearclcwp=0.2*pi;ws=0.3*pi;%tr=ws-wp;%N=ceil(8*pi/tr)+1%n=[0:1:N-1];wc=(ws+wp)/2;%m=n-(N-1)/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]=freqz(h,[1],1000,'whole');%mag=abs(H);db=20*log10((mag+eps)/max(mag));%delta_w=2*pi/1000;ap=-(min(db(1:1:wp/delta_w+1)))%技术指标验证ap,asAs=-round(max(db(ws/delta_w+1:1:501)))%%%绘图%%%figure(1)subplot(221);stem(n,hd,'.');gridon;title('hd(n)');subplot(222);stem(n,w_ham,'.');gridon;title('哈明窗w_ham(n)');subplot(223);stem(n,h,'.');gridon;title('h(n)');subplot(224);plot(w(1:501)/pi,db(1:501));gridon;title('H(w)');%%%%%%%%%%x=[-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];%含噪的ECG信号y=conv(x,h);%X=fft(x,128);%Y=fft(y,128);%wx=2*pi*(0:N/2-1)/N;%%%%绘图%%%figure(2)subplot(221);plot(x);axis([060-10020]);gridon;title('x(n)');subplot(222);stem(wx/pi,abs(X(1:N/2)),'.');gridon;title('|X(w)|');subplot(223);plot(y);axis([(N-1)/2(N-1)/2+60-10020]);gridon;title('y(n)');subplot(224);stem(wx/pi,abs(Y(1:N/2)),'.');gridon;title('|Y(w)|');%%%绘图%%%要求：绘出结果图形，并观察图2的第1和第3张小图，可见，y(n)和x(n)相比，延迟了点。运行上述程序，得：滤波器阶数N=；实际的ap=；实际的as=；可见，在设计相同技术指标的滤波器时，FIR滤波器的阶数IIR滤波器的阶数。修改窗函数为汉宁窗，得：滤波器阶数N=；实际的ap=；实际的as=；2.文件”noisy.wav”是含有高频噪声的音乐信号，首先对其进行谱分析，观察信号和噪声的频带范围，再设计适当的FIR低通滤波器对其进行消噪处理恢复原信号，将结果存入音乐文件”FIRxiaozao.wav”并用耳机监听消噪前后的音乐信号。%程序:yinyueFIR.mclearclcfs=16000;x=wavread('noisy.wav');%读入声音含噪音乐文件到xNx=length(x);n=0:Nx-1;Ts=1/fs;t=n*Ts;L=ceil(log2(Nx));N=2^L;X=fft(x,N);%对x做频谱Xf=fs/N*(0:N/2-1);figure(1)subplot(211);plot(t,x);%画xsubplot(212);plot(f,abs(X(1:N/2)));%画|X(f)|(N/2点）实验要求：(1)运行上述程序，绘出结果波形并监听含噪音乐信号。(2)从第二张子图可见音乐信号的频带范围在Hz,噪声频带在Hz。(3)根据%后的要求编写程序，将这些程序连在上述程序之后。%%%%%%FIR低通滤波器%%%%fp=;%通带截止频率fp=1500Hzfs1=;%阻带截止频率为fs1=1800Hzwp=2*pi*fp/fs;%技术指标wp和ws1ws1=2*pi*fs1/fs;tr=;%过渡带tr=ws1-wpNFIR=;%滤波器的阶数NFIRn=[0:1:NFIR-1];%滤波器的时间范围nwc=;%理想低通滤波器的截止频率m=;%m=n-(NFIR-1)/2+eps;hd=;%理想低通滤波器hd(n)w_ham=;%哈明窗h=;%加哈明窗之后的实际低通滤波器h(n)[H,w]=freqz(h,[1],N,'whole');%低通滤波器的频响特性mag=;%幅频特性magdb=20*log10((mag+eps)/max(mag));%幅频特性(dB)subplot(221);;gridon;title('hd(n)');%画hd(n)subplot(222);;%画哈明窗w_ham(n)gridon;title('哈明窗w_ham(n)');subplot(223);;%画h(n)gridon;title('h(n)');subplot(224);plot(w(1:N/2)/pi,db(1:N/2));%画dbgridon;title('H(w)');%%%%%%%%%%%%%%低通滤波%%%%%%%%%%%%%%%%y=;%低通滤波消噪y(n)=x(n)*h(n)Y=fft(y,N);ny=0:length(y)-1;ty=ny*Ts;figure(3)subplot(211);plot(ty,y);%画ysubplot(212);plot(f,abs(Y(1:N/2)));%画|Y(f)|(N/2点)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sound(x,fs);%以fs=16000Hz播放xpause(5)%停顿5ssound(y,fs);%以fs=16000Hz播放ywavwrite(y,'FIRxiaozao.wav');%将y存入声音文件"FIRxiaozao.wav"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(4)绘出结果波形并监听已消除高频噪声的音乐信号。4.二维低通滤波器可进行图像平滑。用卷积核对图像f(m,n)进行低通滤波，连接图像中文字的断裂部分。%程序：wenzi.mclearclcf=imread(‘ea.jpg’);%读入文字图像’ea.jpg’到二维数组ff=double(f);%将f转换为双精度型N=11;%卷积核的尺寸Nh=1/(N*N)*ones(N,N);%卷积核hy=imfilter(f,h);%用h对f进行二维低通滤波(平滑),得输出图像yH=fft2(h,32,32);%对卷积核h做二维fftsubplot(211);imshow(f,[]);%显示fsubplot(223);mesh(abs(fftshift(H)));%显示H的幅频谱subplot(224);imshow(y,[]);%显示y实验要求：运行上述程序，绘出结果波形。改变N=5或N=11,问：输出图像有何改变?实验五离散系统的时域及复频域分析一、实验目的1.掌握Matlab编程求解离散时间系统差分方程的两种方法：迭代法和filter函数法。2.掌握利用Z变换对系统进行复频域分析。3.掌握系统零极点的绘制方式及Z反变换的求解方法。4.熟悉Z变换的应用。二、实验环境1.Windows2000以上操作系统2.安装MATLAB6.0以上版本三、实验原理系统的输出y(n)与输入x(n)及单位脉冲响应h(n)的关系:一个N阶常系数线性差分方程表达式:3.系统函数H(z)与X(z),Y(z)的关系为：H(z)=,N阶常系数线性差分方程对应的系统函数H(z)=。4.H(z)的零点定义为：；H(z)的极点定义为：。5.本实验涉及的Matlab函数1)求解差分方程：y=filter(b,a,x)其中：x输入信号；y差分方程的解。b,a输入/输出信号前的系数。2)画零极点图:zplane(b,a,n)其中：b,aH(z)分子/分母系数。3)求z反变换[r,p,k]=residuez(b,a)bH(z)的分子多项式系数；aH(z)的分母多项式系数；输出：r,p,k的含义如下：“residuez”可将H(z)分解为简单有理分式之和：式<1>中p(1),p(2)…p(n)的集合为列向量p;r(1),r(2)…r(n)的集合为列向量r;k(1),k(2)…的集合为行向量k;已知r,p,k即可手工写出h(n)的表达式：h(n)={r(1)[p(1)]n+r(2)[p(2)]n+…r(n)[p(n)]n}u(n)+k(1)δ(n)+k(2)δ(n-1)+……四、实验内容1.已知：系统的差分方程为：y(n)-y(n-1)+0.9y(n-2)=x(n)1)分别利用迭代法和filter函数法求此系统的单位脉冲响应h(n),并绘图。[1]迭代法：%程序：shiyan511.mclearclcN=100;x=[1zeros(1,N-1)];%产生x(n)=δ(n)y(1)=0;%y(-2)=0y(2)=0;%y(-1)=0y(3)=1;%y(0)=1forn=4:N%迭代法求解该差分方程,得y(n)y(n)=x(n)+y(n-1)-0.9*y(n-2);endk=-2:N-3;%时间轴范围kstem(k,y,'.');%画y(n)=h(n)title('系统的单位脉冲响应h(n)');%写标题h(n)波形：[2]filter函数法：%程序：shiyan512.mclearclcN=100;x=[1zeros(1,N-1)];%产生x(n)=δ(n)b=[1];%差分方程系数b,aa=[1,-1,0.9];h=filter(b,a,x);%用filter函数求解系统的单位脉冲响应h(n)k=-2:N-3;stem(k,h,'.');title('系统的单位脉冲响应h(n)');h(n)波形：2)修改输入信号为x(n)=u(n),求系统的阶跃响应y1(n)。[u(n)由：u=ones(1,n)产生]。%程序：shiyan513.my1(n)波形：3)利用由1)产生的h(n),利用卷积法产生系统的阶跃响应：y2(n)=u(n)*h(n)。%程序：shiyan514.my2(n)波形：问：y1(n)和y2(n)是否相同？2.用求解差分方程的方式（迭代法及filter函数法）实现一个数字振荡器（正弦波信号）：h(n)=sin(ωkn),其中ωk=2πfk/fs,频率fk=1000Hz,取样频率fs=10kHz。[提示：h(n)对应的差分方程为：h(n)-ah(n-1)+h(n-2)=bδ(n-1)<1>a=2cos(ωk)=1.618,b=sin(ωk)=0.588]画出h(n)并用“zplane”画该系统的零极点图。%迭代法实现程序shiyan521.mh(n)波形：%filter函数法实现程序shiyan522.mh(n)波形：问：极点位置处于单位圆（内/上/外）？3.求下列序列的z变换并用”zplane”画出零极点分布图。1)x(n)={1,2,1,3}该序列的z变换为：X(z)=1+2z-1+z-2+3z-3%程序shiyan531.mclearclcb=[1