基於MATLAB的FIR濾波器的設計

FIR濾波器設計的整體流程圖

本設計使用fir濾波器對語音信號進行濾波處理，本仿真設計使用matlab作為仿真平台，matlab自帶的信號作為語音原始信號來進行濾波器的仿真。其流程圖表示如下：

總設計流程圖

首先要設計的是fir濾波器，根據fir濾波器的理論形式，fir濾波器（有限長度沖擊響應）是全零點型的濾波器，其數學的實現形式如下：

\[y[n] = a_{0}x[n]+a_{1}[n-1]+...a_{m}x[n-m] \]

其中濾波器最重要的兩個特性為線性相位特性和幅度特性。原始信號的頻率特性曲線如圖1.1 所示，通過進行三次樣條插值處理，將語音信號的采樣率提升到20Khz，進行采樣后的頻譜圖如1.2所示。

圖1.1：原始信號的頻譜圖

圖1.2：提高采樣到20khz頻譜圖

觀察圖1.2后，可以發現信號的能量集中在低頻部分，為了減少無用的高頻分量，本設計設計了fir低通濾波將其濾除。本設計采樣了等紋波來完成低通濾波的設計。設計工具使用MATLAB的濾波器工具FADTOOL，FAD可以根據使用者的參數進行濾波器的快速驗證和設計。根據圖1.2的頻譜圖，確定設計的低通濾波器的fp = 8khz，fs = 8.5khz。設計的濾波器的頻譜圖為圖1.3所示，其階數為71階，單位圓的零極點分布如圖1.4所示.

圖1.3：低通濾波器的頻譜圖

圖1.4：FIR濾波器的零極點分布圖

語音信號通過低通濾波器后，其頻譜響應如圖1.5所示。

圖1.5 通過低通信號語音信號頻譜

可以明顯看出信號大於8khz以上的頻譜被過濾了。

過濾后使用希爾伯特濾波將語音信號的雙邊帶信號變換為單邊帶信號，減少通信帶寬是的使用。這里設計使用matlab的希爾伯特濾波\(\frac{1}{\pi t}\)過濾信號的頻譜如圖1.6所示，通過希爾伯特變化變成單邊帶頻率如圖1.7所示。

圖1.6 語音信號雙邊譜

圖1.7 語音信號單邊譜

由於信號需要進行發送，所有比如講頻率提高，將20khz的采樣率轉換到10M的采樣率。一般使用分階段的采樣將采樣率提高500倍。在每次的提高采樣率相當於在頻譜上進行延擴。為了保證頻譜的形狀，使用fir設置0.2和0.25的濾波器進行處理。圖1.8為進行升采樣5倍的頻譜圖。

圖1.8 語音信號升采樣5倍后的頻譜圖

所以通過截止頻率為0.2的濾波進行過濾，濾波器的頻率響應如圖1.9所示。

圖1.9 截止頻率為0.2的低通濾波器的頻譜圖 圖1.10為截止頻率為0.25的低通濾波器的頻譜圖


圖1.10截止頻率為0.2的低通濾波器的頻譜圖

當通過對語音信號進行500倍的upsample 后，其頻譜圖如下圖所示1.11所示。

圖1.11：upsample后的語音波形的頻譜

最后進行語音信號的頻率搬移，使得信號的可以發送。頻譜搬移后的語音信號頻譜圖如圖1.12所示。

圖1.12：語音信號通過頻譜搬移后的頻譜圖形

設計語音信號通過AWGN信道，相當在語音信號的頻譜上加上一個高斯白噪聲。高斯白噪聲的頻譜如圖1.12所示。

圖1.12 高斯白噪聲的頻譜圖

總結以上的過程，對語音信號進行采樣增加到10M，以滿足發送的條件，並使用高斯白噪聲模擬過程中信道造成的信號噪聲影響。這個時候，接收機要接收到原始的信號並且可以解調原始信號，必須有一個帶通濾波器，將語音信號過濾出來，而拋棄無用的噪聲信號。通過圖1.11的頻譜可以得到帶通濾波器的帶通為0.6-0.75，利用MATLAB的FAD工作選擇最小紋波設計，得到帶通濾波的歸一化系數矩陣，其頻率響應如圖1.13所示。

圖1.13：帶通濾波器的頻率響應圖形

將接受的信號通過帶通濾波后的語音頻譜圖如圖1.14所示。

圖1.14：語音信號通過帶通濾波器的頻率響應圖形

對接受到語音信號進行dowmsample（下邊頻處理），同樣下變頻500，過程和上變頻相反。最后得到聽到的語音信號的頻譜為圖1.15所示。

通過matlab的sound函數聽取原始信號和接受的信號，發現信號的語音基本一致，所以通過該方案進行語音信號的通信和濾波器的仿真是成功可行的。

clear all
close all
clc
%% 
load handel
h = fvtool(y);
set(h, 'Fs', Fs);
sound(y, Fs);
pause(10)

% 第一步，把語音的采樣率提升到20KHz
T = (length(y)-1) / Fs;
Y = interp1( [0:1/Fs:T], y, [0:1/20e3:T], 'spline');
h = fvtool(Y);
set(h, 'Fs', 20e3);
sound(Y, 20e3);
pause(10)


