語音訊號濾波去噪——使用FLATTOPWIN設計的FIR濾波器

阿新 • • 發佈：2019-02-04

摘要 本課程設計主要內容是設計利用視窗設計法選擇FLATTOPWIN窗設計一個FIR
濾波器，對一段含噪語音訊號進行濾波去噪處理並根據濾波前後的波形和頻譜分析濾波效能。本課程設計模擬平臺為MATLAB7.0，開發工具是M語言程式設計，通過課程設計瞭解FIR濾波器設計的原理和步驟，掌握用MATLAB語言設計濾波器的方法，瞭解FLATTOPWIN對FIR濾波器的設計及程式設計方法。首先利用windows自帶的錄音機錄製一段語音訊號，加入一單頻噪聲，對訊號進行頻譜分析以確定所加噪聲頻率，設計濾波器進行濾波去噪處理，比較濾波前後的波形和頻譜並進行分析。由分析結果可知，濾波
後的語音訊號與原始訊號基本一致，即設計的FIR濾波器能夠去除訊號中所加單頻噪聲，達到了設計目的。
關鍵詞

濾波去噪；FIR濾波器；FLATTOPWIN窗；MATLAB

引言

本課程設計主要解決在含噪情況下對語音訊號的濾波去噪處理，處理時採用的是利用視窗設計法選擇FLATTOPWIN窗設計的FIR濾波器[1]。通過對比濾波前後波形圖和濾波前後語音訊號的對比，可以看出濾波器對有用訊號無失真放大具有重大意義。

課程設計目的

熟悉Matlab語言環境，掌握Matlab語言的程式設計規則，利用Flattopwin窗函式設計法來設計符合要求的FIR濾波器來實現語音訊號的濾波去噪。並繪製濾波前後的時域波形和頻譜圖。根據圖形分析判斷濾波器設計的正確性。通過本次課程設計熟悉利用
flattopwin窗函式法設計FIR濾波器的過程[2]。增強自己獨立解決問題的能力，提高自己的動手能力。加深對理論知識聯絡實際問題的理解。為以後的工作奠定堅實的基礎。

課程設計要求

錄製一段語音，繪製波形並觀察其頻譜特點，加入一個帶外單頻噪聲，設計一個滿足指標的濾波器，對該含噪語音訊號進行濾波去噪處理，比較濾波前後的波形和頻譜並進行分析，根據結果和學過的理論得出合理的結論[3]。

課程設計平臺

MATLB是美國MathWorks公司出品的商業教學軟體，用於演算法開發。資料視覺化、資料分析以及數值計算的高階技術計算語言和互動式環境，主要包括MATLAB和Simulink兩大部分。
MATLAB由一系列工具組成。這些工具方便使用者使用MATLAB的函式和檔案，其中許多工具採用的是圖形使用者介面。包括MATLAB桌面和命令視窗、歷史命令視窗、編輯器和偵錯程式、路徑搜尋和用於使用者瀏覽幫助、工作空間、檔案的瀏覽器。隨著MATLAB的商業化以及軟體本身的不斷升級，MATLAB的使用者介面也越來越精緻，更加接近Windows的標準介面，人機互動性更強，操作更簡單。而且新版本的MATLAB提供了完整的聯機查詢、幫助系統，極大的方便了使用者的使用。簡單的程式設計環境提供了比較完備的除錯系統，程式不必經過編譯就可以直接執行，而且能夠及時地報告出現的錯誤及進行出錯原因分析。
MATLAB是一個高階的矩陣/陣列語言，它包含控制語句、函式、資料結構、輸入和輸出和麵向物件程式設計特點。使用者可以在命令視窗中將輸入語句與執行命令同步，也可以先編寫好一個較大的複雜的應用程式（M檔案）後再一起執行。新版本的MATLAB語言是基於最為流行的C++語言基礎上的，因此語法特徵與C++語言極為相似，而且更加簡單，更加符合科技人員對數學表示式的書寫格式。使之更利於非計算機專業的科技人員使用。而且這種語言可移植性好、可拓展性極強，這也是MATLAB能夠深入到科學研究及工程計算各個領域的重要原因[4]。
MATLAB是一個包含大量計算演算法的集合。其擁有600多個工程中要用到的數學運算函式，可以方便的實現使用者所需的各種計算功能。函式中所使用的演算法都是科研和工程
計算中的最新研究成果，而且經過了各種優化和容錯處理。在通常情況下，可以用它來代替底層程式語言，如C和C++ 。在計算要求相同的情況下，使用MATLAB的程式設計工作量會大大減少。MATLAB的這些函式集包括從最簡單最基本的函式到諸如矩陣，特徵向量、快速傅立葉變換的複雜函式。函式所能解決的問題其大致包括矩陣運算和線性方程組的求解、微分方程及偏微分方程的組的求解、符號運算、傅立葉變換和資料的統計分析、工程中的優化問題、稀疏矩陣運算、複數的各種運算、三角函式和其他初等數學運算、多維陣列操作以及建模動態模擬等。

