This site uses cookies to deliver our services and to show you relevant ads and job listings. By using our site, you acknowledge that you have read and understand our Cookie Policy , Privacy Policy , and our Terms of Service . Your use of the Related Sites, including DSPRelated.com, FPGARelated.com, EmbeddedRelated.com and Electronics-Related.com, is subject to these policies and terms.

I’m writing this article in a room with a bunch of other people talking, and while sometimes I wish they would just SHUT UP, it would be better if I could just filter everything out. Filtering is one of those things that comes up a lot in signal processing. It’s either ridiculously easy, or ridiculously difficult, depending on what it is that you’re trying to filter.

I’m going to show you a one-line algorithm that can help in those cases where it’s ridiculously easy. Interested?

Here’s a practical example of a filtering problem: you have a signal you care about, but it has some noise that snuck in through a little hole in the wall, and darnit, you’re stuck with the noise.

import matplotlib.pyplot as plt import numpy as np np.random.seed(123456789) # repeatable results f0 = 4 t = np.arange(0,1.0,1.0/65536) mysignal = (np.mod(f0*t,1) < 0.5)*2.0-1 mynoise = 1.0*np.random.randn(*mysignal.shape) plt.figure(figsize=(8,6)) plt.plot(t, mysignal+mynoise, 'gray', t, mysignal, 'black');

And the question is, can we get rid of the noise?

Here we’ve got a lot of noise: relative to a square-wave signal with peak-to-peak amplitude of 2.0, the peak-to-peak noise waveform is in the 6-7 range with even larger spikes. Yuck.

The signal processing approach is to look at the frequency spectrum.

def spectrum(x): return np.abs(np.fft.ifft(x)) mysignal_spectrum = spectrum(mysignal) mynoise_spectrum = spectrum(mynoise) N1 = 500 f = np.arange(N1) plt.figure(figsize=(8,4)) plt.plot(f,mysignal_spectrum[:N1]) plt.plot(mynoise_spectrum[:N1]) plt.xlim(0,N1) plt.legend(('signal','noise')) plt.xlabel('frequency') plt.ylabel('amplitude')

Now this is some good news for signal processing. The signal spectrum content really stands out above the noise. Square waves have harmonics that are odd multiples of the fundamental frequency: here 4Hz, 12Hz, 20Hz, 28Hz, etc. The harmonic content decays as 1/f. White noise, on the other hand, has a roughly flat spectrum that is spread out across all harmonics, so it is lower in amplitude for the same amount of energy. Kind of like making a peanut butter sandwich — if you put a glob of peanut butter on a slice of bread, it’s pretty thick, but if you spread it over the bread evenly, it’s much thinner, even though it’s the same amount of peanut butter. Mmmm. I hope you like peanut butter; this is making me hungry.

This effect is easier to see on a log scale:

z0 = 1e-10 # darnit, we have to add a small number so log graphs don't barf on zero for N1 in [500, 8000]: f = np.arange(N1) plt.figure(figsize=(8,4)) plt.semilogy(f,z0+mysignal_spectrum[:N1],f,z0+mynoise_spectrum[:N1]) plt.xlim(0,N1) plt.ylim(1e-4,1) plt.legend(('signal','noise')) plt.xlabel('frequency') plt.ylabel('amplitude')

What we’d really like is to keep the low frequencies only; most of the signal energy is only in the lowest harmonics, and if we want to get rid of the noise, there’s no point in keeping the higher frequencies where the noise energy is much higher. Looks like the crossover here is roughly around 500Hz; below 500Hz the signal dominates, and above that point the noise dominates.

So let’s do it! We’ll use a low-pass filter to let the low frequencies pass through and block the high frequencies out.

In an ideal world, we’d use a low-pass filter with a very sharp cutoff, in other words one that lets everything through below 500Hz and nothing through above 500Hz. But in practice, sharp-cutoff filters are challenging to implement. It’s much easier to create a gradual-cutoff filter, and the simplest is a single-pole infinite impulse response (IIR) low-pass filter, sometimes called a exponential moving average filter.

We’re going to use a filter which has a transfer function of \( H(s) = \frac{1}{\tau s + 1} \). This is a first-order system, which I’ve talked about in a previous article, in case you’re interested in more information. But here we’re just going to use it.

We’ll start with the differential equation of this filter, which is \( \tau \frac{dy}{dt} + y = x \).

We can take this continuous-time equation and convert to a discrete-time system with timestep \( \Delta t \), namely \( \tau \frac{\Delta y}{\Delta t} + y = x \), and if we solve for \( \Delta y \) we get \( \Delta y = (x-y)\frac{\Delta t}{\tau} \). Blah blah blah, math math math, yippee skip. Here’s the upshot, if all you care about is some computer code that does something useful:

// // SHORTEST USEFUL ALGORITHM ****EVER**** // y += alpha * (x-y); // // # # # ## ## ##### # // # # # # # # # # # // # # # # # # # # // # # ## ## # # //

That’s it! One line of C code. The input is x . The output is y . The filter constant is alpha . Every timestep we add alpha times the difference between input and output. So simple.

