15

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.

Patrick Roberts
  • 49,224
  • 10
  • 102
  • 153
ChrisG
  • 333
  • 1
  • 3
  • 10
  • Where do you set your passband? Do you have your formula that you can paste? and your `filtfilt()` function. and your `butter()` function – Daniel Lee Aug 19 '16 at 06:39
  • @Daniel Lee: I edited my question. The butter() and filtfilt() functions are from scipy.signal. – ChrisG Aug 19 '16 at 06:44
  • Can you give some examples how you use it, what the observed behavior is like and what the expected behavior would be? Questions seeking debugging help require a clear problem statement. _works fine here_ vs. _doesn't work there_ is not sufficiently clear. – moooeeeep Aug 19 '16 at 06:50
  • @moooeeeep: I added a sample signal and an example how I apply the filter. – ChrisG Aug 19 '16 at 07:07

3 Answers3

27

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()
Konstantin Purtov
  • 799
  • 1
  • 11
  • 19
  • Thanks for your help! Now I am a little confused. What is the relationship bewteen cutoff and order? – ChrisG Aug 19 '16 at 07:10
  • There is no relationship. The "cutoff" means the cut frequency. The "order" determine the accuracy of filter. – Konstantin Purtov Aug 19 '16 at 07:20
  • Thats what I thougt. But if you set the order to '30'. The result changes strongly. And for high orders, the low cutoff frequency do not work. – ChrisG Aug 19 '16 at 07:22
  • It depends on many factors. You can design and try what you want with Matlab filter design toolkit. As I think the LMS filters is the best way for removing the very low frequencies (http://stackoverflow.com/questions/18252160/how-to-apply-an-adaptive-filter-in-python) – Konstantin Purtov Aug 19 '16 at 08:03
  • Hi. Do you have any idea how can i use an mp3 file instead of the sine wave. I want to remove all the frequencies below certain threshold from the mp3 and wav file – It's a trap Oct 31 '17 at 09:24
  • @It'satrap Hi. Yep. You can read mp3 file as a signal, use filter and save it back. The good example to read mp3 and use filter is here https://stackoverflow.com/questions/15311853/plot-spectogram-from-mp3 – Konstantin Purtov Nov 01 '17 at 16:08
9

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.

Sriharsha Madala
  • 275
  • 3
  • 10
3

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()
Cloud Cho
  • 1,594
  • 19
  • 22