Low-pass filter

Hello! Tell me please if it is possible to apply a low-pass filter to the signal using Sage? And if it is possible, what is a way to do it?
I have read that FFT could be used for this, but haven't got the whole idea.

Thanks.

edit retag close merge delete

Sort by » oldest newest most voted

Yes, using the functions in the scipy.signal module. You're probably better off searching the scipy docs and mailing lists for help if you're going to be doing signal processing.


# modified for Sage from an example by Ivo Maljevic
# http://mail.scipy.org/pipermail/scipy-user/attachments/20090917/d2b7a007/attachment-0001.py

import numpy as np
import scipy.signal

t = np.linspace(0,1,1000) # f_s = 1000 Hz

f_1 = 100 # 100 Hz component
f_2 = 300 # 300 Hz component

f = np.linspace(0, 1000,1000)
v = np.sin(2*np.pi*f_1*t) + np.cos(2*np.pi*f_2*t) # f_1 and f_2 freq components added together
V = np.abs(np.fft.fft(v))**2
p = list_plot(zip(f, V),color="blue",plotjoined=True)
p.show()

(b,a) = scipy.signal.butter(12, 200.0/500.0, btype='low') # filter out anything above 200 Hz
v = scipy.signal.lfilter(b,a,v)
V = np.abs(np.fft.fft(v))**2
p2 = list_plot(zip(f,V),color='red',plotjoined=True)
p2.show()


I'd definitely recommend reading up on the different filters before applying them to any important data. N.B. I didn't verify that there were no minor errors in the above (I always have off-by-one errors in this stuff anyway, so my confirmation wouldn't be worth much..)

more

Finally I decided to create a simple as a log hand-made FFT filter. The concept is as follows:

• Perform the direct Fourier transform.
• Zero out the coefficients corresponding to the target frequencies
• Perform the inverse Fourier transform.

I borrowed this idea described here and wrote the following function:

def fft_lowpass_filter(signal, min_period=2):
"""Applies a simple as a log FFT lowpass filter to a signal.

INPUT:

- signal -- a list of data representing the sampled signal;

- min_period -- a threshold for the lowpass filter (the minimal period
of the oscillations which should be left intact) expressed in a number of
samples per one full oscillation.

EXAMPLES:

If you have a sampling frequency equal to 1 Hz (one sample per second)
and want to filter out all the frequencies higher than 0.5 Hz (two samples
per one oscillation period) you should call this function like this::

sage: fft_lowpass_filter(signal, min_period=2)

If you, for example, have a sampling frequency of 1 kHz (1000 samples per
second) and wish to filter out all frequencies hihger than 50 Hz, you should
use the value of min_period = <sampling frequency> / <cutoff frequency> = 1000 Hz / 50 Hz = 20::

sage: fft_lowpass_filter(signal, min_period=20)

"""
import scipy

signal_fft = scipy.fft(signal)
spectrum_array_length = n(len(signal_fft))

i = 0                                          # This is used instead of ' for i in range()'
while i < spectrum_array_length:               # because the 'signal' array can be rather long
if i >= spectrum_array_length/min_period :
signal_fft[i] = 0
i += 1

signal_filtered = scipy.ifft(signal_fft)

fft_filtered_signal = []

filtered_signal_length = n(len(signal_filtered))
i = 0
while i < filtered_signal_length:
fft_filtered_signal.append(float(signal_filtered[i]))
i += 1

norm = max(signal)/max(fft_filtered_signal)              # In case we need
# to renormalize
fft_filtered_signal_length = n(len(fft_filtered_signal)) # the filtered signal
i = 0                                                    # in order to obtain
while i < fft_filtered_signal_length:                    # the same amplitude
fft_filtered_signal[i] *= norm                       # as the inintial
i += 1                                               # signal had.

return fft_filtered_signal

more

You shouldn't filter by zeroing FFT bins. The frequency response *in between* bins will not be what you expect.

( 2012-08-27 17:25:14 -0600 )edit

Do you mean that some kind of smooth window function should be applied instead of this 'step-function' approach in order not to obtain the frequencies that were absent in the initial signal?

( 2012-11-28 05:39:00 -0600 )edit

Stats

Seen: 917 times

Last updated: Aug 19 '11