Now we just have to calculate what alpha is. The approximate answer is \( \alpha = \frac{\Delta t}{\tau} \): we just divide \( \Delta t \) by the time constant \( \tau \). The more accurate answer is \( \alpha = 1 - e^{-\Delta t/\tau} \approx \frac{\Delta t}{\tau} \), with the approximation being more accurate for time constants much longer than the timestep, namely \( \tau \gg \Delta t \).

Let’s see it in action. If we want to have a filter with a cutoff of 500Hz (in other words, pass through frequencies much lower than 500Hz, block frequencies much higher than 500Hz), we use a value of \( \tau \) = 1/500 = 2msec.

def lpf(x, alpha, x0 = None): y = np.zeros_like(x) yk = x[0] if x0 is None else x0 for k in xrange(len(y)): yk += alpha * (x[k]-yk) y[k] = yk return y delta_t = t[1]-t[0] tau = 0.002 alpha = delta_t / tau plt.figure(figsize=(8,6)) y1 = lpf(mysignal+mynoise, alpha, x0=0) plt.plot(t,y1)

Hey, that’s not bad, especially considering that the original signal including noise had a huge amount of noise on it. There is always a tradeoff in filtering: if you want less noise, you can use a lower cutoff frequency, at the cost of more phase lag. Phase lag is a frequency-dependent delay that rounds off sharp edges and plays havoc with feedback control systems, if they aren’t designed to take that phase lag into account. We could use a lower cutoff frequency here, say 100Hz → \( \tau \) = 10 msec:

tau = 0.010 alpha = delta_t / tau plt.figure(figsize=(8,6)) y2 = lpf(mysignal+mynoise, alpha, x0=0) plt.plot(t,y2)

There you go: less noise, more rounding / distortion / phase lag / whatever you want to call it.

Here are the signal spectra for signals shown in the two previous graphs:

N1 = 500 f = np.arange(N1) plt.figure(figsize=(8,4)) plt.semilogy(f,spectrum(y1)[:N1],'.',f,spectrum(y2)[:N1],'+') plt.xlim(0,N1) plt.ylim(1e-4,1) plt.legend(('tau = 2msec','tau = 10msec')) plt.xlabel('frequency') plt.ylabel('amplitude')

The higher time constant filter has a lower cutoff frequency, so the remaining amplitude of the output signal is lower.

Now in practice, this example we’ve used contains a huge amount of noise on the input signal. We always hope for signals that are less noisy than that. Here’s a more reasonable example, with only 10% of the noise we showed earlier:

mynoise2 = 0.1*np.random.randn(*mysignal.shape) alpha=delta_t/0.0005 lpf_out2 = lpf(mysignal+mynoise2, alpha=alpha, x0=0) fig=plt.figure(figsize=(8,10)) ax=fig.add_subplot(3,1,1) ax.plot(t,mysignal+mynoise2,t,mysignal) ax.set_ylabel('original signal') ax2=fig.add_subplot(3,1,2) ax2.plot(t,lpf_out2) ax2.set_ylabel('lpf output') ax3=fig.add_subplot(3,1,3) ax3.plot(t,lpf_out2) ax3.set_ylim(-1.05,-0.95) ax3.set_ylabel('lpf output, zoom')

Here we can go all the way out to a cutoff frequency of 2000 Hz (time constant of 0.5 msec), get fairly nice edges in our signal, but still remove most of the noise. In general the noise content gets reduced in amplitude by a factor of about \( \sqrt{\alpha} \): here we have \( \alpha = \Delta t/\tau = 1.0/f\tau \) where f is the sampling frequency, and for f = 65536Hz we have \( \alpha = 1/32.768 \), so should expect about a reduction to about 1/5.7 of the original noise level. We might have had a peak-to-peak noise of about 0.5 in the original signal, and we’re down to maybe 0.08 now, which is around 1/6 of the original noise level. Measuring peak noise is kind of subjective, but we’re in the ballpark of what theory predicts.

The interesting thing here is the value \( \alpha \) drops if the timestep \( \Delta t \) gets smaller: we can improve the noise filtering if we sample more often. This is called oversampling and it can be used, up to a point, to reduce the impact of noise on signals of interest, by sampling them more than normal and then using low-pass filters to get rid of the noise. Oversampling has diminishing returns, because if we want half the noise amplitude, we have to sample four times as often, and compute the filter four times as often. It also only works on noise. If our unwanted signal is, say, a Rick Astley song, with the same general type of frequency spectrum content of our signal of interest, then we’re totally out of luck if we try to filter it out, and we just have to hope the person playing it has some sense of decency and will turn it off.

Fixed-point implementation of a low-pass filter

The one line of code posted above:

// curtains part, angels sing.... y += alpha * (x-y); // now back to normal

is applicable for floating-point math. Most of the lower-end processors don’t have floating-point math, and have to emulate it, which takes time and CPU power. If you want to work in fixed-point instead, the code to implement a low-pass filter looks like this:

y_state += (alpha_scaled * (x-y)) >> (Q-N); y = y_state >> N;

