Practical Introduction to Digital Filtering

Compensating for Delay Introduced by Filtering

Digital filters introduce delay in your signal. Depending on the filter characteristics, the delay can be constant over all frequencies, or it can vary with frequency. The type of delay determines the actions you have to take to compensate for it. The grpdelay function allows you to look at the filter delay as a function of frequency. Looking at the output of this function allows you to identify if the delay of the filter is constant or if it varies with frequency (in other words, if it is frequency-dependent).

Filter delay that is constant over all frequencies can be easily compensated for by shifting the signal in time. FIR filters usually have constant delay. On the other hand, delay that varies with frequency causes phase distortion and can alter a signal waveform significantly. Compensating for frequency-dependent delay is not as trivial as for the constant delay case. IIR filters introduce frequency-dependent delay.

Compensating for Constant Filter Delay

As mentioned before, you can measure the group of delay of the filter to verify that it is a constant function of frequency. You can use the grpdelay function to measure the filter delay, D, and compensate for this delay by appending D zeros to the input signal and shifting the output signal in time by D samples.

Consider a noisy electrocardiogram signal that you want to filter to remove high frequency noise above 75 Hz. You want to apply an FIR lowpass filter and compensate for the filter delay so that the noisy and filtered signals are aligned correctly and can be plotted on top of each other for comparison.

Fs = 500;                    % Sample rate in Hz
N = 500;                     % Number of signal samples
rng default;
x = ecg(N)'+0.25*randn(N,1); % Noisy waveform
t = (0:N-1)/Fs;              % Time vector

Design a 70th-order lowpass FIR filter with a cutoff frequency of 75 Hz.

Fnorm = 75/(Fs/2);           % Normalized frequency
df3 = designfilt("lowpassfir",FilterOrder=70,CutoffFrequency=Fnorm);

Plot the group delay of the filter to verify that it is constant across all frequencies indicating that the filter is linear phase. Use the group delay to measure the delay of the filter.

grpdelay(df3,2048,Fs)         % Plot group delay

D = mean(grpdelay(df3)) % Filter delay in samples

D = 
35

Before filtering, append D zeros at the end of the input data vector, x. This ensures that all the useful samples are flushed out of the filter, and that the input signal and the delay-compensated output signal have the same length. Filter the data and compensate for the delay by shifting the output signal by D samples. This last step effectively removes the filter transient.

y = filter(df3,[x; zeros(D,1)]);  % Append D zeros to the input data
y = y(D+1:end);                  % Shift data to compensate for delay
plot(t,x,t,y,LineWidth=1.5)
title("Filtered Waveforms")
xlabel("Time (s)")
legend("Original Noisy Signal","Filtered Signal")
grid on
axis tight

Compensating for Frequency-Dependent Delay

Frequency-dependent delay causes phase distortion in the signal. Compensating for this type of delay is not as trivial as for the constant delay case. If your application allows off-line processing, you can remove the frequency-dependent delay by implementing zero-phase filtering using the filtfilt function. filtfilt performs zero-phase filtering by processing the input data in both the forward and reverse directions. The main effect is that you obtain zero-phase distortion, i.e., you filter data with an equivalent filter that has a constant delay of 0 samples. Other effects are that you get a filter transfer function which equals the squared magnitude of the original filter transfer function, and a filter order that is double the order of the original filter.

Consider the ECG signal defined in the previous section. Filter this signal with and without delay compensation. Design a 7th-order lowpass IIR elliptic filter with a cutoff frequency of 75 Hz.

Fnorm = 75/(Fs/2); % Normalized frequency
df4 = designfilt("lowpassiir", ...
    FilterOrder=7, ...
    PassbandFrequency=Fnorm, ...
    PassbandRipple=1, ...
    StopbandAttenuation=60);

Plot the group delay of the filter. The group delay varies with frequency, indicating that the filter delay is frequency-dependent.

grpdelay(df4,2048,Fs)

Filter the data and look at the effects of each filter implementation on the time signal. Zero-phase filtering effectively removes the filter delay.

y1 = filter(df4,x);    % Nonlinear phase filter - no delay compensation
y2 = filtfilt(df4,x);  % Zero-phase implementation - delay compensation
plot(t,x)
hold on
plot(t,y1,"r",LineWidth=1.5)
plot(t,y2,LineWidth=1.5)
title("Filtered Waveforms")
xlabel("Time (s)")
legend("Original Signal","Nonlinear Phase IIR Output", ...
  "Zero-Phase IIR Output")
xlim([0.25 0.55])
grid on

Zero-phase filtering is a great tool if your application allows for the non-causal forward/backward filtering operations, and for the change of the filter response to the square of the original response.

Filters that introduce constant delay are linear phase filters. Filters that introduce frequency-dependent delay are non-linear phase filters.

Removing Unwanted Spectral Content from a Signal

Filters are commonly used to remove unwanted spectral content from a signal. You can choose from a variety of filters to do this. You choose a lowpass filter when you want to remove high frequency content, or a highpass filter when you want to remove low frequency content. You can also choose a bandpass filter to remove low and high frequency content while leaving an intermediate band of frequencies intact. You choose a bandstop filter when you want to remove frequencies over a given band.