基本理論

FIR濾波器

FIR(Finite Impulse Response)濾波器：有限長單位衝激響應濾波器，又稱為非遞迴型濾波器，是數字訊號處理系統中最基本的元件，它可以在保證任意幅頻特性的同時具有嚴格的線性相頻特性，同時其單位抽樣響應是有限長的，因而濾波器是穩定的系統。因此，FIR濾波器在通訊、影象處理、模式識別等領域都有著廣泛的應用。
工作原理：在進入FIR濾波器前，首先要將訊號通過A/D器件進行模數轉換，把模擬訊號轉化為數字訊號；為了使訊號處理能夠不發生失真，訊號的取樣速度必須滿足奈奎斯特定理，一般取訊號頻率上限的4-5倍做為取樣頻率；一般可用速度較高的逐次逼進式A/D轉換器，不論採用乘累加方法還是分散式演算法設計FIR濾波器，濾波器輸出的資料都是一串序列，要使它能直觀地反應出來，還需經過數模轉換，因此由FPGA構成的FIR濾波器的輸出須外接D/A模組。FPGA有著規整的內部邏輯陣列和豐富的連線資源，特別適合於數字訊號處理任務，相對於序列運算為主導的通用DSP晶片來說，其並行性和可擴充套件性更好，利用FPGA乘累加的快速演算法，可以設計出高速的FIR數字濾波器[5]。工作原理圖如圖1

圖一，FIR濾波器工作原理圖

視窗設計法

視窗設計法是一種通過截斷和計權的方法使無限長非因果序列成為有限長脈衝響應序列的設計方法。通常在設計濾波器之前，應該先根據具體的工程應用確定濾波器的技術指標。在大多數實際應用中，數字濾波器常常被用來實現選頻操作，所以指標的形式一般為在頻域中以分貝值給出的相對幅度響應和相位響應。
視窗設計法步驟如下：
（1）根據過渡頻寬及阻帶衰減要求，選擇窗函式的型別並估計視窗長度N。窗函式的型別可根據最小阻帶衰減AS獨立選擇[6]。
（2）根據待求濾波器的理想頻率響應求出理想單位脈衝響應hd(n)。
（3）由效能指標確定窗函式W(n)和長度N。
（4）求得實際濾波器的單位脈衝響應h(n)， h(n)即為所設計FIR濾波器係數向量b(n)。
這裡寫圖片描述（2.1）

常見的窗函式效能表如下表2.1所示：
表2.1 常見窗函式效能表

FLATTOPWIN窗

w=Flattopwin (L) 返回L-點Flattopwin視窗中列向量。Flattopwin窗的濾波器的過渡頻寬為19.6π/M，最小阻帶衰減108db。
這裡寫圖片描述
時間波形和幅度譜如下圖2.2、圖2.3：
圖2.2 時間波形
圖2.3 幅度譜

設計步驟

設計流程圖

