Sometimes, when working with scientific data, you have noisy data that you need to extract low-frequency components from.
Imagine, for example, that for a project you have recorded some audio clips that have a high-frequency "hiss" artifact from your recording equipment. Using a convolution, you can eliminate the hiss while not destroying the underlying signal that you're interested in.
For my current research project on an adaptive optics instrument, we needed to smooth a signal as part of our troubleshooting process to ensure we had the pattern we expected at low frequencies.
For this, we used IPython (with NumPy, SciPy, Matplotlib and friends), and AstroPy (an up-and-coming library providing implementations of common functionality for astronomers).
Open up IPython
The IPython notebook makes a lot of things easier, from keeping track of what you've tried to writing blog posts like this one. Open it up, and let's get started.
Use the magic command
%pylab inline to bring the usual suspects (
plt, etc.) into the namespace and tell IPython you want inline plots.
Import the convolution functionality from AstroPy.
from astropy.convolution import convolve, Box1DKernel # n.b. this overrides pylab's convolve()
Create a noisy signal to smooth
N = 1000 # number of samples we're dealing with dt = 1.0 / 100.0 # 100 samples / sec
Create some nice random noise. By default, random noise is in the range [0.0, 1.0], so shift it down by 0.5 such that it's equal parts positive and negative.
noise_ts = 3 * (np.random.rand(N) - 0.5) # center at 0 plot(timesteps, noise_ts, 'b.') ylim(-8, 8)
Now, we create our sine wave that we're going to mess up with noise.
A = 5.0 frequency = 0.5 # Hz omega = 2 * np.pi * frequency timesteps = np.linspace(0.0, N*dt, N) signal = A * np.sin(omega * timesteps)
plot(timesteps, signal) ylim(-8, 8)
Add the noise to the signal and get a much noisier wave.
noisy_signal = noise_ts + signal
plot(timesteps, noisy_signal, 'b.') ylim(-8, 8)
Smooth the noisy signal with
Boxcar smoothing is equivalent to taking your signal %%x[t]%% and using it to make a new signal %%x'[t]%% where each element is the average of w adjacent elements. In other words, for w = 5, element %%x'%% will be given by
It's not perfect, but it's frequently good enough. Here we're using AstroPy's
convolve function with a "boxcar" kernel of width w = 11 to eliminate the high frequency noise.
smoothed_signal = convolve(noisy_signal, Box1DKernel(11))
plot(timesteps, noisy_signal, 'b.', alpha=0.5, label="Noisy") plot(timesteps, smoothed_signal, 'r', label="Smoothed") ylim(-8, 8)
Let's make the figure we had at the beginning of this post, too.
figure(figsize=(12,4)) ax = subplot(121) title("Noisy Signal") plot(timesteps, noisy_signal) xlim(0, 4) subplot(122, sharey=ax) title("Nice, Smooth Signal") plot(timesteps, smoothed_signal, 'r', label="Smoothed") xlim(0, 4)
Different shaped kernels can provide useful behavior. Convolution can also be performed in two dimensions. For example, if you want to smooth an image, you can use the
Box2DKernel or any of the other kernels available in AstroPy. (If you are familiar with Photoshop, the
Gaussian2DKernel is analogous to the useful "Gaussian Blur" filter.)
The documentation for this part of AstroPy is particularly good, and I highly recommend reading it!