[SOLVED] scipy butterworth, different orders for low and high frequencies

scipy butterworth, different orders for low and high frequencies

The original application of the butter function for bandpass filtering from scipy applies the same order (N) to lower and higher frequency cut-offs. any suggestions to apply different orders for low and high frequencies?

I thought of combing two separate low-cut and high-cut filters, but was not sure about the implementation.

Solution

import numpy as np
from scipy.signal import butter, filtfilt, freqz
import matplotlib.pyplot as plt


def butterworth_bandpass_filter(signal, lowcut, highcut, fs, order_low=4, order_high=4):
    """
    Apply a Butterworth bandpass filter to the signal and return the filtered signal and frequency response.

    Parameters:
    signal (numpy array): The input signal.
    lowcut (float): The low cutoff frequency in Hz.
    highcut (float): The high cutoff frequency in Hz.
    fs (float): The sampling frequency of the signal.
    order_low (int): The order of the low cutoff Butterworth filter.
    order_high (int): The order of the high cutoff Butterworth filter.

    Returns:
    filtered_signal (numpy array): The filtered signal.
    freqs (numpy array): The frequency bins of the filter response.
    response (numpy array): The frequency response of the filter.
    """
    # Normalize the frequencies by the Nyquist frequency
    nyquist = 0.5 * fs
    low = lowcut / nyquist
    high = highcut / nyquist

    # Design low cutoff filter
    b_low, a_low = butter(order_low, low, btype="high")
    # Design high cutoff filter
    b_high, a_high = butter(order_high, high, btype="low")

    # Apply the low cutoff filter
    low_filtered = filtfilt(b_low, a_low, signal)

    # Apply the high cutoff filter
    filtered_signal = filtfilt(b_high, a_high, low_filtered)

    # Calculate frequency response
    w, h_low = freqz(b_low, a_low, worN=1024)
    _, h_high = freqz(b_high, a_high, worN=1024)

    # Combined frequency response
    combined_response = h_low * h_high
    freqs = w * nyquist / np.pi  # Convert from normalized frequency to Hz

    return filtered_signal, freqs, combined_response, h_low, h_high


# Example usage
if __name__ == "__main__":
    # Create a sample signal: a mix of sine waves
    fs = 1000  # Sampling frequency
    t = np.linspace(0, 1, fs, endpoint=False)  # Time vector
    freq1 = 5  # Frequency of the first sine wave
    freq2 = 50  # Frequency of the second sine wave
    freq3 = 200  # Frequency of the third sine wave
    signal = (
        0.5 * np.sin(2 * np.pi * freq1 * t)
        + 0.5 * np.sin(2 * np.pi * freq2 * t)
        + 0.5 * np.sin(2 * np.pi * freq3 * t)
    )

    # Set filter parameters
    lowcut = 20  # Low cutoff frequency in Hz
    highcut = 100  # High cutoff frequency in Hz
    order_low = 7  # Order of the low cutoff filter
    order_high = 17  # Order of the high cutoff filter

    # Apply the bandpass filter
    filtered_signal, freqs, combined_response, h_low, h_high = (
        butterworth_bandpass_filter(signal, lowcut, highcut, fs, order_low, order_high)
    )

    # Plot the original and filtered signals
    plt.figure(figsize=(12, 6))
    plt.subplot(2, 1, 1)
    plt.plot(t, signal, label="Original Signal")
    plt.title("Original Signal")
    plt.xlabel("Time [s]")
    plt.ylabel("Amplitude")
    plt.legend()

    plt.subplot(2, 1, 2)
    plt.plot(t, filtered_signal, label="Filtered Signal", color="orange")
    plt.title("Filtered Signal")
    plt.xlabel("Time [s]")
    plt.ylabel("Amplitude")
    plt.legend()

    plt.tight_layout()
    plt.show()

    # Plot the combined frequency response
    plt.figure(figsize=(10, 5))
    plt.plot(freqs, np.abs(combined_response), label="Combined Frequency Response")
    plt.plot(freqs, np.abs(h_low), '--', label="L-cut Frequency Response")
    plt.plot(freqs, np.abs(h_high), '--', label="H-cut Frequency Response")
    plt.title("Combined Butterworth Bandpass Filter Frequency Response")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Gain")
    plt.grid()
    plt.legend()
    plt.xlim(0, fs / 2)  # Show only up to Nyquist frequency
    plt.show()

    # Export the frequency response to a CSV file
    np.savetxt(
        "butterworth_frequency_response.csv",
        np.column_stack((freqs, np.abs(combined_response))),
        delimiter=",",
        header="Frequency [Hz], Gain",
        comments="",
    )