Consider an audio signal that has a power-line hum and white noise. The power-line hum is caused by a 60 Hz tone. White noise is a signal that exists across all the audio bandwidth.

Load the audio signal. Specify a sample rate of 44.1 kHz.

Fs = 44100;
y = audioread("noisymusic.wav");

Plot the power spectrum of the signal. The red triangular marker shows the strong 60 Hz tone interfering with the audio signal.

[P,F] = pwelch(y,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,[60 60],[-9.365 -9.365], ...
  {"Original Signal Power Spectrum", "60 Hz Tone"})

You can first remove as much white noise spectral content as possible using a lowpass filter. The passband of the filter should be set to a value that offers a good trade-off between noise reduction and audio degradation due to loss of high frequency content. Applying the lowpass filter before removing the 60 Hz hum is very convenient since you will be able to downsample the band-limited signal. The lower rate signal will allow you to design a sharper and narrower 60 Hz bandstop filter with a smaller filter order.

Design a lowpass filter with passband frequency of 1 kHz and stopband frequency of 1.4 kHz. Choose a minimum-order design.

Fp = 1e3;    % Passband frequency in Hz
Fst = 1.4e3; % Stopband frequency in Hz
Ap = 1;      % Passband ripple in dB
Ast = 95;    % Stopband attenuation in dB
df3 = designfilt("lowpassfir",PassbandFrequency=Fp, ...
    StopbandFrequency=Fst,PassbandRipple=Ap, ...
    StopbandAttenuation=Ast,SampleRate=Fs);

Analyze the filter response.

opts = filterAnalysisOptions;
opts.FrequencyRange = "freqvector";
opts.FrequencyVectorUnits = "Hz";
opts.FrequencyVector = F;
opts.FrequencyScale = "log";
filterAnalyzer(df3,AnalysisOptions=opts)

Filter the data and compensate for delay.

D = mean(grpdelay(df3)); % Filter delay
ylp = filter(df3,[y; zeros(D,1)]);
ylp = ylp(D+1:end);

Look at the spectrum of the lowpass filtered signal. The frequency content above 1400 Hz has been removed.

[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,Flp,Plp, ...
    {"Original Signal","Lowpass Filtered Signal"})

From the power spectrum plot above, you can see that the maximum non-negligible frequency content of the lowpass filtered signal is at 1400 Hz. By the sampling theorem, a sample rate of $2\times 1400=2800$ Hz would suffice to represent the signal correctly, you however, are using a sample rate of 44100 Hz which is a waste since you will need to process more samples than those necessary. You can downsample the signal to reduce the sample rate and reduce the computational load by reducing the number of samples that you need to process. A lower sample rate will also allow you to design a sharper and narrower bandstop filter, needed to remove the 60 Hz noise, with a smaller filter order.

Downsample the lowpass filtered signal by a factor of 10 to obtain a sample rate of Fs/10 = 4.41 kHz. Plot the spectrum of the signal before and after downsampling.

Fs = Fs/10;
yds = downsample(ylp,10);
[Pds,Fds] = pwelch(yds,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,Fds,Pds, ...
    {"Signal Sampled at 44100 Hz", "Downsampled Signal, Fs = 4410 Hz"})

Now remove the 60 Hz tone using an IIR bandstop filter. Let the stopband have a width of 4 Hz centered at 60 Hz. Choose an IIR filter to achieve a sharp frequency notch, small passband ripple, and a relatively low order.

df4 = designfilt("bandstopiir", ...
    PassbandFrequency1=55, ...
    StopbandFrequency1=58, ...
    StopbandFrequency2=62, ...
    PassbandFrequency2=65, ...
    PassbandRipple1=1, ...
    StopbandAttenuation=60, ...
    PassbandRipple2=1, ...
    DesignMethod="ellip", ...
    SampleRate=Fs);

Analyze the magnitude response.

opts = filterAnalysisOptions;
opts.FrequencyRange = "freqvector";
opts.FrequencyVectorUnits = "Hz";
opts.FrequencyVector = Fds(Fds>F(2));
opts.FrequencyScale = "log";
filterAnalyzer(df4,AnalysisOptions=opts)

Perform zero-phase filtering to avoid phase distortion. Use the filtfilt function to process the data.

ybs = filtfilt(df4,yds);

Finally, upsample the signal to bring it back to the original audio sample rate of 44.1 kHz which is compatible with audio sound cards.

yf = interp(ybs,10);
Fs = Fs*10;

Take a final look at the spectrum of the original and processed signals. The high frequency noise floor and the 60 Hz tone have been attenuated by the filters.

[Pfinal,Ffinal] = pwelch(yf,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal, ...
    {"Original Signal","Final Filtered Signal"})

If your computer has a sound card, you can listen to the signal before and after processing. As mentioned above, the end result is that you have effectively attenuated the 60 Hz hum and the high-frequency noise in the audio file.

% To play the original signal, uncomment the next two lines
% hplayer = audioplayer(y, Fs);
% play(hplayer)

% To play the noise-reduced signal, uncomment the next two lines
% hplayer = audioplayer(yf, Fs);
% play(hplayer);