根據設計的要求，首先自己錄製一段語音訊號，修改語音檔案格式，對語音訊號加入噪聲干擾，再利用Flattopwin窗設計合理的FIR濾波器。最後用濾波器對干擾後的語音訊號進行濾波去噪。具體設計流程圖如下圖3.1所示：
![這裡寫圖片描述]

圖3.1設計流程圖

錄製語音訊號

從電腦上錄製一段語音訊號，並命名為“cf.wav”，修改語音檔案的格式，並放在E盤目錄下。在MATLAB軟體中呼叫wavread函式可採集到語音訊號。程式碼如下：

[x,fs,bits]=wavread('e:\cf.wav');

%fs是生成該波形檔案的取樣頻率，bits是波形檔案沒樣本的編碼位數

得到原始語音訊號時域波形圖如圖3.2
這裡寫圖片描述圖3.2 原始語音訊號時域波形圖

然後對語音號進行快速傅立葉變換，得到訊號的頻譜特性，並將原始音樂訊號的波形圖與加干擾後的波形圖進行比較。部分程式碼如下：

y=x+0.1*sin(fn*2*pi*t);    %給原始訊號加噪聲
X=abs(fft(x));             %對原始訊號進行傅立葉變換
Y=abs(fft(y));             %對加噪聲後訊號進行傅立葉變換

加噪聲前後的時域與頻譜圖如圖3.3
這裡寫圖片描述
圖3.3 加噪聲前後的時域圖與頻譜圖比較

濾波器設計

在本次的課程設計中我所採用的是利用Flattopwin窗函式來設計FIR濾波器。其中主要用了freqz_m.m和ideal_lp.m兩個自編函式。
設計的濾波器圖如圖3.4
這裡寫圖片描述
圖3.4 濾波器圖形

訊號濾波處理

濾波器設計完成後，在MATLAB平臺上用函式fftfilt實現濾波。
得到的濾波前後語音訊號的時域波形圖和頻譜圖對比圖如圖3.5、3.6
這裡寫圖片描述
圖3.5 濾波前後語音訊號的時域波形圖和頻譜圖

圖3.6 濾波前後語音訊號的比較

結果分析

在MATLAB中，對原始的語音訊號加噪音，利用FLATTOPWIN視窗設計FIR濾波器，通過所設計的濾波器對加噪聲後的訊號進行處理。再通過sound(x,fs,bits)函式對濾波後的語音訊號進行回放，可以聽到濾波之後的訊號和原始訊號一樣清晰。具體程式碼如下：

sound(x,fs,bits);        %播放原始語音訊號
sound (y_fil,fs,bits);   % 播放濾波後的音樂訊號

所得結果證明了用Flattopwin窗設計的FIR濾波器和音樂訊號去噪設計是成功的。

出現的問題及解決方法

在本次課程設計中我遇到的問題如下：
1、對Flattopwin窗不是很瞭解，對用Flattopwin窗函式設計FIR濾波器的設計步驟很生疏。
2、語音檔案的格式有問題，不知道如何修改。
3、不知道如何呼叫函式，對MATLAB的使用不熟悉。
4、對濾波器一些引數以及通帶阻帶等概念不是很清楚。
5、在採用Flattopwin窗函式設計的FIR濾波器時得不到理想的濾波器，因而訊號的恢
復不是特別理想。

針對以上問題，相應的解決方案如下：

1、自己上網查閱資料，以及檢視之前《數字訊號處理》教材和向圖書館借閱資料，掌
握利用Flattopwin窗函式設計FIR濾波器的方法和步驟。
2、詢問同學，找同學幫忙。
3、利用MATLAB這款軟體自帶的help功能，以及檢視老師所給提供的課設資料和上
網看閱資料。
4、請教值班老師以及在網上查閱資料。
5、通過不斷設定引數的值，最終達到最理想的值，設計出理想的濾波器，使訊號得到理想恢復。

結束語

