Python High Pass Filter

Question

I implemented an high pass filter in python using this code:

from scipy.signal import butter, filtfilt
import numpy as np

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)
    return y

rawdata = np.loadtxt('sampleSignal.txt', skiprows=0)
signal = rawdata
fs = 100000.0

cutoff = 100
order = 6
conditioned_signal = butter_highpass_filter(signal, cutoff, fs, order)

I am applying this filter on a 100 kHz voltage signal and it works fine for cutoff frequencies >= 60 Hz. But it doesn't work below. I would like to cut off all the frequencies below 10 Hz. Any hints where my mistake is? What I observed is the lower the order of the filter the lower the cutoff frequency could be.

The sample Signal can be found here.

Answer 1

I add function to verify above Konstantin Purtov's code, so you could see relationship between order and cut off frequency. The code is mostly from Warren Weckesser.

import scipy
import sys  
from scipy import signal
from scipy import pi
from scipy.io.wavfile import write
import matplotlib.pyplot as plt
import numpy as np    
from scipy.signal import butter, lfilter, freqz   


# ----- ----- ----- -----    
def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y


# ----- -----    
# (1)
def foo(sel):
    if (sel == 1):
        # Filter requirements.
        order = 6
        fs = 300.0  # sample rate, Hz
        cutoff = 10  # desired cutoff frequency of the filter, Hz

        # Get the filter coefficients so we can check its frequency response.
        b, a = butter_highpass(cutoff, fs, order)

        # Plot the frequency response.
        w, h = freqz(b, a, worN=8000)
        plt.subplot(2, 1, 1)
        plt.plot(0.5 * fs * w / np.pi, np.abs(h), 'b')
        plt.plot(cutoff, 0.5 * np.sqrt(2), 'ko')
        plt.axvline(cutoff, color='k')
        plt.xlim(0, 0.5 * fs)
        plt.title("High Filter Frequency Response")
        plt.xlabel('Frequency [Hz]')
        plt.grid()

        # Demonstrate the use of the filter.
        # First make some data to be filtered.
        T = 0.5  # seconds
        n = int(T * fs)  # total number of samples
        t = np.linspace(0, T, n, endpoint=False)
        # "Noisy" data.  We want to recover the 20 Hz signal from this.
        data = np.sin(1.2 * 2 * np.pi * t) + 1.5 * np.cos(5 * 2 * np.pi * t) + 0.5 * np.sin(20.0 * 2 * np.pi * t)

        # Filter the data, and plot both the original and filtered signals.
        y = butter_highpass_filter(data, cutoff, fs, order)

        plt.subplot(2, 1, 2)
        plt.plot(t, data, 'b-', label='data')
        plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
        plt.xlabel('Time [sec]')
        plt.grid()
        plt.legend()

        plt.subplots_adjust(hspace=0.35)
        plt.show()
    else:
        print ('Please, choose among choices, thanks.')


# ----- -----
def main():
    sel = int (sys.argv[1])
    foo(sel)


# ----- ----- ----- ----- ----- -----
if __name__ == '__main__':
    main()

Answer 2

Since my reputation is low I am unable to comment on your question - "What is the relationship between cutoff and filter order?" This is not an answer to your original question.

For an FIR filter, for a given cutoff frequency, the slope of the impulse response plot (|H(f)| vs f) is steeper for a higher order filter. So, to achieve higher attenuation for the undesired frequency range, you increase the filter order. But what happens when the filter order is so high that the impulse response is an ideal box function? You will see an inter-symbol interference (ISI in digital communications) like effect. The intensity of this effect increases when the ratio of cutoff frequency to the sampling frequency gets smaller (think of the relation between the width of box function in frequency domain and the width of the main lobe of the sinc function).

I first observed this when I tried to implement a very narrow band low-pass IIR filter on a TI DSP microcontroller. The TI library, implemented the filter as a cascaded bi-quad structure to handle the well known truncation effects. This still did not solve the problem because the problem is not due to truncation alone. The way I solved this problem was to use an anti-aliasing filter followed by down-sampling the input signal, followed my desired low-pass IIR filter.

I understand that you are implementing a HPF, which is an LPF translated in frequency domain. Hope this answers some of your questions. Let me know if down-sampling works for you.

Answer 3

I hope this can help you:

import numpy as np
import pandas as pd
from scipy import signal
import matplotlib.pyplot as plt
def sine_generator(fs, sinefreq, duration):
    T = duration
    nsamples = fs * T
    w = 2. * np.pi * sinefreq
    t_sine = np.linspace(0, T, nsamples, endpoint=False)
    y_sine = np.sin(w * t_sine)
    result = pd.DataFrame({ 
        'data' : y_sine} ,index=t_sine)
    return result

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y

fps = 30
sine_fq = 10 #Hz
duration = 10 #seconds
sine_5Hz = sine_generator(fps,sine_fq,duration)
sine_fq = 1 #Hz
duration = 10 #seconds
sine_1Hz = sine_generator(fps,sine_fq,duration)

sine = sine_5Hz + sine_1Hz

filtered_sine = butter_highpass_filter(sine.data,10,fps)

plt.figure(figsize=(20,10))
plt.subplot(211)
plt.plot(range(len(sine)),sine)
plt.title('generated signal')
plt.subplot(212)
plt.plot(range(len(filtered_sine)),filtered_sine)
plt.title('filtered signal')
plt.show()