Differentiating a Signal

The MATLAB diff function differentiates a signal with the drawback that you can potentially increase the noise levels at the output. A better option is to use a differentiator filter that acts as a differentiator in the band of interest, and as an attenuator at all other frequencies, effectively removing high frequency noise.

As an example, analyze the speed of displacement of a building floor during an earthquake. Displacement or drift measurements were recorded on the first floor of a three story test structure under earthquake conditions and saved in the quakedrift.mat file. The length of the data vector is 10e3, the sample rate is 1 kHz, and the units of the measurements are cm.

Differentiate the displacement data to obtain estimates of the speed and acceleration of the building floor during the earthquake. Compare the results using diff and an FIR differentiator filter.

load quakedrift.mat 
Fs  = 1000;                 % Sample rate
dt = 1/Fs;                  % Time differential
t = (0:length(drift)-1)*dt; % Time vector

Design a 50th-order differentiator filter with a passband frequency of 100 Hz, which is the bandwidth over which most of the signal energy is found. Set the stopband frequency of the filter to 120 Hz.

df5 = designfilt("differentiatorfir", ...
    FilterOrder=50, ...
    PassbandFrequency=100, ...
    StopbandFrequency=120, ...
    SampleRate=Fs);

The diff function can be seen as a first order FIR filter with response $$H(Z)=1-Z^{-1}$$. Use Filter Analyzer to compare the magnitude response of the 50th-order differentiator FIR filter and the response of the diff function. Clearly, both responses are equivalent in the passband region (from 0 to 100 Hz). However, in the stopband region, the 50th-order filter attenuates components while the diff response amplifies components. This effectively increases the levels of high frequency noise.

fa = filterAnalyzer(df5, ...
    MagnitudeMode="zerophase", ...
    FilterNames="FIR_Differentiator_50th_Order");
addFilters(fa,[1 -1],1, ...
    SampleRate=Fs, ...
    FilterNames="Response_of_Diff_Function");

Differentiate using the diff function. Add zeros to compensate for the missing samples due to the diff operation.

v1 = diff(drift)/dt;
a1 = diff(v1)/dt;

v1 = [0; v1];
a1 = [0; 0; a1];

Differentiate using the 50th order FIR filter and compensate for delay.

D = mean(grpdelay(df5)); % Filter delay
v2 = filter(df5,[drift; zeros(D,1)]);
v2 = v2(D+1:end);
a2 = filter(df5,[v2; zeros(D,1)]);
a2 = a2(D+1:end);
v2 = v2/dt;
a2 = a2/dt^2;

Plot a few data points of the floor displacement. Plot also a few data points of the speed and acceleration as computed with diff and with the 50th order FIR filter. Notice how the noise has been slightly amplified in the speed estimates and largely amplified in the acceleration estimates obtained with diff.

helperFilterIntroductionPlot2(t,drift,v1,v2,a1,a2)

Integrating a Signal

A leaky integrator filter is an all-pole filter with transfer function $$H(Z)=1/[1-c{Z}^{-1}]$$ where $$c$$ is a constant that must be smaller than 1 to ensure stability of the filter. It is no surprise that as $$c$$ approaches one, the leaky integrator approaches the inverse of the diff transfer function. Apply the leaky integrator to the acceleration and speed estimates obtained in the previous section to get back the speed and the drift respectively. Use the estimates obtained with the diff function since they are noisier.

Use a leaky integrator with $$a=0.999$$. Plot the magnitude response of the leaky integrator filter. The filter acts as a lowpass filter effectively eliminating high-frequency noise.

filterAnalyzer(1,[1 -0.999],SampleRate=Fs)

Filter the velocity and acceleration with the leaky integrator. Multiply by time differential.

v_original = v1;
a_original = a1;

d_leakyint = filter(1,[1 -0.999],v_original);
v_leakyint = filter(1,[1 -0.999],a_original);

d_leakyint = d_leakyint * dt;
v_leakyint = v_leakyint * dt;

Plot the displacement and speed estimates and compare to the original signals.

helperFilterIntroductionPlot3(t,drift,d_leakyint,v_original,v_leakyint)

You can also integrate a signal using the cumsum and cumtrapz functions. Results will be similar to those obtained with the leaky integrator.

Conclusions

In this example you learned about linear and nonlinear phase filters and you learned how to compensate for the phase delay introduced by each filter type. You also learned how to apply filters to remove unwanted frequency components from a signal, and how to downsample a signal after limiting its bandwidth with a lowpass filter. Finally, you learned how to differentiate and integrate a signal using digital filter designs. Throughout the example you also learned how to use analysis tools to look at the response and group delay of your filters.

For more information on filter applications, see the Signal Processing Toolbox™ documentation. For more information on how to design digital filters see the Practical Introduction to Digital Filter Design example.

References

Proakis, J. G., and D. G. Manolakis. Digital Signal Processing: Principles, Algorithms, and Applications. Englewood Cliffs, NJ: Prentice-Hall, 1996.

Orfanidis, S. J. Introduction To Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1996.

Appendix

The following helper functions are used in this example: