去除干擾蜂鳴音

一、課程設(shè)計(jì)題目去除干擾蜂鳴音1?目的：掌握信號(hào)時(shí)頻域分析方法，正確理解采樣定理，準(zhǔn)確理解濾波器的概念。2.內(nèi)容：提供一個(gè)包含某人說話語音片段的聲音文件，但該語音信號(hào)被一個(gè)包含有幾個(gè)諧波分量的蜂鳴信號(hào)十?dāng)_了。用Matlab的wavread命令讀取該聲音文件。注意，該命令可以同時(shí)得到聲音文件的采樣率和采樣位寬，請(qǐng)查閱Matlab的幫助文件。用快速傅立葉變換(FFT)計(jì)算并畫出聲音信號(hào)的頻譜，列寫出蜂鳴信號(hào)的諧波頻率。思考如何將這些蜂鳴音去除將去除了蜂鳴音的語音片段播放出來，仔細(xì)聆聽并寫下語音片段中人物所說的話。注意：由于只能播放實(shí)信號(hào)，因此記得提取信號(hào)的實(shí)部。Matlab命令：wavread,wavplay,fft,fftshift,firl,filter,plot,figure.二、 設(shè)計(jì)思路用waveread()函數(shù)讀取音頻和其采樣率和采樣位寬，對(duì)讀取的音頻信號(hào)使用fft()函數(shù)進(jìn)行快速傅立葉變換并繪出得到的頻譜。觀察頻譜分析噪聲(蜂鳴信號(hào))的諧波頻率分布，選擇合適的濾波模式將噪聲信號(hào)的諧波濾去，便可以得到去除噪聲后的人聲。設(shè)計(jì)濾波器的頻域特性便成了除去噪聲并留下原聲的關(guān)鍵，我們注意到所學(xué)的采樣定理以及一維sinc函數(shù)(辛格函數(shù)疽”偵=誓，然而漢寧窗可以看作是3個(gè)矩形時(shí)間窗的頻譜之和，或者說是3個(gè)Sinc(x)型函數(shù)之和，而括號(hào)中的兩項(xiàng)相對(duì)于第一個(gè)譜窗向左、右各移動(dòng)了n/T，從而使旁瓣互相抵消，消去高頻十?dāng)_和漏能。它適用于非周期性的連續(xù)信號(hào)。利用它的選擇特性使用fir1()建立一個(gè)漢寧窗函數(shù)，并用filter()函數(shù)進(jìn)行濾波，去除噪聲部分。最后用play()函數(shù)播放音頻檢查效果。三、 設(shè)計(jì)過程音頻的讀取和分析先將原始音頻文件讀入，[audioO,Fs,nbits]=wavread('C:\Users\Administrator\Desktop\signal\')；%按路徑讀取音頻存入audioO變量，并用Fs變量記錄采樣率，nbits變量記錄采樣位寬。其中，F(xiàn)s=11025 #采樣率為11025Hznbits=32 #采樣帶為32pO=audioplayer(audioO,Fs);%將audioO載入音頻播放器play(pO)；%并進(jìn)行播放subplot(2,1,1);%分屏繪圖plot(audio0);%繪制原始音頻時(shí)域圖，如下圖所示title('時(shí)域');%標(biāo)注題目[k]=fft(audio0,Fs);%對(duì)原始音頻進(jìn)行快速傅立葉變換subplot(2,1,2);plot(abs(k));%繪制原始音頻頻域圖title('頻域');頻域圖如圖1下圖所示圖1原聲時(shí)域頻域圖像此時(shí)，在時(shí)域中雜亂分布的聲音信號(hào)變換到頻域中將噪聲諧波所分布的頻域中顯得尤為清晰簡潔，這樣就不難發(fā)現(xiàn)干擾信號(hào)主要分布在0--1000HZ以內(nèi)，由于聽到的噪聲比人聲大得多，我們又畫出了，聲音強(qiáng)度的時(shí)域和頻域圖像，figure;subplot(2,1,1);plot(audio0);plot(20*log10(abs(audio0)/max(abs(audio0))));%繪制原始音頻聲音強(qiáng)度的頻域圖ylabel('分貝/dB');xlabel(時(shí)間/ms');title('時(shí)域');[k]=fft(audio0,Fs);subplot(2,1,2);plot(abs(k));plot(20*log10(abs(abs(k))/max(abs(abs(k)))));%繪制原始音頻聲音強(qiáng)度的頻域圖ylabel('分貝/dB');xlabel('頻率/Hz');title('頻域');figure;plot(abs(k));ylabel('振幅/A')；xlabel('頻率/Hz');set(gca,'XLim',[01000]);set(gca,'XTick',[0:20:1000]);得到如圖2所示的聲音強(qiáng)度的時(shí)域和頻域圖像2DA0SDCQ1MDQ12DD0set(gca,'XTick',[0:20:1000]);得到如圖2所示的聲音強(qiáng)度的時(shí)域和頻域圖像2DA0SDCQ1MDQ12DD0圖2原聲聲音強(qiáng)度時(shí)域頻域圖像從頻域圖中我們?nèi)匀话l(fā)現(xiàn)了四個(gè)異常尖峰，再次明顯的證實(shí)了干擾信號(hào)主要分布在0--1000HZ以內(nèi)；為了跟清晰地和觀察干擾信號(hào)的頻域分布情況，我們進(jìn)一步繪制這一范圍的圖像，進(jìn)行局部深入研究。plot(abs(k));%重繪圖像set(gca,'XLim',[01000]);%更改顯示范圍為0-1000set(gca,'XTick',[0:20:1000]);%更改坐標(biāo)步長為20圖3噪聲局部頻域放大圖觀察圖像可發(fā)現(xiàn)，干擾信號(hào)的諧波頻率為：220Hz、440Hz、660Hz、880Hz，在放大后發(fā)現(xiàn)噪聲信號(hào)為4個(gè)幾乎對(duì)稱的三角波，而非單位沖擊，結(jié)合以上幾個(gè)圖，我們分析得到干擾信號(hào)主要分布在01000Hz以內(nèi)，而人聲是大部分分布于大于1000Hz區(qū)間的，由此我們想到了對(duì)頻率具有選擇特性的濾波器，且此處需要設(shè)計(jì)一個(gè)高通濾波器，以將位于01000Hz以內(nèi)的噪聲濾掉，留下大于1000Hz人聲信號(hào)。濾波器的設(shè)計(jì)由于干擾信號(hào)的諧波頻率為：220Hz、440Hz、660Hz、880Hz，而人類說話的頻率大概在3003400Hz，而干擾信號(hào)非常大，需要一個(gè)濾波器來實(shí)現(xiàn)將大約高于1000Hz的信號(hào)保留，低于1000Hz的信號(hào)濾掉，觀察分貝圖，發(fā)現(xiàn)大部分噪聲分布在40dB以內(nèi)，因此阻帶最小衰減不應(yīng)小于40dB表1常用窗函數(shù)的特點(diǎn)窗函數(shù)主瓣寬旁瓣峰值/dB阻帶最小衰減/dB矩形窗4qr/N-13-21三角形窗M/N-25-25漢寧窗阮小-31-44哈明窗Stt/N-41-53kaiser窗12M-57-80根據(jù)上表顯示各種窗函數(shù)的參數(shù)特點(diǎn)，選擇hanning窗濾波，利用其可以使旁瓣互相抵消頻域特性，據(jù)此可設(shè)定合適的參數(shù)設(shè)計(jì)一個(gè)hanning窗函數(shù)高通濾波器。fp=1000;fs=900;%通帶頻率fp，阻帶頻率fswp=2*fp*pi/Fs;ws=2*fs*pi/Fs;%歸一化邊界頻率wc=(wp+ws)/2/pi;%歸一化中心頻率wdp=wpws;#過渡帶寬N=ceil(12*pi/wdp);%由窗函數(shù)主瓣寬和過渡帶寬，求得窗函數(shù)最小長度N=N+mod(N,2);%高通濾波器N必為奇數(shù)HPfir=fir1(N,wc,'high',hanning(N+1));%設(shè)計(jì)高通hanning窗濾波器HPfir該濾波器的主要參數(shù)為：通帶邊界為1000Hz，阻帶邊界為900Hz，阻帶衰減不小于40dB。接下來在時(shí)域和頻域直觀地展示其濾波特性，再根據(jù)濾波效果對(duì)其參數(shù)做微調(diào)，figure;%新建圖像subplot(2,1,1);plot(HPfir);%繪制濾波器時(shí)域圖像title('濾波器時(shí)域');plot(abs(fft(HPfir)));%繪制濾波器頻域圖像title('濾波器頻域');得到圖4所示的濾波器時(shí)域和頻域特性圖。圖4基于漢寧函數(shù)的高通濾波器時(shí)域頻域圖由圖4可見該濾波器的截止頻率大約在900Hz--1000Hz之間，完全符合設(shè)計(jì)的目的，濾波器設(shè)計(jì)完成之后，對(duì)原始音頻信號(hào)進(jìn)行濾波處理：audio1=filter(HPfir,1,audio0);%使用filter函數(shù)對(duì)原聲做一維數(shù)字濾波p1=audioplayer(audio1,Fs);figure;%新建濾波后的圖像subplot(2,1,1);plot(audio1);%繪制濾波后的時(shí)域圖像title('濾波后時(shí)域');[k0]=fft(audio1,Fs);%對(duì)濾波后的信號(hào)做快速傅里葉變換subplot(2,1,2);plot(abs(k0));%繪制濾波后的頻域圖像title('濾波后頻域');此時(shí)，我們得到了如圖5所示的濾波后的時(shí)域頻域圖像圖5濾波后時(shí)域頻域圖最后由于處理后的聲音信號(hào)幅度較小，聽不清晰，需要對(duì)音頻信號(hào)進(jìn)行增幅處理。audio1=audio1*10;%增幅處理p1=audioplayer(audio1,Fs);play(pl);%播放處理后的音頻至此，設(shè)計(jì)結(jié)束，我們獲得了去除噪音后較為清晰的、完整的人聲信號(hào)。四、源代碼[audio。，F(xiàn)s,nbits]=wavread('C:\Users\xufanyun\Desktop\signal\');p0=audioplayer(audio0,Fs);subplot(2,1,1);plot(audio0);ylabel('振幅/A');xlabel(時(shí)間/ms');title('時(shí)域');[k]=fft(audio0,Fs);subplot(2,1,2);plot(abs(k));ylabel（'振幅/A'）;xlabel('頻率/Hz');title('頻域');figure;subplot(2,1,1);plot(audio0);plot(20*log10(abs(audio0)/max(abs(audio0))));ylabel('分貝/dB');xlabel(時(shí)間/ms');title('時(shí)域');[k]=fft(audio0,Fs);subplot(2,1,2);plot(abs(k));plot(20*log10(abs(abs(k))/max(abs(abs(k)))));ylabel('分貝/dB');xlabel('頻率/Hz');title('頻域');figure;subplot(2,1,1);plot(abs(k));ylabel('振幅/A');xlabel('頻率/Hz');set(gca,'XLim',[01000]);set(gca,'XTick',[0:20:1000]);fp=1000;fs=900;wp=2*fp*pi/Fs;ws=2*fs*pi/Fs;wc=(wp+ws)/2/piwdp=wp-ws;N=ceil(8*pi/wdp)N=N+mod(N,2);HPfir=fir1(N,wc,'high',hanning(N+1));figure;plot(HPfir);title('濾波器時(shí)域');subplot(2,1,2);plot(abs(fft(HPfir)));title('濾波器頻域');audio1=conv(audio0,HPfir);audio1=filter(HPfir,1,audio0)p1=audioplayer(audio1,Fs);figure;subplot(2,1,1);plot(audio1);title('濾波后時(shí)域');[k0]=fft(audio1,Fs);subplot(2,1,2);plot(abs(k0));title('濾波后頻域');xlabel('f(Hz)');audio1=audio1*10;p1=audioplayer(audio1,Fs);wavwrite(audio1,Fs,nbits,'C:\Users\Administrator\Desktop\signal\'五、結(jié)論我們將包含有幾個(gè)諧波分量的蜂鳴信號(hào)干擾了的人聲信號(hào)讀入到MATLAB當(dāng)中，對(duì)讀取的音頻信號(hào)做出其時(shí)