数字信号实验

实验步骤1、编写nonrec.m函数文件，实现FIR滤波。这里给定，，求。2、编写rec.m函数文件，实现纯递归IIR滤波。这里给定，，，，求单位取样响应。3、用nonrec.m和rec.m函数编制dfilter.m函数文件，构成完整的一般IIR滤波器程序，并完成下列信号滤波这里给定系统函数，计算。4、用help查看内部函数filter.m，了解调用格式，重做3，并和你编写的dfilter.m的结果进行比较。实验内容1、编程实现FIR滤波2、编写时实现纯递归IIR滤波；差分方程：系统函数：3、调用库函数filter.m实现IIR滤波实验数据1.nonrec.m函数文件functiony=nonrec(h,x)m=length(h);n=length(x);fork=1:(m+n-1)S=0;ifk=nn1=k;form1=1:mifn1=1S=S+h(m1)*x(n1);y(k)=S;n1=n1-1;endendelsen2=n;form1=(1+k-n):mifn2=1S=S+h(m1)*x(n2);y(k)=S;n2=n2-1;endendendendendh=ones(1,8);x=0:15;y=nonrec(h,x)y=0136101521283644526068768492847565544229152.rec.m函数文件functiony=rec(x,a)La=length(a);Lx=length(x);y(1)=x(1);fork=2:100A=0;ifk=Lxifk=La+1form=1:k-1A=A+a(k-m)*y(m);y(k)=A+x(k);endelseforc=1:LaA=A+a(c)*y(k-c);y(k)=x(k)+A;endendelselwfree.cnelseford=1:LaA=A+a(d)*y(k-d);y(k)=A;endendendendend3.dfilter.m函数文件functiony=dfilter(x,a,b)y1=rec(x,a);y=nonrec(b,y1);a1=2*0.95*cos(pi/8);a2=-0.95*0.95;a=[a1a2];x=[1];y=rec(x,a);a=[0.6,-0.5]b=[1-31];n=0:63;x=sin(3*pi*n/5);y1=dfilter(x,a,b);y2=filter(b,a,x);subplot(2,1,1),stem(y1),gridon,title('dfilter.m计算结果')subplot(2,1,2),stem(y2),gridon,title('filter.m计算结果')内部函数filter.m的调用格式是y=filter(b,a,X)编写的dfilter.m调用格式是y=dfilter(x,a,b)实验总结通过本此上机实践熟悉和掌握了数字信号出来的因果性数字系统的时域实现，进一步了解和认识了数字信号的MATLAB实现及用MATLAB实现数字信号处理。176本文来自网学()，转载请注明出处：实验一离散傅里叶变换的性质及应用一、实验目的1、了解DFT的性质及应用。2、熟悉MATLAB编程的特点。二、实验内容1、用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换X(ejw),并比较三种程序计算机运行时间。（1）用forloop语句的M函数文件dft1.m，用循环变量逐点计算X(k)；（2）编写用MATLAB矩阵运算的M函数文件dft2.m,完成下列矩阵运算；（3）调用FFT库函数，直接计算X(k)；（4）分别利用上述三种不同方式编写的DFT程序计算序列x(n)的傅立叶变换X(ejw)，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。M函数文件如下：dft1.m:function[Am,pha]=dft1(x)N=length(x);w=exp(-j*2*pi/N);fork=1:Nsum=0;forn=1:Nsum=sum+x(n)*w^((k-1)*(n-1));endAm(k)=abs(sum);pha(k)=angle(sum);enddft2.m:function[Am,pha]=dft2(x)N=length(x);n=[0:N-1];k=[0:N-1];w=exp(-j*2*pi/N);nk=n'*k;wnk=w.^(nk);Xk=x*wnk;Am=abs(Xk);pha=angle(Xk);dft3.m:function[Am,pha]=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk);源程序及运行结果：（1）x=[ones(1,8),zeros(1,248)];t=cputime;[Am1,pha1]=dft1(x);t1=cputime-tn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(1)subplot(2,1,1),plot(w,Am1,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha1,'r');grid;title('PhasePart');xlabel('frequencyinradians');ylabel('argX[exp(jw)]/radians');运行时间：t1=2.25s（2）x=[ones(1,8),zeros(1,248)];t=cputime;[Am2,pha2]=dft2(x);t2=cputime-tn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(2)subplot(2,1,1),plot(w,Am2,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha2,'r');grid;title('PhasePart');xlabel('frequencyinradians');ylabel('argX[exp(jw)]/radians');运行时间：t2=0.77s（3）x=[ones(1,8),zeros(1,248)];t=cputime;[Am3,pha3]=dft3(x);t3=cputime-t;n=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(3)subplot(2,1,1),plot(w,Am3,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha3,'r');grid;title('PhasePart');xlabel('frequencyinradians');ylabel('argX[exp(jw)]/radians');运行时间：t3=0从以上运行结果可以看出，调用FFT库函数直接计算X(k)速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。因此FFT算法大大提高了DFT的实用性。2、利用DFT实现两序列的卷积运算，并研究DFT点数与混叠的关系。给定x(n)=nR16(n)，h(n)=R8(n)，用FFT和IFFT分别求线性卷积和混叠结果输出，并用函数stem(n,y)画出相应图形。源程序如下：N=16;x=[0:N-1];h=ones(1,8);线性卷积：Xk1=fft(x,23);%做23点fftHk1=fft(h,23);Yk1=Xk1.*Hk1;y1=ifft(Yk1);n=0:22;figure(1)stem(n,y1)混叠：Xk2=fft(x);Hk2=fft(h,16);%做16点fftYk2=Xk2.*Hk2;y2=ifft(Yk2);n=0:15;figure(2)stem(n,y2)运行结果：线性卷积：混叠结果：将两个信号的DFT相乘，再做反变换，相当于时域做循环卷积。而循环卷积是线性卷积以N为周期周期性延拓的结果，其中N为离散傅立叶变换的点数。若x1(n)是N1点序列，x2(n)是N2点序列，则其线性卷积为N1+N2-1点，做N=max(N1,N2)点的傅立叶变换，则结果中将有N2-1点是混叠的。题目中N1=16,N2=8,由程序运行结果可以看出，第二张图的前7（N2-1）个点为混叠结果，后9个点是线性卷积的结果。所以，要用DFT来做线性卷积，DFT的点数必须大于或等于线性卷积的结果。本题中取N1+N2-1=23为DFT的点数，即可得到正确的线性卷积结果。3、研究高密度频谱与高分辨率频谱设有连续信号Xa(t)=cos(2π×6.5×10^3t)+cos(2π×7×10^3t)+cos(2π×9×10^3t)以采样频率fs=32kHz对信号Xa(t)采样,分析下列三种情况的幅频特性。（1）采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。（2）采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。（3）采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。源程序如下：fs=32000;%采样频率N=16;%采集16点n=0:N-1;%做16点DFTxa=cos(2*pi*6.5*10^3*n/fs)+cos(2*pi*7*10^3*n/fs)+cos(2*pi*9*10^3*n/fs);[Am1,pha1]=dft3(xa);n=[0:(length(xa)-1)];w=(2*pi/length(xa))*n;figure(1)plot(w,Am1,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');x=[xa(1:16),zeros(1,240)];%补零到256[Am2,pha2]=dft2(x);%做256点DFTn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(2)plot(w,Am2,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');N0=256;%采集256点n0=0:N0-1;%做256点DFTxa0=cos(2*pi*6.5*10^3*n0/fs)+cos(2*pi*7*10^3*n0/fs)+cos(2*pi*9*10^3*n0/fs);[Am3,pha3]=dft3(xa0);n=[0:(length(xa0)-1)];w=(2*pi/length(xa0))*n;figure(3)plot(w,Am3,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');运行结果：N=16点DFT幅频特性高密度频谱：高分辨率频谱：观察运行结果可知，N=16点时的DFT，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。将16点DFT补零到256点，做256点DFT,由结果可以看出，频率间隔缩小（由2?/16减小为2???256），得到了比较平滑的连续曲线，此为高密度频谱。它是由16点经