4
$\begingroup$

The question is : What is wrong with averaging as low pass filter ?

The details : I want to lowpass filter a signal to downsample it. The constraints are : I have no RAM available and I work in streaming, therefore I cannot use a weighting window to filter. Consequently, I have to do either an average on the last X samples, or a classic filter like :

signal_t+1 = signal_t * (1-factor) + new_measurement * factor

with factor a value adapted to the sampling. What is more, I know there are signals around a frequency just above the downsampling frequency.

Is it better to use the average or to use the 1-order filter described above ?

More details and tests : I work with real samples, not complex. Sampling frequency = 200 Hz, frequency after downsampling = 4Hz.

Let's compare a 1-order filter, with a cutoff frequency of 1Hz to limitate aliasing, with two averaging filters. The first averaging filter is the average of 50 samples to downsample from 200Hz to 4Hz. The second averaging filter is an average on 66 samples to get as much rejection as with the 1-order filter.

comparison of the 4 filters

I think that the filter "average on 66 samples" is the best : The 50-samples average filter has a rejection higher than the 1-order filter, and I really need to limitate aliasing. But the 66-samples filter has a cutoff frequency larger than the 1-order filter and a rejection at least equal after 2Hz. What is more, the 66-samples filter cuts the frequencies around 3Hz and I know there will be signals on 3Hz that will make my signal noisy after the 4Hz downsampling - I work with real samples, not complex.

Then, I would chose the 66-samples averaging filter.

Is it something wrong in the logic above ? Would you have suggestions ?

Here is the code :

import numpy as np
import matplotlib.pyplot as plt

nfft = 8192
fsampl = 200
dtsampl = 1 / 200
fcutoff = 1
dtcutoff = 1 / (2 * np.pi * fcutoff)
faverage1 = 4
faverage2 = 3

print("generate signal")
raw_signal = np.zeros(nfft)
raw_signal[int(nfft/2)]=1

print("filter signal")
filtered_signal = 0.*np.copy(raw_signal)
for i in range(1, np.size(raw_signal)):
    filtered_signal[i] = filtered_signal[i-1]+ (raw_signal[i-1]-filtered_signal[i-1]) * dtsampl / dtcutoff

