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.
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()