ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 28 Nov 2012 12:39:00 +0100Low-pass filterhttps://ask.sagemath.org/question/7865/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.Wed, 26 Jan 2011 07:46:34 +0100https://ask.sagemath.org/question/7865/low-pass-filter/Answer by DSM for <p>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? <br/>
I have read that FFT could be used for this, but haven't got the whole idea.</p>
<p>Thanks.</p>
https://ask.sagemath.org/question/7865/low-pass-filter/?answer=12012#post-id-12012Yes, 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.
<pre><code>
# 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()
</code></pre>
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..)
Wed, 26 Jan 2011 08:37:41 +0100https://ask.sagemath.org/question/7865/low-pass-filter/?answer=12012#post-id-12012Answer by v_2e for <p>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? <br/>
I have read that FFT could be used for this, but haven't got the whole idea.</p>
<p>Thanks.</p>
https://ask.sagemath.org/question/7865/low-pass-filter/?answer=12591#post-id-12591 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](http://www.swharden.com/blog/2009-01-21-signal-filtering-with-python/) 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_signalFri, 19 Aug 2011 10:37:56 +0200https://ask.sagemath.org/question/7865/low-pass-filter/?answer=12591#post-id-12591Comment by endolith for <p>Finally I decided to create a simple as a log hand-made FFT filter. The concept is as follows:</p>
<ul>
<li>Perform the direct Fourier transform.</li>
<li>Zero out the coefficients corresponding to the target frequencies</li>
<li>Perform the inverse Fourier transform.</li>
</ul>
<p>I borrowed this idea described <a href="http://www.swharden.com/blog/2009-01-21-signal-filtering-with-python/">here</a> and wrote the following function:</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/7865/low-pass-filter/?comment=19158#post-id-19158You shouldn't filter by zeroing FFT bins. The frequency response *in between* bins will not be what you expect.Tue, 28 Aug 2012 00:25:14 +0200https://ask.sagemath.org/question/7865/low-pass-filter/?comment=19158#post-id-19158Comment by v_2e for <p>Finally I decided to create a simple as a log hand-made FFT filter. The concept is as follows:</p>
<ul>
<li>Perform the direct Fourier transform.</li>
<li>Zero out the coefficients corresponding to the target frequencies</li>
<li>Perform the inverse Fourier transform.</li>
</ul>
<p>I borrowed this idea described <a href="http://www.swharden.com/blog/2009-01-21-signal-filtering-with-python/">here</a> and wrote the following function:</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/7865/low-pass-filter/?comment=18642#post-id-18642Do 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?Wed, 28 Nov 2012 12:39:00 +0100https://ask.sagemath.org/question/7865/low-pass-filter/?comment=18642#post-id-18642