averaged_signal_1 = np.copy( raw_signal )
num_average = fsampl / faverage1
for i in range(1, np.size(raw_signal)):
    averaged_signal_1[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

averaged_signal_2 = np.copy( raw_signal )
num_average = fsampl / faverage2
for i in range(1, np.size(raw_signal)):
    averaged_signal_2[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

print("get FFT")
raw_signal_fft = np.fft.rfft(raw_signal)
filtered_signal_fft = np.fft.rfft(filtered_signal)
averaged_signal1_fft = np.fft.rfft( averaged_signal_1)
averaged_signal2_fft = np.fft.rfft( averaged_signal_2)

raw_signal_fftlog = 20*np.log10(np.abs(raw_signal_fft))
filtered_signal_fftlog = 20*np.log10(np.abs(filtered_signal_fft))
averaged_signal1_fftlog = 20*np.log10(np.abs(averaged_signal1_fft))
averaged_signal2_fftlog = 20*np.log10(np.abs(averaged_signal2_fft))

print("plot")
freq_axe= np.linspace( 0, fsampl / 2, int( nfft / 2 ) + 1 )
plt.plot(freq_axe, raw_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal_fftlog, 'r-')
plt.plot(freq_axe, averaged_signal1_fftlog, '-')
plt.plot(freq_axe, averaged_signal2_fftlog, '-')
plt.plot([0.02, 100], [-3, -3], 'k:')
plt.legend(['raw signal (dirac)', '1-order 1Hz filter', 'average on 50', 'average on 66'],  loc="upper right")
plt.ylabel('amplitude (dB)')
plt.xlabel('frequencies (Hz)')
plt.xlim([0,8])
plt.ylim([-30,2])

plt.show()
$\endgroup$
  • 2
    $\begingroup$ You had to use 66 samples of average for your moving average filter to match the performances of the 1st order LP. That's a 66 times longer delay. To quote The Scientist & Engineer's Guide to Digital Signal Processing, "the moving average is an exceptionally good smoothing filter (the action in the time domain), but an exceptionally bad low-pass filter (the action in the frequency domain)" $\endgroup$ – Florent Aug 24 '17 at 23:29
  • $\begingroup$ I agree with @Florent Ecochard about the delay versus performance of the moving average filter and would not recommend using it as an anti-aliasing filter (nor would I recommend using a first order anti-aliasing filter). You could potentially use the "zeros" in the magnitude response of the filter to suppress specific frequencies very effectively. $\endgroup$ – user883521 Aug 25 '17 at 10:31
  • $\begingroup$ Thank you for your quick answers ! @user883521 what would you recommend to handle both the delay and the need of anti-aliasing at 3Hz ? I have no need to work in the frequency domain, but I would like to minimize delay if possible. $\endgroup$ – Mac Aug 25 '17 at 11:32
  • $\begingroup$ The issue with the 1st order filter is the attenuation it can achieve at 3 Hz. In your case with a cutoff at 1 Hz, the content at 3 Hz is roughly attenuated by a factor 1/3. If this performance is acceptable then a first order filter is an easy, acceptable solution. $\endgroup$ – user883521 Aug 25 '17 at 12:27
  • $\begingroup$ How much of a problem is the delay? What are the exact requirements with respect to the resampling and attenuation? Is real-time (streaming?) processing needed? You could select a higher order filter, or possibly a higher order linear phase filter and correct for the delay or implement a resampling filter, etc. I do not directly have experience with doing this in realtime so there are probably better options that I've missed. $\endgroup$ – user883521 Aug 25 '17 at 12:40
1
$\begingroup$

The question was : What is best between mooving average and 1st order filter when no RAM is available ?

The solution would be : None of them, except for a specific need : rejection of a particular frequency. Otherwise, it is better to use filter of higher orders than 1st order.

  • Mooving average :
    • Drawbacks : delay and poor rejection.
    • Assets : huge rejection at several particular frequencies
  • 1-order filter
    • Drawbacks : rejection which is poor whatever frequency
    • Assets : no delay
  • Higher order filters. Higher order filters include : cascaded 1st order filters (here under an exemple with 5 cascaded 1st order filters), butterworth filters, etc.
    • Assets : better rejection than 1st order filter and no need of much RAM. Small delay.

View of the result : enter image description here

Thank you @user883521 and @Florent Ecochard for your help !

Here is the Pyhton code for a 5th order filter (for other orders, cf links here under. for other filters, see signal.lfilter)

import numpy as np
import matplotlib.pyplot as plt

nfft = 8192
fsampl = 200
dtsampl = 1 / 200

fcutoff_1st_order = 1
dtcutoff_1storder = 1 / (2 * np.pi * fcutoff_1st_order)
fcutoff_5thorder =2.635
faverage2 = 3

print("generate signal")
raw_signal = np.zeros(nfft)
raw_signal[int(nfft/2)]=1

print("filter signal")
filtered_signal = 0.*np.copy(raw_signal)
for i in range(1, np.size(raw_signal)):
    filtered_signal[i] = filtered_signal[i-1]+ (raw_signal[i-1]-filtered_signal[i-1]) * dtsampl / dtcutoff_1storder
filtered_signal5= np.copy(raw_signal)
x = np.exp( -2 * np.pi * fcutoff_5thorder / fsampl )
for i in range(5, np.size(raw_signal)):
    filtered_signal5[i] = np.power( 1 - x, 5 ) * raw_signal[i - 1] \
                          + 5 * np.power( x, 1 ) * filtered_signal5[i - 1] + \
                          -10 * np.power( x, 2 ) * filtered_signal5[i - 2] + \
                          +10 * np.power( x, 3 ) * filtered_signal5[i - 3] + \
                          -5 * np.power( x, 4 ) * filtered_signal5[i - 4] + \
                          +np.power( x, 5 ) * filtered_signal5[i-5]
averaged_signal_66 = np.copy( raw_signal )
num_average = fsampl / faverage2
for i in range(1, np.size(raw_signal)):
    averaged_signal_66[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

print("get FFT")
raw_signal_fft = np.fft.rfft(raw_signal)
filtered_signal_fft = np.fft.rfft(filtered_signal)
filtered_signal5_fft = np.fft.rfft(filtered_signal5)
averaged_signal66_fft = np.fft.rfft( averaged_signal_66 )
raw_signal_fftlog = 20*np.log10(np.abs(raw_signal_fft))
filtered_signal_fftlog = 20*np.log10(np.abs(filtered_signal_fft))
filtered_signal5_fftlog = 20*np.log10(np.abs(filtered_signal5_fft))
averaged_signal66_fftlog = 20 * np.log10( np.abs( averaged_signal66_fft ) )

print("plot")
freq_axe= np.linspace( 0, fsampl / 2, int( nfft / 2 ) + 1 )
plt.plot(freq_axe, raw_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal5_fftlog, '-')
plt.plot( freq_axe, averaged_signal66_fftlog, '-' )
plt.plot([0.02, 100], [-3, -3], 'b--')
plt.plot([2, 2],[-50, 0],  'k:')
plt.plot([3, 3],[-50, 0],  'k--')
plt.legend(['raw signal (dirac)', '1-order 1Hz filter', '5-order 1Hz filter', 'mooving average 66 samples',
            '-3dB', 'aliasing freq', 'unwanted signals frequency'],  loc="upper right")
plt.ylabel('amplitude (dB)')
plt.xlabel('frequency (Hz)')
plt.xlim([0,8])
plt.ylim([-30,2])

plt.show()

useful links : http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch14.pdf http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch15.pdf http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch19.pdf

$\endgroup$
  • $\begingroup$ In addition to the fifth order filter you use above you could also use SciPy's signal.iirfilter or signal.iirdesign to design a digital filter with some more direct control over the properties. See the help of signal.lfilter on how to implement it recursively. $\endgroup$ – user883521 Aug 25 '17 at 17:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.