本次的課程設計，我的題目是《語音訊號濾波去噪——使用Flattopwin設計的FIR濾波》。在本次課程設計之前，我對Flattopwin窗函式完全沒有了解，因此在看到這個題目時，我是一頭霧水。但是通過自己翻閱資料和詢問同學，我掌握了用Flattopwin窗函式設計FIR濾波器的方法步驟，瞭解了窗函式的基本設計流程。經過這三週的課程設計，我學會了很多東西，不僅僅是技術，更多的是自學能力的提高以及處理問題的方法。
我們通訊工程專業是個實踐性很強的專業，而我們在校大部分的學習時間都是花在理論學習上面，實踐的機會很少。因而我對很多所學的理論知識如何跟實踐聯絡的概念很模糊，這次的課程設計給了我這個機會，加深了我對理論聯絡實際的理解，增強了自己獨立分析問題和解決問題的能力，開闊了自己的思維。
還有讓我看到了自己的不足，自己對本專業的相關知識掌握的還很少，還有很多知識都沒掌握，還讓我認識到解決問題的方法、途徑很多，做事要開闊自己的思維，看待問題要從多個角度看。
在此我要感謝學校為我們提供這次課程設計的機會，感謝老師對我的悉心指導，也感謝同學對我的幫助。這次的課程設計讓我理論聯絡實際，不僅鞏固了我們的理論知識，還提高了我的動手能力，在這次課程設計中我所學到的知識是我的財富，讓我終身受益。

程式碼

語音訊號濾波去噪設計源程式清單

% 程式名稱：cxf.m
% 程式功能：利用FLATTOPWIN設計的FIR濾波對語音訊號進行濾波去噪
% 程式作者： 曹雪鋒
% 最後修改日期：2017年3月8日
[x,fs,bits]=wavread('e:\cf.wav');  % 輸了引數為檔案的全路徑和檔名，輸出的第一個引數是每個樣本的值，fs是生成該波形檔案時的取樣率，bits是波形檔案每樣本的編碼位數。
sound(x,fs,bits);      % 按指定的取樣率和每樣本編碼位數回放
N=length(x);       % 計算訊號X的長度
fn=2200;            % 單頻噪聲頻率
t=0:1/fs:(N-1)/fs;         % 計算時間範圍，樣本數除以取樣頻率
x=x(:,1)';               % 將雙聲道轉為單聲道
y=x+0.1*sin(fn*2*pi*t);     % 加噪聲
sound(y,fs,bits);         % 應該可以明顯聽出有尖銳的單品嘯叫聲
X=abs(fft(x));  Y=abs(fft(y));    % 對原始訊號和加噪聲訊號進行fft變換，取幅度值
X=X(1:N/2); Y=Y(1:N/2);       % 擷取前半部分
deltaf=fs/N;            % 計算頻譜的譜線間隔
f=0:deltaf:fs/2-deltaf;     % 計算頻譜頻率範圍
figure(1)
subplot(2,2,1);plot(t,x);xlabel('時間t');
ylabel('幅度');title('原始語音訊號');
axis([0,4,-1.5,1.5]);

subplot(2,2,2);plot(f,X);xlabel('頻率f');ylabel('幅度譜');
title('原始語音訊號幅度譜');
axis([-10,6000,0,700]);
subplot(2,2,3);plot(t,y);xlabel('時間');ylabel('幅度');
title('加干擾後的語音訊號');
axis([0,4,-1.5,1.5]);
subplot(2,2,4);plot(f,Y);xlabel('頻率 f');ylabel('幅度譜');
title('加干擾後的語音訊號幅度譜');
axis([-10,6000,0,700]);

fpd=2100;fsd=2150;fsu=2250;fpu=2300;Rp=1;As=30;      % 阻帶濾波器設計指標
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;    %將Hz為單位的模擬頻率換算為rad為單位的數字頻率 
wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;
M=ceil(6.2*pi/dw)+1;      %計算flattopwin窗設計濾波器時需要的階數
n=0:M-1;         % 定義時間範圍
w_fla=(Flattopwin(M));    % 產生M階的flattopwin窗
hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); 
% 呼叫自編函式計算理想阻帶濾波器的脈衝響應
h_bs=w_fla'.*hd_bs;       % 用視窗法計算實際濾波器脈衝響應
[db,mag,pha,grd,w]=freqz_m(h_bs,1);  % 呼叫自編函式計算濾波器的頻率特徵

