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?
