数字信号处理实验报告

实验报告实验一：1．用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换;比较三种程序计算机运行时间。(1)编制用forloop语句的M函数文件dft1.m，用循环变量逐点计算X(k)；(2)编写用Matlab矩阵运算的M函数文件dft2.m，完成下列矩阵运算；(3)调用FFT库函数，直接计算X(k)；(4)编写M程序文件，分别得到序列x(n)的傅里叶变换，并画出相应的幅频和相频特性；比较用三种不同方式编写的DFT程序各自的计算机运行时间。%%%主程序：x=[ones(1,8),zeros(1,248)];t=cputime;[Am1,pha1]=dft1(x);t1=cputime-t;n=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(1);subplot(2,1,1);plot(w,Am1,'b');title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha1,'r');grid;xlabel('frequencyinradians');ylabel('arg(X[exp(jw)]/radians|');%%%%子程序dft1.mfunction[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);End%%%%子程序dft2.mfunction[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.mfunction[Am,pha]=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk);运行结果：t1=0.1563t2=0.2969t3=02．用ＤＦＴ实现两序列的卷积运算，并研究ＤＦＴ点数与混迭的关系给定：x(n)＝n*R16（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)3．高密度频谱与高分辨率频谱的研究对连续信号)1000*9*2cos()1000*7*2cos()1000*5.6*2cos()(ttttxa以fs=32kHz对xa(t)采样：1)做N=16点的采样和补零到256点的频谱；2)做N=256点的采样，并做出其FFT幅度特性；做出分析。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))|');实验2IIR数字低通Butterworth滤波器的设计一．编写设计数字滤波器所需函数：1．确定滤波的阶数N和3dB截止频率：[N,Wc]=Buttord(Wp,Ws,Rp,Rs)函数returnstheorderNofthelowestorderdigitalButterworthfilterthatlosesnomorethanRpdBinthepassbandandhasatleastRsdBofattenuationinthestopband.WpandWsarethepassbandandstopbandedgefrequencies，andnormalizedfrom0to1(where1correspondstopiradians).Forexample，Lowpass:Wp=0.1，Ws=0.22．确定滤波器的系统函数[B,A]=BUTTER(N,Wc)DesignsanNthorderlowpassdigitalButterworthfilterandreturnsthefiltercoefficientsinlengthN+1vectorsB(numerator)andA(denominator).Thecoefficientsarelistedindescendingpowersofz.3．Freqz：数字滤波器的频率响应[H,w]=freqz(B,A,n)4．Zplane:画零极点图zplane(B,A)5．Filter：利用设计的滤波器对信号进行滤波y=filter(B,A,x)三．实验内容和要求1．设计满足下列指标要求的数字低通滤波器：通带起始频率Wp=0.2π，通带允许起伏αp=-3dB阻带起始频率Ws=0.4π，阻带衰减αs=-30dB2．连续信号s1(t)=sin(2*π*100*t)，以fs=2000Hz对s1(t)采样(采100个点)；此信号在传输过程中受到一噪声信号s2(t)的影响，其幅度约为s1(t)的1/10，s2(t)用randn(1,N)模拟。试观察此信号（s1(t)+s2(t)）通过以上设计的滤波器后的波形。3．要求：①画出设计数字滤波器的幅频特性曲线；②画出s1(t)，s1(t)+s2(t)，以及s1(t)+s2(t)通过滤波器后的波形。脉冲响应不变法：T=1;wp=0.2*pi/T;ws=0.4*pi/T;rp=-3;rs=-30;[N,wc]=buttord(wp,ws,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=impinvar(B,A);[Hz,W]=freqz(Bz,Az);plot(W/pi,20*log10(abs(Hz)));gridon;axis([01-18010]);xlabel('w/π');ylabel('幅度/dB)');fs=2000;N=100;n=0:N-1;s1=sin(2*pi*100*n/fs);s2=0.1*randn(1,N);s=s1+s2;y=filter(Bz,Az,s);figure,subplot(3,1,1),plot(n,s1);subplot(3,1,2),plot(n,s);subplot(3,1,3),plot(n,y);双线性变换法：T=1;Fs=1/T;wpz=0.2;wsz=0.4;wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=-3;rs=-30;[N,wc]=buttord(wp,ws,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs);[Nd,wdc]=buttord(wpz,wsz,rp,rs);[Bdz,Adz]=butter(Nd,wdc);[Hz,W]=freqz(Bdz,Adz);plot(W/pi,20*log10(abs(Hz)));gridon;axis([01-18010]);xlabel('w/π');ylabel('幅度/dB)');fs=2000;N=100;n=0:N-1;s1=sin(2*pi*100*n/fs);s2=0.1*randn(1,N);s=s1+s2;y=filter(Bdz,Adz,s);figure,subplot(3,1,1),plot(n,s1);subplot(3,1,2),plot(n,s);subplot(3,1,3),plot(n,y);实验三、四语音信号的滤波一、实验内容1.录制一段自己的语音信号（单通道，采样率8kHz）；2.画出采样后的语音信号的时域波形和频谱图，回放该语音信号；3.给语音信号加噪声（自己选择噪声），并画出加噪语音信号的时域波形和频谱图；回放加噪后的语音信号。4.IIR滤波：用butterworth滤波器完成滤波，画出滤波器的幅频特性图，butterworth滤波器的指标自己根据实际情况设定，并回放滤波后的语音信号。closeall;clearall;[y,fs,bits]=wavread('t4.wav');%读入声音文件,y是内容,fs是采样频率,bits是比特数sound(y,fs,bits);%回放声音Y1=fft(y);%计算频谱figure,subplot(2,1,1),plot(y),title('原始语音信号');%做原始语音信号的时域图形subplot(2,1,2);plot(abs(Y1)),title('原始语音信号的频谱');%做频域图形yn=awgn(y,20);%加入高斯白噪声，信噪比为20db%n=length(y);%t=0:1/fs:(n-1)/fs;%s=sin(2*2000*pi*t);%s=s';%yn=y+s;%加入正弦噪声sound(yn,fs,bits);%回放声音Yn=fft(yn);figure,subplot(2,1,1),plot(yn),title('加噪声后的语音信号');subplot(2,1,2);plot(abs(Yn)),title('加噪声后的信号频谱');[Bz,Az]=butt();%滤波器y2=filter(Bz,Az,yn);sound(y2,fs,bits);%回放声音Y2=fft(y2);figure,subplot(2,1,1),plot(y2),title('滤波后的语音信号');subplot(2,1,2);plot(abs(Y2)),title('滤波后的信号频谱');function[Bz,Az]=butt()wpz=0.3;wsz=0.5;rp=1;rs=40;%0.3,0.6[N,WN]=buttord(wpz,wsz,rp,rs);[Bz,Az]=butter(N,WN);[Hz,W]=freqz(Bz,Az);fi