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_signal
Asked: 2011-01-26 07:46:34 +0200
Seen: 1,367 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.
