The matched filter is one of the most useful tools available to the Digital Signal Processing (DSP) engineer. Matched filtering maximizes the Signal to Noise Ratio (SNR) of the filter output when a signal matching a known pattern is passed through it. This process is extremely useful for “finding” signal patterns, possibly hidden within a noisy signal (low SNR). One application is in radar applications where the received signal may be orders of magnitude lower in amplitude than the transmitted waveform. In such case, a matched filter may be used to improve the SNR to a level which allows for proper detection of the target.
Without delving into the fundamental details of signal processing, suffice to say that the core of the matched filter can be constructed from a time-reversed (and conjugated if the signal is complex) version of a known signal. After construction of the the core, the matched filter is then defined by the convolution of the measured signal and the core.
For practical purposes, it is convenient to remember that convolution in the time domain is multiplication in the frequency domain. As a result, we can perform a Fourier Transform on both parts of the equation on the right hand side to obtain the frequency domain components of each. Following that, we multiply them, then perform the Inverse Fourier Transform to return the time domain result. In DSP, we use the Discrete Fourier Transform (DFT) and its inverse to complete the calculations. An efficient implementation of the DFT is the Fast Fourier Transform (FFT), which is often constructed using the radix-2 Cooley–Tukey algorithm. The frequency domain method is usually much faster than time-domain convolution.
It is important to remember that the matched filter output does not preserve the shape of the original signal in the time domain, but will simply maximize the SNR where the known signal and the filter input overlap. When a “matching” input signal is introduced, the shape of the matched filter output will correspond to the autocorrelation of the known signal delayed to have a maximum peak at time t0. However, in the frequency domain, the filter output is simply a scaled version of the known signal.
With that background, I constructed an example of matched filtering in Matlab to illustrate the utility of this DSP tool. In this case, I chose to use a Linear Frequency Modulated (LFM) signal as an example. An LFM signal sweeps the frequency for a specified period of time and bandwidth. LFM waveforms are used extensively in radar systems, especially where the peak transmitted power is limited by hardware and smaller power levels are transmitted for longer times to increase the average power of a pulse. One type of radar that uses LFM waveforms are Frequency Modulated Continuous Wave (FMCW) radars. In FMCW systems, the matched filter output resembles a sinc function with range resolution proportional to the bandwidth of the transmitted LFM pulse.
The Matlab code to generate the LFM waveform is listed below. In this case, I chose a complex 100 μs pulse with a swept bandwidth of 1 MHz, centered at DC (-0.5 MHz to 0.5 MHz). It should be noted that I chose to represent the analysis using quadrature signals (In-phase and Quadrature, IQ). There is an excellent description of quadrature signals by Dr. Lyons at https://www.dsprelated.com/showarticle/192.php
Fs = 5e6;
L = 100e-6;
t = 0:1/Fs:L-1/Fs;
A = 1;
% FM chirp waveform
f0 = -500e3;
f1 = 500e3;
c = (f1 - f0)/L;
waveform = A*exp(1j*2*pi*((c/2)*t.^2 + f0*t));
The LFM pulse varies in frequency linearly from the beginning of the pulse to the end, with the power spectral density plot showing the expected flat frequency domain response over 1 MHz bandwidth.
Autocorrelation, in a discrete time sense, is the correlation of a signal with a delayed copy of itself as a function of sample delay. The autocorrelation of the LFM pulse is interesting since it shows the characteristics of a compressed pulse. This reveals what the output of our matched filter will look like.
Remember that our LFM pulse is 100 μs long, but the autocorrelation response is a sinc function, which is compressed in time and adds gain over the original signal. This is known as pulse compression, with the -3 dB width of the response given by:
In our case, T’ is 1 μs. The compression ratio is then T/T’ which is 100. The gain associated with the matched filter in this case is also significant. After pulse compression, the power of the received signal can be considered as being amplified by T∆f, with a gain of approximately 20 dB in our case. In a future blog post, I will show how this applies to FMCW radar, specifically how bandwidth affects range resolution and filter gain.
For this case, I chose a signal of 1 ms duration and placed the LFM waveform at a random location within the signal. Noise was then added to obscure the waveform. Notice that the signal is almost indistinguishable from the noise in the image on the right.
The following code was used to construct the matched filter coefficients:
% Matched filter
% Time reversal and conjugate (since complex)
mf = conj(fliplr(waveform));The coefficients are a time-reversed and conjugated version of the original waveform. Both the original waveform and the filter coefficients are shown below for comparison.
In order to perform the matched filtering, I zero-padded the matched filter coefficients to match the length of the signal input. The FFT is performed on both the filter coefficients and the signal, then the results are multiplied together, after which the inverse FFT is performed. Here is the code:
% The matched filter is a time domain convolution of the time-reversed
% original waveform and the signal of interest. Since convolution in the
% time domain is multiplication in the frequency domain, we can take the
% product of the DFTs of both and then perform the inverse DFT to get the
% desired output. First, we need to pad the matched filter with zeros in
% order that the size of the signal and filter are the same
f = fft([mf, zeros(1,length(s)-length(mf))]);
p = fft(s);
res1 = ifft(p.*f);The resulting matched filter output clearly shows the embedded signal that was obscured within the noise.
From the plot on the left, we can identify the time compressed result with the shape of the autocorrelation function, which is a sinc response. The image on the right shows the result has a 3 dB bandwidth of approximately 1 μs and a SNR gain of approximately 20 dB, which is what we expected from the analysis above.
There is an additional aspect that should be explored in our matched filter example. In order to reduce the sidelobes associated with the matched filter result, we can scale the filter coefficients by a window function. Windowing the data also reduces an effect known as spectral leakage in the DFT. The end result is a better performing filter, with the downsides of lower pulse compression gain and a wider (3 dB bandwidth) filter output. This could have some implications in radar applications regarding Constant False Alarm Rate (CFAR) and range resolution, but that is out of scope for this blog post. A good article regarding windowing can be found here https://download.ni.com/evaluation/pxi/Understanding%20FFTs%20and%20Windowing.pdf
In that article, it gives some guidelines regarding the window choice:
- If the signal contains strong interfering frequency components distant from the frequency of interest, choose a smoothing window with a high side lobe roll-off rate.
- If the signal contains strong interfering signals near the frequency of interest, choose a window function with a low maximum side lobe level.
- If the frequency of interest contains two or more signals very near to each other, spectral resolution is important. In this case, it is best to choose a smoothing window with a very narrow main lobe.
- If the amplitude accuracy of a single frequency component is more important than the exact location of the component in a given frequency bin, choose a window with a wide main lobe.
- If the signal spectrum is rather flat or broadband in frequency content, use the uniform window, or no window.
- In general, the Hanning (Hann) window is satisfactory in 95 percent of cases. It has good frequency resolution and reduced spectral leakage. If you do not know the nature of the signal but you want to apply a smoothing window, start with the Hann window.
I chose a Taylor window for this example. I tried others as well, but the performance was similar for this particular instance. The window function is the image on the left and the windowed matched filter coefficients are on the right.
The autocorrelation of the windowed waveform reveals the reduced sidelobe level and wider 3 dB bandwidth of the expected matched filter response.
The block diagram for the matched filter shows that some computational efficiency can be gained by calculating the DFT of the product of the matched filter coefficients and the window coefficients (zero-padded to match the length of the input signal) ahead of time and stored in memory. This eliminates one of the DFT calculations for each filter execution. However, if the length of the input signal length changes, the filter DFT will need to be recalculated and stored.
The Matlab code to perform this is listed below:
% The matched filter processing with Taylor window
w1 = taylorwindow(Fs*L,15,-65);
f = fft([mf.*w1, zeros(1,length(s)-length(mf))]);
p = fft(s);
res2 = ifft(p.*f);The resulting matched filter output of the windowed matched filter clearly shows the embedded signal that was obscured within the noise. With less noise present, the reduced sidelobes would be more apparent.
When the non-windowed and windowed matched filter outputs are normalized and compared, the lower sidelobes can be seen, albeit with a slight reduction in power and wider 3 dB bandwidth.
There are many sources of information regarding the matched filter, and some of my preferred references are:
G, L. R. (2011). Understanding Digital Signal Processing.
Budge, M. C., & German, S. R. (2020). Basic Radar Analysis, Second Edition.
Jankiraman, M. (2018). FMCW Radar Design.
The complete Matlab code for this example:
clear all
close all
% Random number generator seed
rng(2)
% Directory to store images
DirectoryPath = uigetdir();
Fs = 5e6;
L = 100e-6;
t = 0:1/Fs:L-1/Fs;
A = 1;
% FM chirp waveform
f0 = -500e3;
f1 = 500e3;
c = (f1 - f0)/L;
waveform = A*exp(1j*2*pi*((c/2)*t.^2 + f0*t));
figure
plot(t,real(waveform),t,imag(waveform),'Linewidth',1)
grid on
title('Symmetric LFM Waveform')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure1.png']);
saveas(gcf, whereToStore);
% Power Spectral Density
n = 2^nextpow2(length(waveform));
y = fftshift(fft(waveform,n));
% scale
psdx = (1/(Fs*n)) * abs(y).^2;
f = -Fs/2:Fs/n:Fs/2-Fs/n;
figure
plot(f/1e6,10*log10(psdx/max(psdx)),'Linewidth',1)
grid on
xlim([(-Fs/2)/1e6 (Fs/2)/1e6])
ylim([-50 5])
title('Normalized Power Spectral Density')
ylabel('Power/Frequency (dB/Hz)')
xlabel('Frequency (MHz)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure2.png']);
saveas(gcf, whereToStore);
% Normalized Autocorrelation
[ac,lags] = xcorr(waveform,'normalized');
figure
plot(lags/Fs,real(ac),lags/Fs,imag(ac))
grid on
title('Normalized Autocorrelation of LFM Waveform')
xlabel('Lag (s)')
ylabel('Autocorrelation')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure3.png']);
saveas(gcf, whereToStore);
% Construct a time series with a signal and noise
t = 0:1/Fs:L*10-1/Fs;
s = zeros(1,length(t));
% Pick a random spot for the signal
randpos = randi(length(s)-length(waveform),1);
% Insert the waveform
s(randpos:randpos+length(waveform)-1) = 1e-5*waveform;
figure
plot(t,real(s),t,imag(s),'Linewidth',1)
grid on
ylim([-5e-5 5e-5])
title('Signal with Embedded Waveform')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure4.png']);
saveas(gcf, whereToStore);
% Add noise
noise = 1e-5*complex(randn(1,length(t)), randn(1,length(t)));
s = s + noise;
figure
plot(t,real(s),t,imag(s),'Linewidth',1)
grid on
ylim([-5e-5 5e-5])
title('Signal with Embedded Waveform and Noise')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure5.png']);
saveas(gcf, whereToStore);
% Matched filter
% Time reversal and conjugate (since complex)
mf = conj(fliplr(waveform));
figure
plot(0:1/Fs:L-1/Fs,real(mf),0:1/Fs:L-1/Fs,imag(mf),'Linewidth',1)
grid on
title('Matched Filter Coefficients')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure6.png']);
saveas(gcf, whereToStore);
% The matched filter is a time domain convolution of the time-reversed
% original waveform and the signal of interest. Since convolution in the
% time domain is multiplication in the frequency domain, we can take the
% product of the DFTs of both and then perform the inverse DFT to get the
% desired output. First, we need to pad the matched filter with zeros in
% order that the size of the signal and filter are the same
f = fft([mf, zeros(1,length(s)-length(mf))]);
p = fft(s);
res1 = ifft(p.*f);
figure
plot(t,real(res1),t,imag(res1),'Linewidth',1)
grid on
title('Matched Filter Result')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure7.png']);
saveas(gcf, whereToStore);
figure
plot(t,10*log10(abs(res1).^2),'Linewidth',1)
grid on
title('Matched Filter Result')
xlim([0 L*10])
ylim([-75 -40])
xlabel('time (s)')
ylabel('Power (dB)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure8.png']);
saveas(gcf, whereToStore);
% Now use a Taylor window for the matched filter
w1 = taylorwindow(Fs*L,15,-65);
figure
plot(w1,'Linewidth',1)
grid on
title('Taylor Window Function nbar = 15 SLL = -65')
xlabel('Samples')
ylabel('Amplitude')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure9.png']);
saveas(gcf, whereToStore);
figure
plot(0:1/Fs:L-1/Fs,real(mf.*w1),0:1/Fs:L-1/Fs,imag(mf.*w1),'Linewidth',1)
grid on
title('Matched Filter Coefficients (Taylor Window)')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure10.png']);
saveas(gcf, whereToStore);
% Autocorrelation of the windowed waveform
[ac,lags] = xcorr(waveform.*w1,'normalized');
figure
plot(lags/Fs,real(ac),lags/Fs,imag(ac))
grid on
title('Normalized Autocorrelation of LFM Waveform (Taylor Window)')
xlabel('Lag (s)')
ylabel('Autocorrelation')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure11.png']);
saveas(gcf, whereToStore);
% The matched filter processing with Taylor window
f = fft([mf.*w1, zeros(1,length(s)-length(mf))]);
p = fft(s);
res2 = ifft(p.*f);
figure
plot(t,real(res2),t,imag(res2),'Linewidth',1)
grid on
title('Matched Filter Result (Taylor Window)')
xlabel('time (s)')
ylabel('Amplitude')
legend('In-phase (I)','Quadrature (Q)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure12.png']);
saveas(gcf, whereToStore);
figure
plot(t,10*log10(abs(res2).^2),'Linewidth',1)
grid on
title('Matched Filter Result (Taylor Window)')
xlim([0 L*10])
%xlim([.45e-3 .55e-3])
ylim([-75 -40])
xlabel('time (s)')
ylabel('Power (dB)')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure13.png']);
saveas(gcf, whereToStore);
% Normalize to the unwindowed response and compare
res1norm = res1/max(res1);
res2norm = res2/max(res1);
figure
plot(t,10*log10(abs(res1norm).^2),t,10*log10(abs(res2norm).^2),'Linewidth',1)
grid on
title('Matched Filter Result Comparison (Taylor Window)')
%xlim([0 L*10])
xlim([.4726e-3 .5124e-3])
ylim([-30 0])
xlabel('time (s)')
ylabel('Normalized Power (dB)')
legend('Rectangular Window','Taylor Window')
set(gcf,'position',[100,100,800,600])
whereToStore=fullfile(DirectoryPath,['figure14.png']);
saveas(gcf, whereToStore);
The taylorwindow() function I used:
function w = taylorwindow(L, nbar, sll)
rho = 10^(-sll/20);
A = acosh(rho)/pi;
ss = nbar*nbar/(A*A + (nbar-0.5)*(nbar-0.5));
T = [];
for m = 1:(nbar-1)
T(m) = taylorcoef(m,nbar,A,ss);
end
M = (L-1)/2;
w = [];
for k = 0:M
n = k-M;
sum = 0;
for m = 1:(nbar-1)
sum = sum + T(m)*cos(2.0*n*pi*m/L);
end
w(k+1) = 1 + 2*sum;
w(L-k) = w(k+1);
end
function F = taylorcoef(m,nbar,A,ss)
% Compute the Taylor window coefficients
prod1 = 1;
for n = 1:(nbar-1)
prod1 = prod1*(1 - m*m/ss/(A*A + (n-0.5)*(n-0.5)));
end
prod2 = 1;
for n = 1:(nbar-1)
if (n ~= m)
prod2 = prod2*(1 - m*m/n/n);
end
end
F = -0.5*((-1)^m)*prod1/prod2;