figure(2)
subplot(2,2,1);plot(w,db);title('濾波器幅度響應圖');
xlabel('w/pi');ylabel('db');
axis([0,0.8,-60,10]);
line([0,2],[-As,-As],'color','y','linestyle','--','LineWidth',2);
line([0,2],[-Rp,-Rp],'color','g','linestyle','--','LineWidth',2);
line([wsd,wsd],[-30,10],'color','r','linestyle','--','LineWidth',2);
line([wsu,wsu],[-30,10],'color','r','linestyle','--','LineWidth',2);
subplot(2,2,2);plot(w,mag);title('濾波器幅度響應圖');
xlabel('w/pi');ylabel('幅度mag');
axis([0,1,-0.5,1.5]);
subplot(2,2,3);plot(w,pha);title('濾波器幅度響應圖');
xlabel('w/pi');ylabel('相位pha');
axis([0,1,-4,4]);
subplot(2,2,4);stem(n,h_bs);title('濾波器幅度響應圖');
xlabel('w/pi');ylabel('db');
axis([0,1500,0,1]);

y_fil=fftfilt(h_bs,y);  %用設計好的濾波器對y進行濾波
Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2);  %計算頻譜取前一半
figure(3)
subplot(3,2,1);plot(t,x);xlabel('時間t');ylabel('幅度');


title('原始語音訊號');
subplot(3,2,2);plot(f,X);xlabel('頻率f');ylabel('幅度譜');
title('原始語音訊號幅度譜');
axis([0,4000,0,600]);
subplot(3,2,3);plot(t,y);xlabel('時間t');ylabel('·幅度');
title('加干擾後語音訊號');
subplot(3,2,4);plot(f,Y);xlabel('頻率f');ylabel('幅度譜');
title('加干擾後語音訊號幅度譜');
axis([0,4000,0,600]);
subplot(3,2,5);plot(t,y_fil);xlabel('時間t');ylabel('幅度');
title('濾波後語音訊號');
subplot(3,2,6);plot(f,Y_fil);xlabel('頻率f');ylabel('幅度譜');
title('濾波後語音訊號幅度譜');
axis([0,4000,0,600]);
figure(4)
subplot(2,1,1);plot(f,20*log10(Y));grid on; 
subplot(2,1,2);plot(f,20*log10(Y_fil));grid on; 
sound(x,fs,bits);   % 播放原始訊號
sound (y_fil,fs,bits);  %播放濾波後的訊號

[C,B,A]=dir2fs(h_bs);

函式FREQZ_M.M定義：

function [db,mag,pha,grd,w] = freqz_m(b,a);
% Modified version of freqz subroutine
% ------------------------------------
% [db,mag,pha,grd,w] = freqz_m(b,a);
%  db = Relative magnitude in dB computed over 0 to pi radians
% mag = absolute magnitude computed over 0 to pi radians 
% pha = Phase response in radians over 0 to pi radians
% grd = Group delay over 0 to pi radians
%   w = 501 frequency samples between 0 to pi radians
%   b = numerator polynomial of H(z)   (for FIR: b=h)
%   a = denominator polynomial of H(z) (for FIR: a=[1])
%
[H,w] = freqz(b,a,1000,'whole');
    H = (H(1:1:501))'; w = (w(1:1:501))';
  mag = abs(H);
   db = 20*log10((mag+eps)/max(mag));
  pha = angle(H);
%  pha = unwrap(angle(H));
  grd = grpdelay(b,a,w);
%  grd = diff(pha);
%  grd = [grd(1) grd];
%  grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
%  grd = median(grd)*500/pi;

函式IDEAL_LP.M定義：

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);

該文為筆者課程設計論文，如有不妥，還請指出