数字处理课程设计.doc

x,fs,bits=wavread(hexian.wav); sound(x,fs,bits); N=length(x); fn=3100; t=0:1/fs:(N-1)/fs; x=x;y=x+sin(fn*2*pi*t); sound(y,fs,bits); X=abs(fft(x); Y=abs(fft(y);X=X(1:N/2); Y=Y(1:N/2);deltaf=fs/N; f=0:deltaf:fs/2-deltaf; subplot(2,2,1);plot(t,x);xlabel(时间(单位:s);ylabel(幅度);title(原始和弦信号);gridsubplot(2,2,2);plot(f,X);xlabel(频率(单位:Hz);ylabel(幅度谱);title(原始和弦信号幅度谱图);gridsubplot(2,2,3);plot(t,y);xlabel(时间(单位:s);ylabel(幅度);title(加噪后的和弦信号);gridsubplot(2,2,4);plot(f,Y);xlabel(频率(单位:Hz);ylabel(幅度谱);title(加噪后的和弦信号幅度谱图);grid fpd=3000;fsd=3050;fsu=3150;fpu=3200;Rp=1;As=56;fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu);wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(6.6*pi/dw)+1; n=0:M-1;w_parzen=(parzenwin(M); hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); h_bs=w_parzen.*hd_bs; db,mag,pha,grd,w=freqz_m(h_bs,1); Rp=-maxAs=-subplot(2,2,1);plot(w,db);title(滤波器幅度响应图);xlabel(w/pi);ylabel(db);grid on;axis(0.5,1,-100,50)subplot(2,2,2);plot(w,mag);title(滤波器幅度响应图);xlabel(w/pi);ylabel(幅度mag);grid on;axis(0,1.5,-0.5,2)subplot(2,2,3);plot(w,pha);title(滤波器相位响应图);xlabel(w/pi);ylabel(相位pha);grid on;axis(0,2,-4,4)subplot(2,2,4);plot(n,h_bs);title(滤波器脉冲响应图);xlabel(n);ylabel(h(n);grid on; axis(0,1500,0,1.5)subplot(3,2,1);plot(t,x);title(原始信号时间);xlabel(时间t);ylabel(幅度);grid on;axis(0,40,-1,1)subplot(3,2,2);plot(f,X);title(原始信号幅度谱X);xlabel(频率f);ylabel(幅度);grid on;axis(2500,5000,0,1000)subplot(3,2,3);plot(t,y);title(加噪声后信号时间x1);xlabel(时间t);ylabel(幅度);grid on;axis(0,40,-2,2)subplot(3,2,4);plot(f,Y);title(加噪声后的幅度谱X1);xlabel(频率f);ylabel(幅度);grid on;axis(2500,5000,0,1000) subplot(3,2,5);plot(t,y_fil);title(滤波后信号时间y);xlabel(时间t);ylabel(幅度);grid on;axis(0,40,-1,1)subplot(3,2,6);plot(f,Y_fil);title(滤波后信号幅度谱Y);xlabel(频率f);ylabel(幅度);grid on;axis(2500,5000,0,1000) sound (y_fil,fs,bits);x,fs,bits=wavread(hexian.wav)fs = 22050bits =16 Rp=-(min(db(1:1:wcd/delta_w+1);Rp =4.4096As=-round(max(db(wcu/delta_w+1:1:501);As=52函数定义：function hd=ideal_lp(wc,M);%Ideal Lowpass filter computation%-%hd=ideal_lp(wc,M)% hd=ideal impulse response between 0 to M-1% wc=cutoff frequency in radians% M=length of the ideal filter%alpha=(M-1)/2;n=0:1:(M-1);m=n-alpha+eps;hd=sin(wc*m)./(pi*m);1. freqz_m函数定义：function db,mag,pha,grd,w = freqz_m(b,a);2. % Modified version of freqz subroutine3. % -4. % db,mag,pha,grd,w = freqz_m(b,a);5. %db = Relative magnitude in dB computed over 0 to pi radians6. % mag = absolute magnitude computed over 0 to pi radians 7. % pha = Phase response in radians over 0 to pi radians8. % grd = Group delay over 0 to pi radians9. % w = 501 frequency samples between 0 to pi radians10. % b = numerator polynomial of H(z) (for FIR: b=h)11. % a = denominator polynomial of H(z) (for FIR: a=1)12. %13. H,w = freqz(b,a,1000,whole);14. H = (H(1:1:501); w = (w(1:1:501);15. mag = abs(H);16. db = 20*log10(mag+eps)/max(mag);17. pha = angle(H);18. %pha = unwrap(ang