《数字信号处理》课程研究性学习报告

《数字信号处理》课程研究性学习报告指导教师薛健时间2014.6【目的】(1)掌握IIR和FIR数字滤波器的设计和应用；(2)掌握多速率信号处理中的基本概念和方法；(3)学会用Matlab计算小波分解和重建。(4)了解小波压缩和去噪的基本原理和方法。【研讨题目】一、(1)播放音频信号yourn.wav，确定信号的抽样频率，计算信号的频谱，确定噪声信号的频率范围；(2)设计IIR数字滤波器，滤除音频信号中的噪声。通过实验研究sP,，sP,AA的选择对滤波效果及滤波器阶数的影响，给出滤波器指标选择的基本原则，确定你认为最合适的滤波器指标。（3）设计FIR数字滤波器，滤除音频信号中的噪声。与（2）中的IIR数字滤波器，从滤波效果、幅度响应、相位响应、滤波器阶数等方面进行比较。【设计步骤】【仿真结果】【结果分析】由频谱知噪声频率大于3800Hz。FIR和IIR都可以实现滤波，但从听觉上讲，人对于听觉不如对图像（视觉）明感，没必要要求线性相位，因此，综合来看选IIR滤波器好一点，因为在同等要求下，IIR滤波器阶数可以做的很低而FIR滤波器阶数太高，自身线性相位的良好特性在此处用处不大。【自主学习内容】MATLAB滤波器设计【阅读文献】老师课件，教材【发现问题】(专题研讨或相关知识点学习中发现的问题)：过渡带的宽度会影响滤波器阶数N【问题探究】通过实验，但过渡带越宽时，N越小，滤波器阶数越低，过渡带越窄反之。这与理论相符合。【仿真程序】信号初步处理部分：[x1,Fs,bits]=wavread('yourn.wav');sound(x1,Fs);y1=fft(x1,1024);f=Fs*(0:511)/1024;figure(1)plot(x1)title('原始语音信号时域图谱');xlabel('timen');ylabel('magnituden');figure(2)freqz(x1)title('频率响应图')figure(3)subplot(2,1,1);plot(abs(y1(1:512)))title('原始语音信号FFT频谱')subplot(2,1,2);plot(f,abs(y1(1:512)));title(‘原始语音信号频谱')xlabel('Hz');ylabel('magnitude');IIR：fp=2500;fs=3500;wp=2*pi*fp/FS;ws=2*pi*fs/FS;Rp=1;Rs=15;Ts=1/Fs;wp=2*pi*fp/FS;ws=2*pi*fs/FS;wp1=2/Ts*tan(wp/2);ws1=2/Ts*tan(ws/2);t=0:1/11000:(size(x1)-1)/11000;Au=0.03;d=[Au*cos(2*pi*5000*t)]';x2=x1+d;[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');[Z,P,K]=buttap(N);[Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);[bz,az]=bilinear(b,a,Fs);%[H,W]=freqz(bz,az);figure(4)plot(W*Fs/(2*pi),abs(H))gridxlabel('频率/Hz')ylabel('频率响应幅度')title('Butterworth')f1=filter(bz,az,x2);figure(5)subplot(2,1,1)plot(t,x2)title('滤波前时域波形');subplot(2,1,2)plot(t,f1);title('滤波后时域波形');sound(f1,FS);FIR[x1,Fs,bits]=wavread('I:/dsp_2014_project3/yourn');fp=2500;fs=3500;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Rs=50;M=ceil((Rs-7.95)/(ws-wp)/2.285);M=M+mod(M,2);beta=0.1102*(Rs-8.7);w=kaiser(M+1,beta);wc=(wp+ws)/2;alpha=M/2;k=0:M;hd=(wc/pi)*sinc((wc/pi)*(k-alpha));h=hd.*w';f1=filter(h,[1],x1);[mag,W]=freqz(h,[1]);figure(1)plot(W*Fs/(2*pi),abs(mag));grid;xlabel('频率/Hz');ylabel('频率响应幅度');title('Kaiser´窗设计Ⅰ型线性相位FIR低通滤波器');figure(2)subplot(2,1,1)plot(t,x1)title('滤波前时域波形');subplot(2,1,2)plot(t,f1);title('滤波后时域波形');sound(f1,Fs);二、(1)音频信号kdqg24k.wav抽样频率为24kHz，用y=wavread('kdqg24k');sound(y,16000);播放该信号。试用频域的方法解释实验中遇到的现象;(2)设计一数字系统，使得sound(y,16000)可播放出正常的音频信号；讨论滤波器的频率指标、滤波器的的类型(IIR,FIR)对系统的影响。【仿真结果】【结果分析】24K的信号用16K播放，频谱会被拉宽，无法正常播放，通过2倍内插，通过滤波器，然后3倍抽取，得到的信号用16K播放器就能正常播放。【自主学习内容】功能：对时间序列进行重采样。格式：1.y=resample(x,p,q)采用多相滤波器对时间序列进行重采样，得到的序列y的长度为原来的序列x的长度的p/q倍，p和q都为正整数。此时，默认地采用使用FIR方法设计的抗混叠的低通滤波器。2.y=resample(x,p,q,n)采用chebyshevIIR型低通滤波器对时间序列进行重采样，滤波器的长度与n成比例，n缺省值为10.3.y=resample(x,p,q,n,beta)beta为设置低通滤波器时使用Kaiser窗的参数，缺省值为5.4.y=resample(x,p,q,b)b为重采样过程中滤波器的系数向量。5.[y,b]=resample(x,p,q)输出参数b为所使用的滤波器的系数向量。说明：x--时间序列p、q--正整数，指定重采样的长度的倍数。n--指定所采用的chebyshevIIR型低通滤波器的阶数，滤波器的长度与n成比列。beta--设计低通滤波器时使用Kaiser窗的参数，缺省值为5.【阅读文献】PPt课本【发现问题】(专题研讨或相关知识点学习中发现的问题)：采样频率与播放频率之间不是整数倍关系【问题探究】此时内插和抽取结合实现正常播放【仿真程序】fs=24000;x1=wavread('I:/dsp_2014_project3/kdqg24k.wav');sound(x1,16000);y1=fft(x1,1024);f1=fs*(0:511)/1024;f2=fs/2*3*(0:511)/1024;figure(1)subplot(2,1,1);plot(f1,abs(y1(1:512)));title('原始语音信号24K正常播放频谱')xlabel('Hz');ylabel('magnitude');subplot(2,1,2);plot(f2,abs(y1(1:512)));title('原始语音信号16K播放频谱')xlabel('Hz');ylabel('magnitude');y=resample(x1,2,3);sound(y,16000);y2=fft(y,1024);figure(2)subplot(2,1,1);plot(f1,abs(y1(1:512)));title('原始语音信号24K正常播放频谱')xlabel('Hz');ylabel('magnitude');subplot(2,1,2);plot(f2,abs(y2(1:512)));title('原始语音信号16K经过设计的数字系统后播放频谱')xlabel('Hz');ylabel('magnitude');三、对连续信号x(t)=40t2(1t)4cos(12t)[0t1]+40(t1)4(2t)2cos(40t)[1t2]+0.1n(t)在区间[0，2]均匀抽样1024点得离散信号x[k]，其中n(t)是零均值方差为1的高斯噪声(1)画出信号x(t)的波形；(2)计算并画出db7小波的5级小波变换系数；(3)通过观察小波系数，确定阈值化处理的阈值；(4)对小波系数进行阈值化处理，画出去噪后的信号波形，求出最大误差和均方误差；(5)对近似系数和小波系数均进行阈值化处理，画出去噪后的信号波形，求出最大误差和均方误差；(6)用Haar小波基，重复(3)-(5)；(7)讨论所得结果。【仿真结果】（1）050010001500200025003000350040004500-1.5-1-0.500.511.5signalwithnoise（2）00.20.40.60.811.21.41.61.82-6-4-202468Waveletcoefficients（3）T=4.3000t=96.84%(4)db7小波基050010001500200025003000350040004500-1-0.500.51Ronconstructedsignalwithnoniose050010001500200025003000350040004500-1-0.500.51ErrorEmabs=0.9166(5)Haar小波基T=3.100;t=96.18%00.20.40.60.811.21.41.61.82-8-6-4-202468Waveletcoefficients050010001500200025003000350040004500-1-0.500.51Ronconstructedsignalwithnoniose050010001500200025003000350040004500-2-1012ErrorEmabs=1.1341(6)db14小波基T=4.6，t=94.39%00.20.40.60.811.21.41.61.82-8-6-4-202468Waveletcoefficients050010001500200025003000350040004500-1-0.500.51Ronconstructedsignalwithnoniose050010001500200025003000350040004500-1-0.500.51ErrorEmabs=0.8965【结果分析】选用dbp系列小波基对信号去噪时，p值越大，去噪的效果越好，最大重建误差也小【自主学习内容】MATLAB小波变换【阅读文献】PPT和课本【发现问题】(专题研讨或相关知识点学习中发现的问题)：选用dbp系列小波基对信号去噪时，p值越大，去噪的效果越好，最大重建误差也小【问题探究】不同dp小波基不同，会对结果产生不同影响【仿真程序】N=4096;k=linspace(0,2,N);nt=randn(size(k));x=40*k.^2.*(1-k).^4.*cos(12*pi*k).*(0k&k1)+40.*(k-1).^4.*(2-k).^2.*cos(80*pi*k).*(1k&k2)+0.1*nt;figure;plot(x);title('signalwithnoise');dwtmode('per');[C,L]=wavedec(x,6,'db14')figure;plot(k,C);title('Waveletcoefficients');M=0;fork=1:4096;T=4.6;ifabs(C(1,k))=T;C(1,k)=0;endifC(1,k)~=0;M=M+1;endendA1=C.*C;U1=0;A2=x.*x;U2=0;