psambit9791 / jdsp

A Java Library for Digital Signal Processing
https://jdsp.dev
MIT License
240 stars 45 forks source link

Result of Butterworth Bandpass and Bandstop Inconsistent with Python Scipy Library #76

Closed SnowCharmQ closed 3 months ago

SnowCharmQ commented 3 months ago

I used jdsp (Kotlin) for signal bandpass and bandstop processing. And I implemented a util function like the documentation as shown below:

    fun bandpass(
        signal: DoubleArray,
        sampleRate: Int,
        lowFreq: Double,
        highFreq: Double,
        order: Int
    ): DoubleArray {
        val nyquist = sampleRate / 2
        val low = lowFreq / nyquist
        val high = highFreq / nyquist
        val butterworth = Butterworth(sampleRate.toDouble())
        return butterworth.bandPassFilter(signal, order, low, high)
    }

However, when I used Python Scipy library for comparison, I found that their results were inconsistent. The Python Scipy implementation is shown below:

def bandpass(data, sampling_rate, low_freq, high_freq, order):
    nyquist = 0.5 * sampling_rate
    low = low_freq / nyquist
    high = high_freq / nyquist
    b, a = scipy.signal.butter(order, [low, high], btype='bandpass')
    data_bandpass = scipy.signal.lfilter(b, a, data)
    return data_bandpass

Can you give me an explanation or just let me know how to reproduce the result from Python using jdsp?

psambit9791 commented 3 months ago

@SnowCharmQ For JDSP, you do not need to divide the cut-off frequencies by the Nyquist frequency. Simply passing the cutoff frequencies to the function is sufficient.

So you need to use:

    fun bandpass(
        signal: DoubleArray,
        sampleRate: Int,
        lowFreq: Double,
        highFreq: Double,
        order: Int
    ): DoubleArray {
        val nyquist = sampleRate / 2
        val low = lowFreq
        val high = highFreq
        val butterworth = Butterworth(sampleRate.toDouble())
        return butterworth.bandPassFilter(signal, order, low, high)
    }
SnowCharmQ commented 3 months ago

@SnowCharmQ For JDSP, you do not need to divide the cut-off frequencies by the Nyquist frequency. Simply passing the cutoff frequencies to the function is sufficient.

So you need to use:

    fun bandpass(
        signal: DoubleArray,
        sampleRate: Int,
        lowFreq: Double,
        highFreq: Double,
        order: Int
    ): DoubleArray {
        val nyquist = sampleRate / 2
        val low = lowFreq
        val high = highFreq
        val butterworth = Butterworth(sampleRate.toDouble())
        return butterworth.bandPassFilter(signal, order, low, high)
    }

Thank you for your kind reply. I have already obtained the correct result.