【正文】 。title(39。原信號頻域39。)。xlabel(39。頻率(Hz)39。)。ylabel(39。幅度39。)。subplot(4,1,3)。plot(td,xd)。title(39。延時信號時域圖39。)xlabel(39。時間(s)39。)。ylabel(39。幅度39。)。subplot(4,1,4)。plot([fs/2:dFs: fs/2dFs],abs(Xd))。title(39。延時信號頻域圖39。)xlabel(39。頻率(Hz)39。)。ylabel(39。幅度39。)。figure(2)。subplot(2,1,1)。plot(td,x1)。title(39。混響信號時域圖39。)。xlabel(39。時間(s)39。)。ylabel(39。幅度39。)。figure(2)。subplot(2,1,2)。plot([fs/2:dFs: fs/2dFs],abs(Xh2))。title(39?；祉懶盘栴l域圖39。)。xlabel(39。頻率(Hz)39。)。ylabel(39。幅度39。)。sound(x)。sound(xd)。sound(x1)。第三題、數(shù)字濾波器的設(shè)計及實(shí)現(xiàn)：（1）調(diào)用信號產(chǎn)生函數(shù)mstg產(chǎn)生由三路抑制載波調(diào)幅信號相加構(gòu)成的復(fù)合信號st，觀察st的時域波形和幅頻特性曲線；（2）要求將st中的三路調(diào)幅信號分離，通過觀察st的幅頻特性曲線，分別確定可以分離st中三路抑制載波單頻調(diào)幅信號的三個濾波器（低通濾波器、帶通濾波器、高通濾波器）的通帶截止頻率和阻帶截止頻率，阻帶最小衰減為60dB；（3）編程調(diào)用MATLAB濾波器設(shè)計函數(shù)分別設(shè)計這三個數(shù)字濾波器，并繪圖顯示其幅頻特性曲線；（4）調(diào)用濾波函數(shù)filter，用所設(shè)計的三個濾波器分別對復(fù)合信號st進(jìn)行濾波，分離出st中的三路不同載波頻率的調(diào)幅信號，并繪圖顯示濾波后信號的時域波形和頻譜，觀察分離效果。：（1）產(chǎn)生三路調(diào)幅信號；（2）三路信號相加組成復(fù)合信號；（3）分別設(shè)計低通、帶通、高通濾波器，對三路信號濾波；（4）繪出相應(yīng)的時域、頻域圖；（5）觀察分離效果；：（1）時域圖與頻域圖（第一路）：第二路：第三路：：N=1600。%N為信號st的長度。Fs=10000。T=1/Fs。Tp=N*T。%采樣頻率Fs=10kHz，Tp為采樣時間 t=0:T:(N1)*T。k=0:N1。f=k/Tp。fc1=Fs/10。%第1路調(diào)幅信號的載波頻率fc1=1000Hz, fm1=fc1/10。%第1路調(diào)幅信號的調(diào)制信號頻率fm1=100Hz fc2=Fs/20。%第2路調(diào)幅信號的載波頻率fc2=500Hz fm2=fc2/10。%第2路調(diào)幅信號的調(diào)制信號頻率fm2=50Hz fc3=Fs/40。%第3路調(diào)幅信號的載波頻率fc3=250Hz, fm3=fc3/10。%第3路調(diào)幅信號的調(diào)制信號頻率fm3=25Hz xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t)。%產(chǎn)生第1路調(diào)幅信號 xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t)。%產(chǎn)生第2路調(diào)幅信號 xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t)。%產(chǎn)生第3路調(diào)幅信號 st=xt1+xt2+xt3。%三路調(diào)幅信號相加 fxt=fft(st,N)。%計算信號st的頻譜 figure(1)。subplot(2,1,1)。plot(t,st)。title(39。三路調(diào)幅信號時域圖39。)grid。xlabel(39。t/s39。)。ylabel(39。s(t)39。)。axis([0,Tp/2,min(st),max(st)])。title(39。(a)s(t)的波形39。)。subplot(2,1,2)。stem(f,abs(fxt)/max(abs(fxt)),39。.39。)。grid。title(39。(b)s(t)的頻譜39。)。axis([0,Fs/5,0,])。xlabel(39。f/Hz39。)。ylabel(39。幅度39。)。Fs=10000。T=1/Fs。%采樣頻率%調(diào)用信號產(chǎn)生函數(shù)mstg產(chǎn)生由三路抑制載波調(diào)幅信號相加構(gòu)成的復(fù)合信號st%低通濾波器設(shè)計與實(shí)現(xiàn)fp=280。fs=450。wp=2*fp/Fs。ws=2*fs/Fs。rp=。rs=60。[N,wp]=ellipord(wp,ws,rp,rs)。[B,A]=ellip(N,rp,rs,wp)。y1t=filter(B,A,st)。[Hb,Wc]=freqz(B,A)。figure(2)。subplot(2,1,1)。plot(Fs*Wc/(2*pi),20*log10(abs(Hb)))。title(39。低通濾波器頻域響應(yīng)圖39。)subplot(2,1,2)。plot(t,y1t)。title(39。第一路低頻信號時域圖39。)%帶通濾波器設(shè)計與實(shí)現(xiàn)fpl=440。fpu=560。fsl=275。fsu=900。wp=[2*fpl/Fs,2*fpu/Fs]。ws=[2*fsl/Fs,2*fsu/Fs]。rp=。rs=60。[N,wp]=ellipord(wp,ws,rp,rs)。[B,A]=ellip(N,rp,rs,wp)。y2t=filter(B,A,st)。[Hb1,Wc1]=freqz(B,A)。figure(3)。subplot(2,1,1)。plot(Fs*Wc1/(2*pi),20*log10(abs(Hb1)))。title(39。帶通濾波器頻域響應(yīng)圖39。)subplot(2,1,2)。plot(t,y2t)。title(39。第二路中頻信號時域圖39。)%高通濾波器設(shè)計與實(shí)現(xiàn) fp=890。fs=600。wp=2*fp/Fs。ws=2*fs/Fs。rp=。rs=60。%DF指標(biāo)（低通濾波器的通、阻帶邊界頻）[N,wp]=ellipord(wp,ws,rp,rs)。%調(diào)用ellipord計算橢圓DF階數(shù)N和通帶截止頻率wp [B,A]=ellip(N,rp,rs,wp,39。high39。)。%調(diào)用ellip計算橢圓帶通DF系統(tǒng)函數(shù)系數(shù)向量B和A y3t=filter(B,A,st)。%濾波器軟件實(shí)現(xiàn) [Hb2,Wc2]=freqz(B,A)。figure(4)。subplot(2,1,1)。plot(Fs*Wc2/(2*pi),20*log10(abs(Hb2)))。title(39。高通濾波器頻域響應(yīng)圖39。)subplot(2,1,2)。plot(t,y3t)。title(39。第三路高頻信號時域圖39。)第四題、心電信號的處理：（1）在MATLAB軟件平臺下，給原始的心電信號疊加上噪聲或干擾，干擾類型分為如下幾種：白噪聲、工頻干擾（50Hz）、諧波干擾（二次、三次諧波為主，分別為100Hz、150Hz）繪出疊加噪聲后的心電信號時域和頻譜圖，在視覺上與原始心電信號圖形對比，分析頻域基本特征變化。（2）給定濾波器的規(guī)一化性能指標(biāo)（參考指標(biāo)，實(shí)際中依據(jù)每個同學(xué)所疊加噪聲情況而定），例如：通帶截止頻率wp=*pi, 阻通帶截止頻率ws=*pi。通帶最大衰減Rp=1 dB。阻帶最小衰減Rs=15 dB（3）采用窗函數(shù)法設(shè)計各型FIR濾波器（低通、高通、帶通、帶阻中的至少2種類型），來對疊加干擾前后的心電信號進(jìn)行濾波處理，繪出濾波器的頻域響應(yīng)及濾波后信號的時域波形和頻譜，并對濾波前后的信號進(jìn)行對比，分析信號的變化；在相同的性能指標(biāo)下比較各方法的濾波效果，并從理論上進(jìn)行分析或解釋；（4）采用雙線性變換法利用不同的原型低通濾波器（Butterworth型與切比雪夫I型）來設(shè)計各型IIR濾波器（低通、高通、帶通、帶阻中的至少2種類型）繪出濾波器的頻域響應(yīng)，并用這些數(shù)字濾波器對含噪心電信號分別進(jìn)行濾波處理；比較不同方法下設(shè)計出來的數(shù)字濾波器的濾波效果，并從理論上進(jìn)行分析或解釋；（5）心電信號波形觀察、頻譜觀察，對濾波后的心電信號觀察其時域、頻域特征變化。：（1）讀取原始心電信號；（2）心電信號加50HZ工頻信號；（3）心電信號加入諧波；（4）心電信號加入白噪聲；（5）設(shè)計相應(yīng)的濾波器進(jìn)行不同噪聲的濾除；（6）繪出相應(yīng)的時域頻域圖；：（1）各時域頻域圖：原始心電信號時域圖以及加入個噪聲后的時域圖：原始心電信號頻域圖以及加入個噪聲后的時頻域圖：（2）具體代碼實(shí)現(xiàn)：%原心電信號xl=load(39。E:39。)。x=xl(:,2)。N=size(x,1)。fs=1000。t=(0:length(x)1)/fs。f=fs/N*(0:(N+1)/21)+1。X=fft(x,N)。%加50hz工頻干擾 y=x+sin(2*50*pi*t)39。Y=fft(y,N)。%加100hz、150hz諧波干擾y1=x+sin(2*100*pi*t)39。+cos(2*150*pi*t)39。Y1=fft(y1,N)。%加白噪聲 yb=awgn(x,20)。Yd=fft(yb,N)。figure(1)。subplot(4,1,1)。plot(t,x)。title(39。原信號時域圖39。)xlabel(39。時間(s)39。)。ylabel(39。幅度39。)。figure(2)。subplot(4,1,1)plot(f,abs(X(1:(N+1)/2)))。title(39。原信號頻域圖39。)xlabel(39。頻率(Hz)39。)。ylabel(39。幅度39。)。figure(1)。subplot(4,1,2)。plot(t,y)。title(39。加50hz工頻信號時域圖39。)xlabel(39。時間(s)39。)。ylabel(39。幅度39。)。figure(2)。subplot(4,1,2)。plot(f,abs(Y(1:(N+1)/2)))。title(39。加50hz工頻信號頻域圖39。)xlabel(39。頻率(Hz)39。)。ylabel(39。幅度39。)。figure(1)subplot(4,1,3)。plot(t,y1)。title(39。加諧波后的時域圖39。)xlabel(39。時間(s)39。)。ylabel(39。幅度39。)。figure(2)subplot(4,1,3)。plot(f,abs(Y1(1:(N+1)/2)))。title(39。加諧波后的頻域圖39。)xlabel(39。頻率(Hz)39。)。ylabel(39。幅度39。)。figure(1)。subplot(4,1,4)。plot(t,y1)。title(39。加白噪后的時域圖39。)。figure(2)。subplot(4,1,4)。plot(f,abs(Y1(1:(N+1)/2)))。title(39。加白噪后的頻域圖39。)。（3）濾除部分：濾除50HZ信號后的時域圖（blackman低通濾波器）：濾除50HZ信號后的頻域圖： blackman低通濾波器濾波：具體代碼：%原心電信號xl=load(39。E:39。)。x=xl(:,2)。N=size(x,1)。fs=1000。t=(0:length(x)1)/fs。f=fs/N*(0:(N+1)/21)+1。%f=fs/N*(0:N1)。X=fft(x,N)。%加50hz工頻干擾 y=x+sin(2*50*pi*t)39。Y=fft(y,N)。%blackman低通濾波器 wp=30/fs*pi。ws=40/fs*pi。B=wswp。M=ceil(12*pi/B)1。bl=fir1(M,(ws+wp)/2/pi,39。low39。,blackman(M+1))。yl=fftfilt(bl,y)。Yl=fft(yl,N)。figure(1)。subplot(3,1,1)。plot(t,x)。title(39。原信號時域圖39。)figure(2)。subplot(3,1,1)plot(f,abs(X(1:(N+1)/2)))。%plot(f,abs(X))。title(39。原信號頻域圖39。)figure(1)。subplot(3,1,2)。plot(t,y)。title(39。加50hz工頻信號時域圖39。)figure(2)。subplot(3,1,2)。plot(f,abs(Y(1:(N+1)/2)))。%plot(f,abs(Y))。title(39。加50hz工頻信號頻域圖39。)figure(1)。subplot(3,1,3)。plot(t,yl)。title(39。濾波后信號時域圖39。)figure(2)subplot(3,1,3)。plot(f,abs(Yl(1:(N+1)/2)))。%plot(f,abs(Yl))。title(39。濾波后信號頻域圖39。)figure(3)[h2,w2]=freqz(bl,1)。plot(w2/(2*pi)*fs,20*log10(abs(h2)))。title(39。blackman低通濾波器頻率響應(yīng)圖39。)；濾除50HZ信號后的時域圖（切比雪夫?yàn)V波器）：濾除50HZ信號后的頻域圖（切比雪夫?yàn)V波器）：切比雪夫?yàn)V波器：具體代碼實(shí)現(xiàn)：%原心電信號xl=load(39。E:39。)。x=xl(:,2)。N=size(x,1)。fs=1000。t=(0:length(x)1)/fs。f=fs/N*(0:(N+1)/21)+1。X=fft(x,N)。%加50hz工頻干擾 y=x+sin(2*50*pi*t)39。Y=fft(y,N)。%切比雪夫?yàn)V波器wp=20/fs*2*pi。ws=30/fs*2*pi。wp1=2*fs*tan(wp/2)。ws1=2*fs*tan(ws/2)。Rp=1。Rs=15。[N1,w]=cheb1ord(wp1,ws1,Rp,Rs,39。s39。)。[B1,A1]=cheby1(N1,1,w,39。low39。,39。s39。)。[B,A]=bilinear(B1,A1,fs)。[Hk,wk]=freqz(B,A)。yf=filter(B,A,y)。Yf=fft(yf,N)。figure(1)。subplot(3,1,1)。plot(t,x)。title(39。原信號時域圖39。)figure(2)。subplot(3,1,1)plot(f,abs(X(1:(N+1)/2)))。title(39。原信號頻域圖39。)figure(1)。subplot(3,1,2)。plot(t,y)。title(39。加50hz工頻信號時域圖39。)figure(2)。subplot(3,1,2)。plot(