基于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