Revision history [back]

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