% 將信號8KHz以上部分濾掉,設計方法詳細參見fdatool
voiceCoef = [-0.00331567291687537,0.00302323989662157,0.00236335999433444,-0.000249895806141854,0.00324405325416548,-0.00268876267273492,0.00360296141490172,-0.00230710811907060,0.00100163740863659,0.00152582296780255,-0.00388560018913911,0.00604697828734485,-0.00697045946214204,0.00639302318173106,-0.00395271528961120,-2.33446800352772e-05,0.00492259451573453,-0.00964894184901816,0.0129924208661240,-0.0137833678820759,0.0112660249790563,-0.00531183793788302,-0.00339120178077228,0.0133543120529965,-0.0224465890588835,0.0282360533532409,-0.0284263495578446,0.0213310613052760,-0.00627693990578004,-0.0161429286985805,0.0440352410987295,-0.0744109723115057,0.103606516887248,-0.127851104964425,0.143880717742118,0.850515547752914,0.143880717742118,-0.127851104964425,0.103606516887248,-0.0744109723115057,0.0440352410987295,-0.0161429286985805,-0.00627693990578004,0.0213310613052760,-0.0284263495578446,0.0282360533532409,-0.0224465890588835,0.0133543120529965,-0.00339120178077228,-0.00531183793788302,0.0112660249790563,-0.0137833678820759,0.0129924208661240,-0.00964894184901816,0.00492259451573453,-2.33446800352772e-05,-0.00395271528961120,0.00639302318173106,-0.00697045946214204,0.00604697828734485,-0.00388560018913911,0.00152582296780255,0.00100163740863659,-0.00230710811907060,0.00360296141490172,-0.00268876267273492,0.00324405325416548,-0.000249895806141854,0.00236335999433444,0.00302323989662157,-0.00331567291687537];
dispFilterResponse(voiceCoef, 1/20e3, '');
voice = filter(voiceCoef, 1, Y);
h = fvtool(voice);
set(h, 'Fs', 20e3);
sound(voice, 20e3);
pause(10)

% 利用濾波法實現單邊帶,需要設計希爾伯特濾波器，1/(pi*t)
d = fdesign.hilbert('N,TW',200, 0.1);
Hd = design(d,'equiripple','SystemObject',true);
Realvoice = circshift(voice, 100);
Imagvoice = filter(Hd.Numerator, 1, voice);

h = fvtool(Realvoice);
set(h, 'FrequencyRange','[-pi, pi)');
set(h, 'Fs', 20e3);
h = fvtool(Realvoice + 1i*Imagvoice);
set(h, 'Fs', 20e3);


Basebandvoice = Realvoice + 1i*Imagvoice;


% 20K-》10M 需要提升500倍采樣率，這里分個階段 500 = 5 * 5 * 5 * 4，需要設計截止頻率分別是0.2和0.25的兩種濾波器
x5Coef = [0.00112692337568268,0.00158692025434125,0.00205896989050378,0.00197709384758320,0.00111473190427317,-0.000472447774337975,-0.00236605559284221,-0.00385845971455964,-0.00419282731353260,-0.00290235997397991,-0.000113082201514355,0.00334505346493630,0.00612323159188132,0.00682609929601474,0.00462600477070206,-0.000226559774923881,-0.00625173793228883,-0.0110992424858532,-0.0123738225850133,-0.00866846592173900,-0.000402088021268619,0.00993919626594012,0.0183840950231772,0.0208076890127266,0.0146172933092511,0.000185101945744734,-0.0185805496678055,-0.0348933246131477,-0.0409283295252545,-0.0303450353032694,-0.000682784674361512,0.0452255563608527,0.0994103424861294,0.150580752726224,0.187122188502142,0.200369871454618,0.187122188502142,0.150580752726224,0.0994103424861294,0.0452255563608527,-0.000682784674361512,-0.0303450353032694,-0.0409283295252545,-0.0348933246131477,-0.0185805496678055,0.000185101945744734,0.0146172933092511,0.0208076890127266,0.0183840950231772,0.00993919626594012,-0.000402088021268619,-0.00866846592173900,-0.0123738225850133,-0.0110992424858532,-0.00625173793228883,-0.000226559774923881,0.00462600477070206,0.00682609929601474,0.00612323159188132,0.00334505346493630,-0.000113082201514355,-0.00290235997397991,-0.00419282731353260,-0.00385845971455964,-0.00236605559284221,-0.000472447774337975,0.00111473190427317,0.00197709384758320,0.00205896989050378,0.00158692025434125,0.00112692337568268];
x4Coef = [0.000856994565694856,0.000591393177382970,-3.85351917650751e-05,-0.00130837344797337,-0.00269292916711479,-0.00331425686614336,-0.00244296037940879,-9.10921433942072e-05,0.00271899015730028,0.00432293869924751,0.00336191370561577,-0.000190182688008892,-0.00462501376328013,-0.00718415092444753,-0.00565091781262915,8.72027154765229e-06,0.00706737732439220,0.0111778371760113,0.00888277600993286,0.000135749254865209,-0.0108425521935169,-0.0173232948547436,-0.0138880550352155,-0.000292556888299184,0.0170984534611941,0.0277598485492297,0.0226701214314515,0.000423692149837658,-0.0297645353604724,-0.0505233937915195,-0.0438668534604438,-0.000511549520956627,0.0737481607287266,0.158297486795985,0.225157564979538,0.250542370102760,0.225157564979538,0.158297486795985,0.0737481607287266,-0.000511549520956627,-0.0438668534604438,-0.0505233937915195,-0.0297645353604724,0.000423692149837658,0.0226701214314515,0.0277598485492297,0.0170984534611941,-0.000292556888299184,-0.0138880550352155,-0.0173232948547436,-0.0108425521935169,0.000135749254865209,0.00888277600993286,0.0111778371760113,0.00706737732439220,8.72027154765229e-06,-0.00565091781262915,-0.00718415092444753,-0.00462501376328013,-0.000190182688008892,0.00336191370561577,0.00432293869924751,0.00271899015730028,-9.10921433942072e-05,-0.00244296037940879,-0.00331425686614336,-0.00269292916711479,-0.00130837344797337,-3.85351917650751e-05,0.000591393177382970,0.000856994565694856];

