【文章內(nèi)容簡介】 分子Az=[1,zeros(1,R1),a]。%無限個回聲濾波器系統(tǒng)函數(shù)分母y2=filter(Bz,Az,y)。Y2=fft(y2,6001)。[h,w]=freqz(Bz,Az)。subplot(3,2,1)。plot(abs(h))。title(39。無限個回聲濾波器幅頻響應39。)。subplot(3,2,2)。 plot(angle(h))。title(39。無限個回聲濾波器相頻響應39。)。subplot(3,2,3:4)。plot(y2)。title(39。無限個回聲濾波器時域圖 39。)。subplot(3,2,5)。plot(abs(Y2))。title(39。無限個回聲濾波器信號幅頻 39。)。subplot(3,2,6)。plot(angle(Y2))。title(39。無限個回聲濾波器信號相頻 39。)。sound(y2,fs,bits)。 全通結構濾波器程序設計[y,fs,bits]=wavread(39。.wav39。)。y=y(:,1)。z=[zeros(3000,1)。y]。y=[y。zeros(3000,1)]。a=。R=5000。%濾波器階數(shù)設置，其值越高，回聲越明顯Bz1=[a,zeros(1,R1),1]。%全通濾波器系統(tǒng)函數(shù)分子Az1=[1,zeros(1,R1),a]。%全通濾波器系統(tǒng)函數(shù)分母[h,w]=freqz(Bz1,Az1)。yy2=filter(Bz1,Az1,y)。YY2=fft(yy2,6001)。subplot(3,2,1)。plot(abs(h))。title(39。全通濾波器幅頻響應39。)。subplot(3,2,2)。 plot(angle(h))。title(39。全通濾波器相頻響應39。)。subplot(3,2,3:4)。plot(yy2)。title(39。全通結構時域圖 39。)。subplot(3,2,5)。plot(abs(YY2))。title(39。全通結構幅頻 39。)。subplot(3,2,6)。plot(angle(YY2))。title(39。全通結構相頻 39。)。sound(yy2,fs,bits)。部分結果展示：數(shù)字濾波器的設計與實現(xiàn)調(diào)用函數(shù)mstg產(chǎn)生st，自動繪圖顯示st的時域波形和幅頻特性曲線調(diào)用ellipord和ellip分別設計三個橢圓濾波器，并繪圖顯示其幅頻響應特性曲線。調(diào)用filter，用三個濾波器分別對信號st進行濾波，分離出三路不同載波頻率的調(diào)幅信號y1(n)、y2(n)和y3(n)繪圖顯示y1(n)、y2(n)和y3(n)的時域波形和幅頻特性曲線結束通過觀察，三路信號的載波頻率分別為250，500，1000Hz,所以低通濾波器的截至頻率250Wc500，帶通濾波器的截至頻率250W1500,500W21000。高通濾波器的截至頻率500Wc1000。 調(diào)用ellipord計算橢圓DF階數(shù)N和通帶截止頻率wp. 調(diào)用ellip計算橢圓帶通DF系統(tǒng)函數(shù)系數(shù)向量B和A.產(chǎn)生信號N=1600。Fs=10000。T=1/Fs。Tp=N*T。 %采樣頻率Fs=10kHz，Tp為采樣時間t=0:T:(N1)*T。k=0:N1。f=k/Tp。T=1/Fs。t=0:T:(N1)*T。k=0:N1。fc1=Fs/10。fm1=fc1/10。fc2=Fs/20。fm2=fc2/10。fc3=Fs/40。fm3=fc3/10。xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t)。xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t)。xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t)。st=xt1+xt2+xt3。fxt=fft(st,N)。 %計算信號st的頻譜subplot(2,1,1)plot(t,st)。grid。xlabel(39。t/s39。)。ylabel(39。s(t)39。)。axis([0,Tp/8,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。)。%三路調(diào)幅信號的載波頻率分別為250Hz、500Hz、1000Hz。低通濾波fp=280。fs=450。%fp必須大于250，fs必須小于500wp=2*fp/Fs。ws=2*fs/Fs。rp=。rs=60。%DF指標（低通濾波器的通、阻帶邊界頻）[N,wp0]=ellipord(wp,ws,rp,rs)。%調(diào)用ellipord計算橢圓DF階數(shù)N和通帶截止頻率wp[B,A]=ellip(N,rp,rs,wp0)。%調(diào)用ellip計算橢圓帶通DF系統(tǒng)函數(shù)系數(shù)向量B和Ay1t=filter(B,A,st)。%濾波器軟件實現(xiàn)figure。subplot(2,1,1)。[H1,w]=freqz(B,A,1000)。%低通的頻率響應m=abs(H1)。plot(w/pi,20*log(m/max(m)))。grid on。title(39。低通濾波損耗函數(shù)曲線39。)。xlabel(39。w39。)。ylabel(39。H139。)。subplot(2,1,2)。plot(t,y1t)。title(39。低通濾波后的波形39。)。xlabel(39。t39。)。ylabel(39。y1t39。)。帶通濾波Fs=10000。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,wp0]=ellipord(wp,ws,rp,rs)。[B,A]=ellip(N,rp,rs,wp0)。y2t=filter(B,A,st)。figure。subplot(2,1,1)。[H2,w]=freqz(B,A,1000)。m=abs(H2)。plot(w/pi,20*log(m/max(m)))。grid on。title(39。帶通濾波損耗函數(shù)曲線39。)。xlabel(39。w39。)。ylabel(39。H239。)。subplot(2,1,2)。plot(t,y2t)。title(39。帶通濾波后的波形39。)。xlabel(39。t39。)。ylabel(39。y2t39。)。高通濾波Fs=10000。fp=550。fs=900。wp=2*fp/Fs。ws=2*fs/Fs。rp=。rs=60。[N,wp0]=ellipord(wp,ws,rp,rs)。[B,A]=ellip(N,rp,rs,wp0,39。high39。)。y3t=filter(B,A,st)。figure。subplot(2,1,1)。[H3,w]=freqz(B,A,1000)。m=abs(H3)。plot(w/pi,20*log(m/max(m)))。grid on。title(39。高通濾波損耗函數(shù)曲線39。)。xlabel(39。w39。)。ylabel(39。H339。)。subplot(2,1,2)。plot(t,y3t)。title(39。高通濾波后的波形39。)。xlabel(39。t39。)。ylabel(39。y3t39。)。部分結果展示：心電信號的處理開始向心電信號中加入白噪聲，50HZ工頻干擾，帶限chirp噪聲使用Butterworth型低通數(shù)字濾波器對信號過濾使用切比雪夫I型數(shù)字低通濾波器對信號過濾使用帶阻濾波器對信號過濾進行效果比較結束一、設計原理（1）、巴特沃斯低通濾波器巴特沃斯低通濾波器的平方幅度響應為其中，n為濾波器的階數(shù)，為低通濾波器的截止頻率。 該濾波器具有 一些特殊的性質(zhì)：① 對所有的n,都有當時, ； ②對所有的n,都有當時, ；③是的單調(diào)遞減函數(shù),即不會出現(xiàn)幅度響應的起伏； ④當時,巴特沃斯濾波器趨向于理想的低通濾波器； ⑤在處平方幅度響應的各級導數(shù)均存在且等于0，因此在該點上取得最大值，且具有最大平坦特性。 圖1展示了2階、4階、8階巴特沃斯低通濾波器的幅頻特性?？梢婋A數(shù)n越高，其幅頻特性越好，低頻檢測信號保真度越高，過渡帶變窄，即衰減加劇，但半功率點不變。圖1 巴特沃斯低通濾波器的幅頻特性（2） 切比雪夫I型數(shù)字低通濾波器（1）確定數(shù)字低通濾波器的技術指標：通帶截止頻率ωp、通帶衰減ap、阻帶截止頻率ωs、阻帶衰減as切比雪夫濾波器的振幅平方特性如圖2所示：圖2 切比雪夫濾波器的振幅平方特性(2)將數(shù)字低通濾波器的技術指標轉換成模擬低通濾波器的技術指標。如果采用脈沖響不變法，邊界頻率的轉換關系為:如果采用雙線性變換法，邊界頻率的轉換關系為(3) 按照模擬低通濾波器的技術指標設計模擬低通濾波器。(4) 利用雙線性變換法將模擬濾波器Ha(s)，從s平面轉換到z平面，得到數(shù)字低通濾波器系統(tǒng)函數(shù)H(z)。（5）數(shù)字低通技術指標為：ωp=, a p=1dB。 ωs=, a s=40Db（6）模擬低通的技術指標為：歸一化截止角頻率wp=2pi*Fs/Ft。 ws=2pi*Fs/Ft(7)利用模擬切比雪夫濾波器設計數(shù)字濾波器。通帶截止頻率為：wp=*pi。 阻帶截止頻率為：ws=*pi。通帶最大衰減為：Rp=1。阻帶最大衰減為：As=15。設定周期為1s；模擬低通濾波器的生成：[b,a]=cheby1(n,1,Wn,39。low39。,39。s39。)。滿足設計指標的最小階數(shù)和截止頻率：Wn[n,Wn]=cheb1ord(OmegaP,OmegaS,1,40,39。s39。)。最后實現(xiàn)輸入輸出、幅頻特性、相頻特性的圖形。（3） IIR數(shù)字濾波器的性質(zhì)無限長沖激響應( IIR) 數(shù)字濾波器是數(shù)字濾器的一種，數(shù)字濾波器根據(jù)其沖激響應函數(shù)的時域特性,還包括有限長沖激響應( FIR)數(shù)字濾波器。IIR數(shù)字濾波器的特征是：①具有無限持續(xù)時間沖激響應；②需要用遞歸模型來實現(xiàn)，這可以從其差分方程得出，也可以從其系數(shù)函數(shù)為：得出。 數(shù)字巴特沃思濾波器屬于IIR濾波器，該類濾波器具有特定的性質(zhì)和設計方法。目前比較成熟的IIR數(shù)字濾波器設計方法有兩種：1）直接法目前所用的方法主要是：零極點累試法、頻域幅度平方誤差最小法和時域單位脈沖響應逼近法。直接法的最大優(yōu)點在于可以設計任意幅頻特征的濾波器。2)間接法，目前所用的方法主要是：沖激響應不變法、階躍響應不變法和雙線性法。它們都是借助于 已經(jīng)成熟的現(xiàn)有低通濾波器原型進行設計，即對數(shù)字低通數(shù)字濾波器，先將數(shù)字低通濾波器的技術指標按希望的設計方法轉換為模擬低通濾波器的技術指標，再按指定的模擬低通濾波器的類型設計模擬濾波器H(s)，然后，將模擬濾波器的系統(tǒng)函數(shù)H(s)從s平面轉換到z平面，得到數(shù)字低通濾波器的系統(tǒng)函數(shù)H(z)；如所設計的數(shù)字濾波器為高通、帶通或帶阻濾波器，則可借助模擬濾波器的頻帶變換轉換為低通模擬濾波器。 由于直接法設計巴特沃思濾波器相對復雜，在不需要任意幅頻特征的情況下，一般采用問接法，同時由于沖激響應不變法和脈沖響應不變法，從s平面轉換到z平面的映射為多值映射，容易造成頻譜混疊，故而本文采用不會產(chǎn)生頻譜混疊的雙線性變換法。 二、 課程設計步驟及結果分析（1）心電數(shù)據(jù)的導入將老師給的心電信號原始數(shù)據(jù)存于桌面，然后我選用桌面上的第一組數(shù)據(jù)a=load(39。C:\Documents and Settings\hnie\桌面\心電信號數(shù)據(jù)\39。)。 （2）繪出心電信號的時域圖和頻譜圖將導入的兩行數(shù)據(jù)分別用t,b來替換，然后通過調(diào)用plot函數(shù)來畫出時域圖，然后通過對1000個心電數(shù)據(jù)的幅值進行FFT運算，再次調(diào)用plot函數(shù)來繪出頻域圖，具體設計如下：figure(1)。subplot(2,1,1)。t=a(1:1000,1)。b=a(1:1000,2)。plot(t,b)。title(39。原始波形圖39。)。xlabel(39。時間（s）39。)。ylabel(39。幅值（A）39。)。y1=fft(a(:,2),1000)。f1=100*(0:999)/1000。 subplot(2,1,2)。plot(f1,abs(y1))。title(39。原始頻譜圖39。)。xlabel(39。頻率（Hz）39。)。ylabel(39。幅度(dB)39。)。圖通過導入的心電信號數(shù)據(jù)發(fā)現(xiàn)在其頻譜圖上的0~20H z和80~100Hz之間的幅值比較大，而在30~70Hz之間的幅值相對較小。（3） 加入噪聲干擾這里我加入的噪聲是：白噪聲，50Hz的工頻噪聲，帶限chirp噪聲。（4）白噪聲通過用s來代表加入白噪聲后的信號，并進行數(shù)字濾波器的頻率響應，對s中的1000個頻率點調(diào)用plot函數(shù)畫出加入白噪聲后的時域圖，再對ws/pi,abshs調(diào)用plot函數(shù)畫出加入白噪聲后的頻譜圖，具體操作如下：q=*rand(1000,1)。s=a(:,2)+q。[hs,ws]=freqz(s,1,1024)。abshs=abs(hs)。figure(2)。subplot(2,1,1)。plot(s(1:1000))。title(39。加入白噪聲后的時域圖39。)。xlabel(‘時間（s）39。)。ylabel(39。幅值（A）39。)。subplot(2,1,2)。plot(ws/pi,abshs)。title(39。加入白噪聲后的頻譜圖39。)。xlabel(39。幅度39。)。ylabel(39。Hz39。)。 加入白噪聲后的時域圖和頻譜圖通過觀察加入白噪聲后的時域圖和頻域圖，將它與未加入白噪聲進行比較，可以