3
$\begingroup$

I have implemented Python signal processing code using a butter bandpass filter and order 3. It is working pretty well.

I am required to implement the same signal processing in C++, so I replicated the filtering logic in C++.

  • Python Parameters:
    • Order: 3
    • Type: bandpass
  • C++ Parameters
    • Order: 2
    • Type: bandpass

Coefficient (a, b) - I generate coefficients using b, a = butter(3, [low_frequency / fs * 2, high_frequency / fs * 2], "bandpass"), and using the same a, b values in Python and C++.

These parameters were working fine for the majority of the data.

Recently for a specific dataset, I am getting 25% less positive predictivity in the C++ signal processed data for a new dataset I am testing. I tried to change the order from 1,2,3,4 in C++, but it's not improving the accuracy of signal testing.

Here is the code snippet for both.

Python

def butter_bandpass_filter(data, lowcut, highcut, fs, order=3, mode=0, verbose=0):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = signal.butter(order, [low, high], btype='band')

    if verbose:
        for tap in a:
            print("%.20f," % tap)
        for tap in b:
            print("%.20f," % tap)

        import pylab as plt
        plt.figure()
        plt.clf()
        w, h = signal.freqz(b, a, worN=8000)
        plt.plot((w / np.pi) * fs/2.0, np.absolute(h), linewidth=2)
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain')
        plt.title('Frequency Response')
        plt.ylim(-0.05, 1.05)
        plt.show()

    if mode == 0:
        y = signal.filtfilt(b, a, data)

x = butter_bandpass_filter(x, 1, 25, 300, order=3, verbose=0)

C++

void basic_feature_extraction(const double *data, const int N, double * basic_feature)
{
    double *fdata = new double[N];
    double *pan_tompkins_array = new double[N];
    // 2nd order filter [1,25] Hz bandpass
    ButterFilter myButter;

    myButter.myfilter.len = 5;
    myButter.myfilter.b = new double[myButter.myfilter.len];
    myButter.myfilter.b[0] = 0.0461318020933179;
    myButter.myfilter.b[1] = 0;
    myButter.myfilter.b[2] = -0.0922636041866359;
    myButter.myfilter.b[3] = 0;
    myButter.myfilter.b[4] = 0.0461318020933179;

    myButter.myfilter.a = new double[myButter.myfilter.len];
    myButter.myfilter.a[0] = 1;
    myButter.myfilter.a[1] = -3.28877608589248;
    myButter.myfilter.a[2] = 4.07514013666218;
    myButter.myfilter.a[3] = -2.27808862062307;
    myButter.myfilter.a[4] = 0.491812237222575;

    myButter.BiDirectionalFilter(data, fdata, N);

void ButterFilter::FiltData(const double * sig, double * sout,const int data_len)
{
    int N = myfilter.len;
    int nfact = 50 * (N - 1);
    if (nfact>data_len)
    {
        nfact = data_len - 1;
    }
    //�ź�����
    double* s_ext = new double[data_len + nfact];   // ���غ�������ź�
    double* tmps = new double[data_len + nfact];    // ��Ϊ��������е�����ź�
    int data_len_ext = signalExtend(sig, s_ext, data_len, nfact);
    for (int i = 0; i < data_len_ext; i++)
    {
        tmps[i] = s_ext[i];
    }
    for (int i = N; i < data_len_ext; i++)
    {
        double tmp = 0;
        tmp += myfilter.b[0] * s_ext[i];
        for (int j = 1; j < N; j++)
        {
            tmp += myfilter.b[j] * s_ext[i - j];
            tmp -= myfilter.a[j] * tmps[i - j];
        }
        tmps[i] = tmp;
        //printf("out[%d]=%f\n",i, tmp);
    }
    //���
    for (int i = 0; i < data_len; i++)
    {
        sout[i] = tmps[i + nfact];
    }
    delete[]tmps;
    delete[]s_ext;
}
void ButterFilter::BiDirectionalFilter(const double * sig, double * sout, const int data_len)
{
    double* tmp1 = new double[data_len];
    double* tmp2 = new double[data_len];
    FiltData(sig, tmp1, data_len);
    dataRollingOver(tmp1, tmp2, data_len);
    FiltData(tmp2, tmp1, data_len);
    dataRollingOver(tmp1, sout, data_len);
    delete[]tmp1;
    delete[]tmp2;
}
void ButterFilter::dataRollingOver(const double * sig,double * sout,const int data_len)
{
    for (int i = 0; i < data_len; i++)
    {
        sout[data_len - i - 1] = sig[i];
    }
}

How can I get C++ and Python to both give the same results?

$\endgroup$
4
  • 1
    $\begingroup$ but a 2nd-order Butterworth filter is not the same as a third-order Butterworth filter. These will, very much by design, never give the same results. $\endgroup$ May 23 at 10:36
  • $\begingroup$ I agee @MarcusMüller: This different order is not giving so I tried with the same order. And other order values eg 1,2 but still getting bad result. $\endgroup$ May 23 at 10:37
  • $\begingroup$ "butter" = Butterworth filter? $\endgroup$ May 24 at 0:36
  • 1
    $\begingroup$ What is "signal.butter"? What is the import for 'signal'? What is the Python library? Is it an interface to MATLAB? There is "signal", but it has nothing to do with signal processing in the sense here. $\endgroup$ May 24 at 0:36

1 Answer 1

10
$\begingroup$

Not really an answer but too long for a comment. Standard software development and debug techniques should solve this

  1. Don't use code on real data before it has passed all unit tests
  2. First write the unit tests (which are the requirements for your application) and then start writing the actual code
  3. If you fail any unit tests, debug it step by step. Start with a unit impulse as signal and impulse response and than gradually increase complexity. Only use test cases with known answers (either by calculating them by hand or use a "known good" reference"
$\endgroup$
3
  • 3
    $\begingroup$ Excellent recommendation. I think this is a proper answer, actually! If you know what the impulse response should look like, at least for a finite amount of impulse (which OP does, having a reference), it's easy to check against that. $\endgroup$ May 23 at 12:21
  • $\begingroup$ Let me try and update you guys. Thanks for assistance $\endgroup$ May 23 at 13:21
  • 1
    $\begingroup$ Also there is no need to re-invent the wheel. There are plenty of C++ IIR implementations out there $\endgroup$
    – Hilmar
    May 23 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.