This chapter provides a quick overview of digital signal processing DSP . DSP is inherent part of many measurement system and at least post processing of data before actual analysis and is therefore included in this course.
The purpose of signal processing is to improve the quality of a signal by suppressing noise. It can also be used to get rid off redundant information to store and transmit data more efficiently.
SciPy has extensive DSP functionality in
scipy.signal
namespace.
Filtering aims to remove unwanted noise from a signal or separate several sources from one measured signal.
Filters are not perfect so in practice filtered data also contains filtered frequencies, but they are greatly attenuated.
Filters can be optimized for either time domain or frequency domain performance, but it is not possible to optimize a filter for both applications.
The frequency response of a filter can be obtained from its impulse response by Fourier transforming it and the step response is the cumulative sum of the impulse response. Chapter 14 of (Smith 1997) provides a good introduction to filters.
Every filter causes some amount of distortion to filtered signal. It is important to be aware of the properties of the filter you use.
The properties of a filter are characterized by the following concepts and Figure 8.1
SciPy filter design and analysis code is in
scipy.signal
namespace.
There are several function to design and analyze filters.
However SciPy doesn’t have good functions for plotting amplitude
response or phase response.
So here are two useful functions that can be
used to plot filter properties.
You’ll need to import
scipy.signal
in order to use them.
from scipy import signal
def plot_freqz(b,a=1):
w,h = signal.freqz(b,a)
h_dB = 20 * log10 (abs(h))
plot(w/pi, h_dB)
ylim([max(min(h_dB), -100) , 5])
ylabel('Magnitude (db)')
xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
title(r'Amplitude response')
def plot_phasez(b, a=1):
w,h = signal.freqz(b,a)
h_Phase = unwrap(arctan2(imag(h),real(h)))
plot(w/pi, h_Phase)
ylabel('Phase (radians)')
xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
title(r'Phase response')
def plot_impz(b, a = 1):
if type(a)== int: #FIR
l = len(b)
else: # IIR
l = 100
impulse = repeat(0.,l); impulse[0] =1.
x = arange(0,l)
response = signal.lfilter(b, a, impulse)
stem(x, response, linefmt='b-', basefmt='b-', markerfmt='bo')
ylabel('Amplitude')
xlabel(r'n (samples)')
title(r'Impulse response')
def plot_stepz(b, a = 1):
if type(a)== int: #FIR
l = len(b)
else: # IIR
l = 100
impulse = repeat(0.,l); impulse[0] =1.
x = arange(0,l)
response = signal.lfilter(b,a,impulse)
step = cumsum(response)
plot(x, step)
ylabel('Amplitude')
xlabel(r'n (samples)')
title(r'Step response')
def plot_filterz(b, a=1):
subplot(221)
plot_freqz(b, a)
subplot(222)
plot_phasez(b, a)
subplot(223)
plot_impz(b, a)
subplot(224)
plot_stepz(b, a)
subplots_adjust(hspace=0.5, wspace = 0.3)
A decibel is dimensionless unit that compares the ratio of numbers in logarithmic scale. It is often used to express the power ratio of two signals.
In signal processing we are usually interested in the amplitude ratio of two signals. Usually before and after filtering.
dB scale | Linear scale |
60 dB | 1000 |
40 dB | 100 |
20 dB | 10 |
0 dB | 1 |
-20 dB | 0.2 |
-40 dB | 0.01 |
-60 dB | 0.001 |
Moving average is a simple lowpass filter. It replaces samples with the average of n samples.
The commonly used variant is the central moving average, here is how you calculate it for 5 points:
where \( y_t \) is \( t \) th of the filtered signal and \( x \) is the original signal.
Moving average is an optimal filter for time domain signals that contain white noise. For noise with specific frequency you can get better performance with FIR or IIR filters.
Figure 8.3 shows the effect of number of points used in filtering a signal. Greater number of points provides smoother result, but it also “cuts the corners” of the signal. Finding the correct number of points for each case is done using experimentation and plotting the original signal and filtered signal to the same plot and zooming close enough to rapid changes helps in the process.
SciPy doesn’t have a builtin implementation of a moving average filter, but it is easy to implement it. A moving average of order \( n \) has an impulse response with \( n \) elements that all have the value of \( 1/n \) .
from scipy import signal
def ma(x, n = 5):
b = repeat(1.0/n, n) #Create impulse response
xf = signal.lfilter(b, 1, x) #Filter the signal
return(xf)
Moving average filters have good performance for time domain signals, but their amplitude response is pretty terrible. You can see from Figure 8.2 that the response is not very sharp and there is a lot of ripple in the passband.
Figure 8.2 was created using following python code using the functions we defined earlier.
b = repeat(1.0/11, 11)
plot_filterz(b)
Another advantage of moving average filter is that it can be calculated using an efficient recursive algorithm. See (Smith 1997) Ch 15 p. 282.
The median filter is similar to moving average, but is uses the moving median instead of average. It is nonlinear filter so you can’t calculate its step or amplitude response. Median filter can be used to suppress heavy non-Gaussian noise in time domain signals e.g. randomly occurring high peaks.
scipy.signal.medfilt
in Python.
A comparison of median filter and moving average filter is shown in Figure 8.3
Let’s see how moving average filters with different order and median filter can handle a noisy ramp signal. We’ll first create a ramp signal with added gaussian noise:
ym = hstack([repeat(0, 100), repeat(10, 100), repeat(0, 100)])
noisy = ym+randn(300)
Then filter it with three diffent MA filters and a median filter and plot the results. The example also demonstrates matplotlib’s object oriented API and axis-sharing between subplots.
#Moving average filters
ym11 = ma(noisy)
ym31 = ma(noisy, 31)
#Median filter from scipy.signal
ymed11 = signal.medfilt(noisy, 11)
#Plots with shared x and y-axis
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2,
sharex=True, sharey=True)
ax1.plot(noisy, 'k')
ax1.set_title("a) Noisy signal")
ax1.set_ylim([-5, 15])
ax2.plot(ym11, 'k')
ax2.set_title("b) Moving average filter n = 11")
ax3.plot(ym31, 'k')
ax3.set_title("b) Moving average filter n = 31")
ax4.plot(ymed11, 'k')
ax4.set_title("b) Median filter n = 31")
#subplots_adjust(hspace=0.4)
Figure 8.4 shows the time and frequency domain properties of 100th order Hamming lowpass filter. The amplitude response has sharp roll off and the phase is linear in the passband. There is mild overshoot in the step response. The roll off can be made almost arbitrarily sharp using even higher orders. If you need to filter frequency specific noise from a time domain signal a linear phase FIR filter is a good choice.
Design 100th order lowpass filter using hamming window, b is the impulse response filter stopband is set to 0.2 which equals \( 0.2 \cdot 0.5 \cdot \) samplerate.
n = 100
#Design filter
b = signal.firwin(n, cutoff = 0.2, window = "hamming")
#Plot properties
plot_filterz(b)
Use the filter and plot results.
x = sin(linspace(0, 500, 1024)) + 0.5*cos(linspace(0, 750, 1024))
+ randn(1,1024)*0.2 + 0.2*cos(linspace(0, 10000, 1024));
xfiltered = signal.lfilter(b, 1, x)
plot(x[:200], label = "Original signal")
plot(xfiltered[:200], "--k", label = "Filtered signal")
legend()
Figure 8.5 shows the properties of 12th order Chebyshev II lowpass filter. The roll of is not as sharp as for the FIR filter in Figure 8.4 and there is more overshoot in the step response. The filter however has only 12 coefficients instead of 100 of the FIR filter making it considerably more efficient computationally.
The properties of 12th order Butterworth filter (which less aggressive IIR filter) are shown in Figure 8.6 . It has smooth roll of and doesn’t overshoot as much making it a better candidate for time domain applications.
The list IIR filters that can be designed with SciPy can be found from
scipy.signal
documentation.
Here are a few examples:
Design 12th order chebyshev lowpass filter with stopband attenuation of 80dB filter stopband is set to 0.2 which means \( 0.2 \cdot 0.5 \cdot \) samplerate ( Figure 8.5 ).
from scipy import signal
[b, a] = signal.cheby2(12, 80, 0.2)
plot_filterz(b,a)
Design 12th order lowpass butterworth filter with stopband at 0.2 in Figure 8.6 .
[b, a] = signal.butter(12, 0.2);
plot_filterz(b, a)
The above commands are similar to the ones found from MATLAB, there is
also another option:
iirdesign
command.
Some examples.
Let’s define 2 filters and plot their amplitude response: Elliptic
bandpass filter and chebyshev I highpass filter.
Parameter
wp
defines the pass band and
ws
defines the stop band.
be, ae = signal.iirdesign(wp = [0.05, 0.3], ws= [0.02, 0.35],
gstop= 60, gpass=1, ftype='ellip')
plot_freqz(be,ae)
bb, ab = signal.iirdesign(wp = 0.3, ws= 0.2,
gstop= 60, gpass=1, ftype='cheby1')
plot_freqz(bb,ab)
Here is an example of using different filters on signal Figure 8.9 a. We can see that the signal is contaminated with noise, but we don’t know if it is white noise or if it has a specific frequency.
We start the analysis by plotting the periodogram of the signal ( Figure 8.10 a). It reveals that there is high frequency noise at around 0.9. (The unit is relative to 0.5 \( \cdot \) sampling rate, 0.9 becomes \( 0.9 \cdot 0.5 \cdot \) sample rate in actual units) and the interesting frequencies are clearly below 0.4. This means we can use a lowpass filter with stopband at 0.4 to remove the noise.
Figures 8.9 and 8.10 show the effect of using 100th order FIR hamming filter to the time domain signal and the PSD. Also the effect of 5 point moving average filter is shown for reference. The FIR filter is clearly the right choice for this signal, it preserves the time domain features nicely ( Beware of small phase shift for time critical applications ) and completely eliminates the high frequency noise. The moving average fails to remove the high frequency noise, but “cuts” the peaks from decreasing the amplitude.
The sampling rate of a digital signal can be changed computationally. There are two basic operations:
Decimation
is used to lower to sample rate.
The operation uses an antialising
filter to make sure there aren’t too high frequency components in
the downsampled signal.
signal.decimate
in SciPy.
Interpolation Increases the sampling rate of the signal It is mainly used for digital to analog conversion to get more continuous results for e.g. audio.
NOTE: Interpolation doesn’t fix too low sample rate.
Detrending means removing the slowly changing trends from a higher frequency signal. It is required for spectral analysis and can also be used to remove the effect of sensor drift from time domain measurements.
Some methods for detrending are: