pythonscipybutterworth

modelling a cut-off frequency using scipy


I'm struggling to code a crossover filter for a loudspeaker. This function is meant to create a dummy input signal that can then be filtered, but when I change the frequencies and order when I call the function, it will no longer plot the filtered frequencies correctly on the graph.

On the scipy website, they say it's recommended to use second-order sections format when filtering, to avoid a numerical error with a transfer function's (ba) format but I am now unsure that's the best solution or if I am overlooking a simple error here.

see:

sosfilt: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfilt.html#scipy.signal.sosfilt

butter: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html#scipy.signal.butter

UPDATE - this appears to be happening for higher frequencies (I used a 3000 Hz cut-off) with both .butter, .sosfilt & .irrfilter, .lfilter commands

from scipy.signal import butter, sosfilt
import matplotlib.pyplot as plt
import numpy as np
def demo(f_s=1000,fc=15,fmin=10,fmax=20,n=10,nf=2):
    # demo filtered signal
    t = np.linspace(0,(15/fc),f_s,False) # sample size should give roughly ten peaks
    # generate signal
    inputf = np.linspace(fmin,fmax,nf) 
    # starts off with only the two frequencies,
    # can add more to check how the filter will respond
    sig = 0   
    for i in inputf:
        sig += np.sin(2*np.pi*i*t)
    # returns the filtered signal over a given time
    fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
    ax1.plot(t, sig)
    ax1.set_title('Initail sinusoidial signal from {} to {} Hz'.format(fmin,fmax))
    ax1.axis([0, (15/fc), -nf, nf])
    sosh = butter(n, fc, 'hp', fs=f_s, output='sos') #highpass
    sosl = butter(n, fc, 'lp', fs=f_s, output='sos') #lowpass
    filteredh = sosfilt(sosh, sig) #filtered high-pass
    filteredl = sosfilt(sosl, sig) #filtered low-pass
    ax2.plot(t, filteredh) # visualised high-pass signal
    ax2.plot(t, filteredl) # visualised low-pass signal
    ax2.set_title('After {} Hz Cut-off (high and low)-pass filter'.format(fc))
    ax2.axis([0, (15/fc), -nf, nf])
    ax2.set_xlabel('Time [seconds]')
    plt.tight_layout()
    plt.show()
demo()

This is the filter response the code needs to match for values that are not preset

This is what I'm currently getting


Solution

  • Tricky. I think you have confusion with the time base. f_s is supposed to be the sampling frequency, and you generate f_s samples, which would always be a full second. But when fc=3000, you only display the X axis as 0 to .005 seconds. Also, when fc=15, you generate f_s time samples running from 0 to 1. But when fc=3000, your time base will only run from 0 to 1/200, even though you're supposedly generating a full second of data.

    I changed the time base to

        t = np.linspace(0,1,f_s,False)
    

    and I think it's closer to what you were expecting.