where alpha_scaled = \( 2^Q \alpha \) for some fixed-point scaling Q, and you have to keep N additional state bits around to accurately represent the behavior of the low-pass filter. The value of N is more or less logarithmic in alpha, e.g. \( N \approx - \log_2 \alpha \), so if you want to use very long time constants, you have to keep around more bits in your filter state. I worked with a software contractor around 2002 or 2003, who was convinced I was wrong about this; if I wanted to filter a 16-bit ADC reading, what was wrong about storing everything in 16 bits? Why would I have to use 32-bit numbers? But you do; try it without adding extra bits: you’ll find that there is a minimum difference in input that is needed in order to get the output to change, and small changes in the input will make the output look stair-steppy. This is an unwanted quantization effect.

Don’t believe me? Let’s try it, first with floating-point, then with fixed-point, on a pair of signals. Signal x1 will be a square wave plus some unwanted high-frequency stuff; signal x2 will be the plain square wave. Here we go:

# Floating-point low-pass filter x1 = mysignal + 0.1*np.sin(800*t) x2 = mysignal*1.0 alpha = 1.0/1024 y1 = lpf(x1,alpha,x0=0) y2 = lpf(x2,alpha,x0=0) fig=plt.figure(figsize=(8,6)) ax1 = fig.add_subplot(2,1,1) ax1.plot(t,x1,'gray',t,x2,'black') ax1.set_ylabel('LPF input') ax2 = fig.add_subplot(2,1,2) ax2.plot(t,y1,'gray',t,y2,'black') ax2.set_ylabel('LPF output');

Now we’ll do a low-pass filter in fixed-point. We’ll plot the filter output (scaled back from integer counts to a floating-point equivalent), along with the error between the fixed-point and floating-point implementations.

# fixed-point low-pass filter def lpf_fixp(x, alpha, Q, N, x0 = None): y = np.zeros_like(x, dtype=int) yk = x[0] if x0 is None else x0 y_state = yk << Q alpha_scaled = np.round(alpha*(1 << Q)).astype(int) for k in xrange(len(y)): y_state += (alpha_scaled * (x[k]-yk)) >> (Q-N) yk = y_state >> N y[k] = yk return y from matplotlib import gridspec Q=12 x1q = np.round(x1*(1 << Q)).astype(int) x2q = np.round(x2*(1 << Q)).astype(int) Nlist = [6,4,2,0] n_axes = len(Nlist) for k,N in enumerate(Nlist): fig = plt.figure(figsize=(8,6)) gs = gridspec.GridSpec(7, 1) ax = fig.add_subplot(gs[:6,0]) y1q = lpf_fixp(x1q, alpha, Q=Q, N=N, x0=0) y2q = lpf_fixp(x2q, alpha, Q=Q, N=N, x0=0) y1qf = y1q*1.0/(1 << Q) # convert back to floating point y2qf = y2q*1.0/(1 << Q) # convert back to floating point ax.plot(t,x2,'#f0f0f0') # show square wave for reference ax.plot(t,y1qf,'red',t,y2qf) ax.set_ylim([-1.1,1.1]) ax.set_xticklabels([]) ax.legend(('sq wave input','y1 output','y2 output'),loc='best') ax2 = fig.add_subplot(gs[6,0]) ax2.plot(t,y1qf-y1,'red',t,y2qf-y2) ax2.set_ylabel('error') fig.suptitle('fixed-point lpf, alpha=%.6f, Q=%d, N=%d' % (alpha,Q,N))

The proper number of extra state bits is \( N = -\log_2 \frac{1}{1024} = 10 \), but we used 6, 4, 2, and 0 in these examples. You’ll notice a couple of things:

the lower N is, the worse the error gets

for N = 6 the error is small but still present

for N = 4, 2 and 0 the error is visually noticeable

there is a DC bias: the output is always lower than the input. (This is true because a shift right represents a rounding down in the integer math.)

the error for the y2 signal is always much higher; this is the output for the plain square wave with no harmonic content.

This last point is kind of interesting. The error caused by not having enough state bits is quantization error, and when high-frequency signals are present, it helps alleviate quantization error, in part because these signals are helping the filter state “bounce” around. This effect is called dithering, and it’s a trick that can be used to reduce the number of state bits somewhat. But to use it, you really need to know what you’re doing, and it also takes extra work (both up-front design time, and also CPU time) on an embedded system to generate the right dithering signal.

Aside from the dithering trick, you’re either stuck with using more bits in your low-pass filters, or you can use a larger timestep \( \Delta t \) so the filter constant \( \alpha \) doesn’t get so small.

Wrapup

Well, aside from the fixed-point quirks, the basic one-pole low-pass filter algorithm is pretty simple.

If you need more filtering than a one-pole low-pass filter can provide, for example you have lots of 1kHz noise on a 3Hz signal, another thing you can do is to cascade two of these one-pole low-pass filters (in other words, filter twice). It’s not as good as an optimally-designed two-pole low-pass filter, but it’s stable and easy to implement.

That’s it for today – happy filtering!

© 2015 Jason M. Sachs, all rights reserved.