Ask Your Question
1

Low-pass filter

asked 2011-01-26 00:46:34 -0600

v_2e gravatar image

updated 2015-01-14 07:55:58 -0600

FrédéricC gravatar image

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 flag offensive close merge delete

2 answers

Sort by » oldest newest most voted
3

answered 2011-01-26 01:37:41 -0600

DSM gravatar image

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
# see thread http://mail.scipy.org/pipermail/scipy-user/2009-September/022565.html

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

edit flag offensive delete link more
0

answered 2011-08-19 03:37:56 -0600

v_2e gravatar image

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
edit flag offensive delete link more

Comments

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

endolith gravatar imageendolith ( 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?

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

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-01-26 00:46:34 -0600

Seen: 799 times

Last updated: Aug 19 '11