自适应滤波算法介绍之匹配滤波器的基本原理和应用示例

以下文章来源于FPGA算法工程师，作者18线工程师

概述

自适应滤波理论在统计信号处理中占据非常重要的地位，在通信、控制、雷达等领域获得广泛应用。自适应滤波器的基本目标，是通过某种方式对参数θ(k)进行调整，使滤波器输出尽可能使得包含参考信号的某个特定目标函数达到最小化。

在实时信号处理中，我们希望滤波器在实现滤波、平滑或预测等任务中，能够跟踪和适应系统或环境的动态变化，这需要滤波器的参数可以随时间做简单的变化或更新，因为复杂的运算不符合实施快速处理的要求。

【自适应滤波算法】系列将介绍常用的匹配滤波器、维纳滤波器、Kalman滤波、LMS算法和RLS算法。参考书籍推荐《现代信号处理》和《自适应滤波算法与实现》。

我们从最基础的匹配滤波器开始。

1. 匹配滤波器的概念

在滤波器设计中，遵循两种最优设计准则：

使滤波器的输出达到最大的信噪比，称为匹配滤波器；

使输出滤波器的均方估计误差为最小，称为Wiener滤波器。

考虑接收信号模型:

线性连续时间滤波器

从接收机的角度看，接收信号y(t)是已知发送信号s(t)经过信道传播后得到，n(t)为信道传播过程中引入的零均值平稳噪声。在实际系统中，我们希望进入基带的接收信号yo(t)信号质量尽可能好，于是通过设计滤波器的冲击响应函数h(t)，使得滤波器输出的信噪比最大。

根据信号与系统的知识，我们可以得到滤波器的输出可以表示为：

滤波器在t＝To时刻的输出信噪比定义为：

根据Parseval定理：

则输出信号可以写作：

于是输出信号在t＝To的瞬时功率：

输出噪声的平均功率：

令Pn(w)是加性噪声n(t)的功率谱密度，则输出噪声的功率谱密度为：

于是，输出噪声的平均功率可以写成：

因此，得到信噪比：

再应用柯西-施瓦茨不等式：

得到：

将上式中等号成立时的滤波器传递函数记作Hopt(w),并利用柯西-施瓦茨不等式取等号的条件，则有：

于是，我们可以得到输出最大信噪比为：

在白噪声情况下，其功率谱密度Pn(w)＝1，故滤波器传递函数可以简化为：

因此有|Ho(w)|=|S*(w)=|S(w)|，当滤波器达到最大输出信噪比时，滤波器的幅频特性|H(w)|与信号s(t)的幅频特性|S(w)|相等。

通过对上式两边做IFFT，得到：

即：匹配滤波器的冲激响应ho(t)是信号s(t)的一镜像信号。

2. 匹配滤波器的性质

匹配滤波器在很多工程问题中有重要应用，其性质概括如下：

性质1：在所有线性滤波器中，匹配滤波器输出的信噪比最大，且SNRmax=2Es/(No)，与输入信号的波形以及加性噪声的分布特性无关。

性质2：匹配滤波器输出信号在t＝To时刻的瞬时功率达到最大。

性质3：匹配滤波器输出信噪比达到最大的时刻To，应该选取等于原信号的持续时间T。

性质4：匹配滤波器对波形相同而幅值不同的时延信号具有适应性。

性质5：匹配滤波器对频移信号不具有适应性。

局限性：在匹配滤波器中，接收机必须已知并存储信号的精确结构或功率谱，并且积分区间必须与信号取非零值的区间同步，但显然在实际系统中，原信号在传输过程中发生传播延迟，并伴随相偏和频偏，信号区间和积分区间无法精确同步而造成误差，匹配滤波难以实现。

3. 匹配滤波器示例

MATLAB示例： 在原点放置一个各向同性天线元件（0;0;0）。然后，在距离发射机约7公里处放置一个RCS（雷达截面积）为1平方米、不动的目标，位置为（5000;5000;10）。将工作（载波）频率设置为10 GHz。模拟单基地雷达的发送和接收,计算发射机到目标的距离和角度。

产生一个时长25us、脉冲频率为10 kHz的矩形脉冲波形。利用雷达方程确定探测目标所需的峰值功率。该目标在发射机工作频率和增益的最大不模糊距离下，RCS为1平方米。信噪比基于非相干检测期望的误检率10-6。

在接收端收集回波，并通过匹配滤波器提升信噪比。

matlab示例代码：

```

antenna = phased.IsotropicAntennaElement('FrequencyRange',[5e9 15e9]);

transmitter = phased.Transmitter('Gain',20,'InUseOutputPort',true);

fc = 10e9;

target = phased.RadarTarget('Model','Nonfluctuating',...

'MeanRCS',1,'OperatingFrequency',fc);

txloc = [0;0;0];

tgtloc = [5000;5000;10];

transmitterplatform = phased.Platform('InitialPosition',txloc);

targetplatform = phased.Platform('InitialPosition',tgtloc);

[tgtrng,tgtang] = rangeangle(targetplatform.InitialPosition,...

transmitterplatform.InitialPosition);

waveform = phased.RectangularWaveform('PulseWidth',25e-6,...

'OutputFormat','Pulses','PRF',10e3,'NumPulses',1);

c = physconst('LightSpeed');

maxrange = c/(2*waveform.PRF);

SNR = npwgnthresh(1e-6,1,'noncoherent');

lambda = physconst('LightSpeed')/target.OperatingFrequency;

Ts = 290;

dbterms = db2pow(SNR - 2*transmitter.Gain);

Pt = (4*pi)^3*physconst('Boltzmann')*Ts/waveform.PulseWidth/target.MeanRCS/(lambda^2)*maxrange^4*dbterms;

transmitter.PeakPower = Pt;

radiator = phased.Radiator('PropagationSpeed',c,...

'OperatingFrequency',fc,'Sensor',antenna);

channel = phased.FreeSpace('PropagationSpeed',c,...

'OperatingFrequency',fc,'TwoWayPropagation',false);

collector = phased.Collector('PropagationSpeed',c,...

'OperatingFrequency',fc,'Sensor',antenna);

receiver = phased.ReceiverPreamp('NoiseFigure',0,...

'EnableInputPort',true,'SeedSource','Property','Seed',2e3);

mfilter = phased.MatchedFilter(...

'Coefficients',getMatchedFilter(waveform),...

'GainOutputPort',true);

wf = waveform();

[wf,txstatus] = transmitter(wf);

wf = radiator(wf,tgtang);

wf = channel(wf,txloc,tgtloc,[0;0;0],[0;0;0]);

wf = target(wf);

wf = channel(wf,tgtloc,txloc,[0;0;0],[0;0;0]);

wf = collector(wf,tgtang);

rx_puls = receiver(wf,~txstatus);

[mf_puls,mfgain] = mfilter(rx_puls);

Gd = length(mfilter.Coefficients)-1;

mf_puls=[mf_puls(Gd+1:end); mf_puls(1:Gd)];

subplot(2,1,1)

t = unigrid(0,1e-6,1e-4,'[)');

rangegates = c.*t;

rangegates = rangegates/2;

plot(rangegates,abs(rx_puls))

title('Received Pulse')

ylabel('Amplitude')

hold on

plot([tgtrng, tgtrng], [0 max(abs(rx_puls))],'r')

subplot(2,1,2)

plot(rangegates,abs(mf_puls))

title('With Matched Filtering')

xlabel('Meters')

ylabel('Amplitude')

hold on

plot([tgtrng, tgtrng], [0 max(abs(mf_puls))],'r')

hold off

```

接收脉冲匹配滤波