数字信号处理matlab作业.docx

1、数字信号处理matlab作业M1-1fs=10;ts=1/fs;n=-2:ts:2;f=cos(6*pi*n);subplot(3,1,1),stem(n,f,filled);g=cos(14*pi*n);subplot(3,1,2),stem(n,g,filled);h=cos(26*pi*n);subplot(3,1,3),stem(n,h,filled);fs=10;ts=1/fs;n=-2:ts:2;f=cos(1*pi*n);subplot(3,1,1),stem(n,f,filled);g=cos(14*pi*n);subplot(3,1,2),stem(n,g,filled);h

2、=cos(2600000*pi*n);subplot(3,1,3),stem(n,h,filled);M2-1(1) x=1 -3 4 2 0 -2 ;h=3 0 1 -1 2 1;L=length(x)+length(h)-1;XE=fft(x,L);HE=fft(h,L);y1=ifft(XE.*HE);y=1 0 -1 0 1 0;g=1 3 9 27 81 243M=length(y)+length(g)-1;YE=fft(y,M);GE=fft(g,M);y2=ifft(YE.*GE);a=0:L-1;subplot(2,1,1);stem(a,real(y1);b=0:M-1;su

3、bplot(2,1,2);stem(b,real(y2);M2-2N=10;k=0:N-1;f=cos(pi*k/20);F=fft(f);plot(k/10,abs(F),o);hold onset(gca,xtick,0,0.25,0.5,0.75,1);set(gca,ytick,0,2,4,6,8);grid onhold offM2-3N=input(lenth of signal:N=);M=input(points of DFT:M=);k = 0:N-1;f = cos(2*pi*100*k/600)+cos(2*pi*150*k/600); %0.15*w=hamming(l

4、ength(f);f=f*w; %加hamming窗F= fft(f, M);L = 0:(M-1);plot(L/M,abs(F)grid on;xlabel(Normalized frequency);ylabel(Magnitude);M2-4M=4;N=64;n=-(N-1)/2:(N-1)/2;x=cos(2*pi*n/15)+0.75*cos(2.3*pi*n/15);X=fft(x,N);subplot(2,1,1);stem(n,fftshift(x);ylabel(xn);xlabel(Time n);subplot(2,1,2);stem(omega,real(fftshi

5、ft(X);ylabel(Xk);xlabel(Frequency(rad);hold;M2-5w=linspace(0,10,1024);plot(w,2./(w.2+1);% sampling points and frequency(rad/s)N=input(抽样点数=);Ws=input(抽样角频率=);Ts=2*pi/Ws;%compute the sampling pointsk=0:N/2;t=k*Ts;f1=exp(-3*t);f1(N/2+1)=2*f1(N/2+1);f2=f1(2:N/2);f=f1 fliplr(f2);F=Ts*real(fft(f);w=k*Ws/

6、N;w1=linspace(0,Ws/2,512);plot(w1,2./(w1.2+1),w,F(1:N/2+1),r);axis(0 10 0 2.2);xlabel(频率(秒/弧度);ylabel(幅度);z=N= num2str(N) Ws= num2str(Ws) 的结果; legend(理论值,z);title(exp(-|t|)的谱);M4-1wp=10;ws=2;Ap=1;As=40;N,Wc=cheb1ord(wp,ws,Ap,As,s);num,den=cheby1(N,Ap,Wc,s);disp(LP 分子多项式);fprintf(%.4en,num);disp(LP 分

7、母多项式);fprintf(%.4en,den);numt,dent=lp2hp(num,den,1);disp(HP 分子多项式);fprintf(%.4en,numt);disp(LP 分母多项式);fprintf(%.4en,dent);M4-2wp=1;ws=3.3182;Ap=1;As=32;w0=sqrt(48);B=2;N,Wc=buttord(wp,ws,Ap,As,s);num,den=butter(N,Wc,s);numt,dent=lp2bp(num,den,w0,B);w=linspace(2,12,1000);h=freqs(numt,dent,w);plot(w,2

8、0*log10(abs(h);grid;xlabel(Frequency in rad/s);ylabel(Gain in dB);M4-3wp=1;ws=3.3182;Ap=1;As=32;w0=sqrt(48);B=2;N,Wc=ellipord(wp,ws,Ap,As,s);num,den=ellip(N,Ap,As,Wc,s);numt,dent=lp2bp(num,den,w0,B);w=linspace(2,12,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h);grid;xlabel(Frequency in rad/s);yla

9、bel(Gain in dB);M4-4Ap=1;As=10;wp1=6;wp2=13;ws1=9;ws2=11;B=ws2-ws1;w0=sqrt(ws1*ws2);wLp1=B*wp1/(w0*w0-wp1*wp1);wLp2=B*wp2/(w0*w0-wp2*wp2);wLp=max(abs(wLp1),abs(wLp2);N,wc=ellipord(wp,ws,Ap,As,s);num,den=ellip(N,Ap,As,Wc,s);numt,dent=lp2bs(num,den,w0,B);w=linspace(5,35,1000);h=freqs(numt,dent,w);plot

10、(w,20*log10(abs(h);w=wp1 ws1 ws2 wp2;set(gca,xtick,w);grid;h=freqs(numt,dent,w);grid;h=freqs(numt,dent,w);A=-20*log10(abs(h)M4-5Wp=0.1*pi;Ws=0.4*pi;Ap=1;As=25;Fs=1;wp=Wp*Fs;ws=Ws*Fs;N=cheb2ord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/2/N);numa,dena=cheby2(N,As,wc,s);numd,dend=impinvar(numa,dena,Fs);w=l

11、inspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w*pi,20*log10(abs(h)/norm)w=Wp Wsh=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log10(abs(h(1);fprintf(Ap=%.4fn,-20*log10(abs(h(2);M4-6M4-7Wp=0.2*pi;Ws=0.4*pi;Ap=1;As=15;Fs=1;wp=Wp*Fs;ws=Ws*Fs;N=cheb1ord(wp,ws,Ap,As,s);wc=wp/(10(0

12、.1*Ap)-1)(1/2/N);numa,dena=cheby1(N,Ap,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w*pi,20*log10(abs(h)/norm)w=Wp Wsh=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log10(abs(h(1);fprintf(Ap=%.4fn,-20*log10(abs(h(2);M4-8Wp=4*pi/44.1;W

13、s=20*pi/44.1;Ap=0.5;As=50;Fs=44100;wp=Wp*Fs;ws=Ws*Fs;N=buttord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h/norm);xlabel(f);ylabel(Gain);w=Wp Ws;

14、h=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log10(abs(h(1);fprintf(As=%.4fn,-20*log10(abs(h(2);Wp=4*pi/44.1;Ws=20*pi/44.1;Ap=0.5;As=50;T=2;Fs=1/T;wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;N,wc=buttord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);numa,dena=butter(N,wc,s);numd,dend=bilinear(numa,dena,Fs);w=linspa

15、ce(0,pi,1024);h=freqz(numd,dend,w);plot(w/pi,20*log10(abs(h);axis(0 1 -50 0);xlabel(f);ylabel(Gain);w=Wp Ws;h=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log10(abs(h(1);fprintf(As=%.4fn,-20*log10(abs(h(2);M4-9Wp=0.8*pi;Ws=0.6*pi;Ap=0.5;As=30;T=2;Fs=1/T;wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;N,wc=buttord(wp,w

16、s,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);numa,dena=butter(N,wc,s);numt,dent=lp2hp(numa,dena,1);numd,dend=bilinear(numt,dent,Fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);plot(w/pi,20*log10(abs(h);axis(0 1 -50 0);xlabel(f);ylabel(Gain);w=Wp Ws;h=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log10(abs(h(1);f

17、printf(As=%.4fn,-20*log10(abs(h(2);M5-1function FIR_LPclearclcwp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;N=ceil(6.6*pi/tr_width)+1;n=0:1:N-1;wc=(ws+wp)/2;%理想LPF截止频hd=ideal_lp(wc,N);%理想低通滤波器计算,hd=0-(N-1)之间的理想脉冲响应,wc截止频率（弧度）,N=理想滤波器的长度wd1=hanning(N);b1=hd.*wd1;wd2=hamming(N);b2=hd.*wd2;wd3=blackman(N);b3=hd.*

18、wd3;wd4=kaiser(N);b4=hd.*wd4;db1,mag1,pha1,w=freqz_m(b1,1);%汉宁窗delta_w=2*pi/1000;rp1=-(min(db1(1:1:wp/delta_w+1);%实际带通波动rp1as1=-round(max(db1(ws/delta_w+1:1:501);%最小带阻衰减as1db2,mag2,pha2,w=freqz_m(b2,1);%海明窗delta_w=2*pi/1000;rp2=-(min(db2(1:1:wp/delta_w+1);%实际带通波动rp2as2=-round(max(db2(ws/delta_w+1:1:

19、501);%最小带阻衰减as2db3,mag3,pha3,w=freqz_m(b3,1);%blackmandelta_w=2*pi/1000;rp3=-(min(db3(1:1:wp/delta_w+1);%实际带通波动rp3as3=-round(max(db3(ws/delta_w+1:1:501);%最小带阻衰减as3db4,mag4,pha4,w=freqz_m(b4,1);%kaiserdelta_w=2*pi/1000;rp4=-(min(db4(1:1:wp/delta_w+1);%实际带通波动rp4as4=-round(max(db4(ws/delta_w+1:1:501);%

20、最小带阻衰减as4figure(1)stem(n,hd);title(理想脉冲响应)axis(0 N-1 -0.1 0.3);xlabel(n);ylabel(hd(n);figure(2)subplot(2,2,1)plot(w,mag1,:b)legend(汉宁窗低通滤波器)subplot(2,2,2)plot(w,mag2,-.g)legend(海明窗低通滤波器)subplot(2,2,3)plot(w,mag3,-r)legend(布来克曼窗低通滤波器)subplot(2,2,4)plot(w,mag4,-c)legend(凯泽窗低通滤波器)figure(3)plot(w,mag1,:

21、b,w,mag2,-.g,w,mag3,-r,w,mag4,-c)legend(汉宁窗低通滤波器,海明窗低通滤波器,布来克曼窗低通滤波器,凯泽窗低通滤波器)figure(4)plot(w/pi,20*log10(mag1),:b,w/pi,20*log10(mag2),-.g,w/pi,20*log10(mag3),-r,w/pi,20*log10(mag4),-c)legend(汉宁窗幅度响应（dB）,海明窗幅度响应（dB）,布来克曼窗幅度响应（dB）,凯泽窗幅度响应（dB）)figure(5)plot(n,b1,:b,n,b2,-.g,n,b3,-r,n,b4,-c)legend(汉宁窗h

22、(n),海明窗h(n),布来克曼窗h(n),凯泽窗h(n)M5-2Wp=0.6*pi;Ws=0.4*pi;Ap=1;As=45;N=ceil(7*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1;w=hamming(N);Wc=(Wp+Ws)/2;k=0:M;hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)+1;h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);plot(omega/pi,20*log10(abs(mag);M5-4N=40;alfa=(40

23、-1)/2;k=0:N-1;w1=(2*pi/N)*k;T1=0.109021; T2=0.59417456;hrs=zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4);hdr=0,0,1,1,0,0; wd1=0,0.2,0.35,0.65,0.8,1;k1=0:floor(N-1)/2); k2=floor(N-1)/2)+1:N-1;angH=-alfa*(2*pi)/N*k1,alfa*(2*pi/N*(N-k2);H=hrs.*exp(j*angH);h=real(ifft(H);d

24、b,mag,pha,grd,w = freqz_m(h,1);Hr,ww,a,L =Hr_Type2(h);figure(1)subplot(2,2,1)plot(w1(1:21)/pi,hrs(1:21),o,wd1,hdr)axis(0,1,-0.1,1.1);title(带通：N=40，T1=0.109021， T2=0.59417456)ylabel(Hr(k);set(gca,XTickMode,manual,XTick,0,0.2,0.35,0.65,0.8,1)set(gca,YTickMode,manual,YTick,0,0.059,0.109,1);grid %绘制带网格的

25、图像subplot(2,2,2); stem(k,h);axis(-1,N,-0.4,0.4)title(脉冲响应); ylabel(h(n); text(N+1,-0.4,n)subplot(2,2,3); plot(ww/pi,Hr,w1(1:21)/pi,hrs(1:21),o);axis(0,1,-0.1,1.1);title(振幅响应)xlabel(频率 （单位：pi）); ylabel(Hr(w)set(gca,XTickMode,manual,XTick,0,0.2,0.35,0.65,0.8,1);set(gca,YTickMode,manual,YTick,0,0.059,0

26、.109,1);gridsubplot(2,2,4); plot(w/pi,db); axis(0,1,-100,10);gridtitle(幅度响应);xlabel(频率 （单位：pi）); ylabel(分贝)set(gca,XTickMode,manual,XTick,0,0.2,0.35,0.65,0.8,1)set(gca,YTickMode,manual,YTick,-60;0);set(gca,YTickLabelMode,manual,YTickLabels,60;0);s,fs,nbits=wavread(sj.wav);%信号de 取样频率为44100HZx=s(:,1);

27、sound(x,fs);L=length(x);f=fs*(0:L-1)/L;t=0:1/fs:(L-1)/fs; %将所加噪声信号的点数调整到与原始信号相同%Au=1d=0.03*abs(max(x)*cos(2*pi*22000*t); %噪声为500和3300Hz的余弦信号%dz=cos(0.5*pi*fs*t);%载波dz=cos(2*pi*11025*t);xd=x.*dz;xz=xd+d;sound(xz,fs); %播放加噪声后的语音信号X=fft(x); %求信号的频谱XD=fft(xd); %信号调制后的频谱XZ=fft(xz);figure(2)subplot(3,1,1)

28、;plot(t,x)title(未加噪的信号); xlabel(time s);ylabel(幅度 );subplot(3,1,2);plot(t,xd)title(调制后的信号); xlabel(time s);ylabel(幅度 );subplot(3,1,3);plot(t,xz)title(调制加噪后的信号); xlabel(time n);ylabel(fuzhi n);figure(3)subplot(3,1,1);plot(f,abs(X);title(原始语音信号频谱);xlabel(频率（单位：Hz）);ylabel(幅度);subplot(3,1,2);plot(f,abs(XD);title(调制后的信号频谱);xlabel(频率（单位：Hz）);ylabel(幅度);subplot(3,1,3);plot(f,abs(XZ);title(加噪后的信号频谱);xlabel(频率（单位：Hz）);ylabel(幅度);y = fftfilt(h,xd);Y=fft(y);sound(3*y,