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.
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..)
Finally I decided to create a simple as a log hand-made FFT filter. The concept is as follows:
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_signalYou shouldn't filter by zeroing FFT bins. The frequency response *in between* bins will not be what you expect.
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?
Asked: 14 years ago
Seen: 1,543 times
Last updated: Aug 19 '11
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