tmp = upsample(Basebandvoice, 5);
tmp = filter(x5Coef, 1, tmp);
tmp = upsample(tmp, 5);
tmp = filter(x5Coef, 1, tmp);
tmp = upsample(tmp, 5);
tmp = filter(x5Coef, 1, tmp);
tmp = upsample(tmp, 4);
tmp = filter(x4Coef, 1, tmp);
h = fvtool(tmp);
set(h, 'Fs', 10e6);

TX = tmp .* exp(1i * 2 * pi * 3.37e6 * [1:length(tmp)]/10e6);
TX = real(TX);
h = fvtool(TX);
set(h, 'Fs', 10e6);

%% 經過AWGN信道
RX = awgn(TX, -10, 'measured');
h = fvtool(RX);
set(h, 'Fs', 10e6);

%% 接收機
% 帶通濾波
bandpassCoef = [-0.000220985157474574	-0.000603199022228535	0.000567700620989989	-0.000144856835556256	-0.000537413993393477	0.000811952546260190	-0.000221093885213475	-0.000818227351198006	0.00123133214160066	-0.000365497144552332	-0.00113169408536944	0.00172960562049404	-0.000556315269703292	-0.00146815254434005	0.00229154024087189	-0.000793782309906573	-0.00180607660804564	0.00288989310575587	-0.00107147573190690	-0.00212115551648486	0.00349042252209103	-0.00137783322201238	-0.00238983933967575	0.00405551842218856	-0.00169658391251432	-0.00259081022572485	0.00454537432747044	-0.00200659683781057	-0.00270824168828992	0.00492493249491361	-0.00228582462373403	-0.00273356399437629	0.00516543626532525	-0.00251194689244855	-0.00266626182302852	0.00524743941405402	-0.00266626182302852	-0.00251194689244855	0.00516543626532525	-0.00273356399437629	-0.00228582462373403	0.00492493249491361	-0.00270824168828992	-0.00200659683781057	0.00454537432747044	-0.00259081022572485	-0.00169658391251432	0.00405551842218856	-0.00238983933967575	-0.00137783322201238	0.00349042252209103	-0.00212115551648486	-0.00107147573190690	0.00288989310575587	-0.00180607660804564	-0.000793782309906573	0.00229154024087189	-0.00146815254434005	-0.000556315269703292	0.00172960562049404	-0.00113169408536944	-0.000365497144552332	0.00123133214160066	-0.000818227351198006	-0.000221093885213475	0.000811952546260190	-0.000537413993393477	-0.000144856835556256	0.000567700620989989	-0.000603199022228535	-0.000220985157474574];
tmp = filter(bandpassCoef, 1, RX);
h = fvtool(tmp);
set(h, 'Fs', 10e6);

% 下變頻
tmp = tmp .* exp(-1i * 2 * pi * 3.37e6 * [1:length(tmp)]/10e6);
h = fvtool(tmp);
set(h, 'Fs', 10e6);  % 注意，這里存在一個鏡像，這個鏡像將由下邊的濾波去掉

% 低通濾波，下采樣到20KHz
tmp = filter(x4Coef, 1, tmp);
tmp = downsample(tmp, 4);
tmp = filter(x5Coef, 1, tmp);
tmp = downsample(tmp, 5);
tmp = filter(x5Coef, 1, tmp);
tmp = downsample(tmp, 5);
tmp = filter(x5Coef, 1, tmp);
tmp = downsample(tmp, 5);

h = fvtool(tmp);
set(h, 'Fs', 20e3);

% 播放音樂
sound(real(tmp)*1e3, 20e3);
pause(10)

%% 畫頻域相應的
function y = dispFilterResponse(coef, Fs, text) 
[h1,w1]=freqz(coef,1);  
plot(w1/pi*Fs/2,20*log10(abs(h1))); 
grid;
xlabel('頻率') ;
ylabel(text) ;
end