I recently translated the code from scipy.signal.bessel to C in order to implement a bessel filter in a piece of analysis software I am working on. Everything went well, I was able to reproduce exactly the filter coefficients, and the filter has the right cutoff and produces numerically exact (ie, my C implementation and the direct use of scipy.signal.filtfilt agree numerically) results up until the cutoff frequency. After the cutoff frequency, though, the C implementation levels off. The image below shows the comparison of the PSD after the two filters (200kHz cutoff, 8th order, forward then back filtered) are applied to the same signal (x - frequency in $Hz$ - y - amplitude in $pA^2/Hz$). Blue is my C code, green is scipy:

enter image description here

My first thought was that it might be numerical noise. I use 64-bit doubles for everything in the C code, and my first guess was that maybe the python implementation uses 128-bit doubles, since the point at which the python code levels off is suspiciously close to 1/2^106. But the C code levels off well above 1/2^53, and none of my digging through scipy code suggests that they are using long doubles.

Another possiblity is that in my C code implementation, I have not implemented Gustaffson's method to ensure that the forward and backward filtered results are the same, which is done in the python version. Could that be the source of the higher noise in high frequencies?

Any ideas what is happening? I am fairly confident that the filter is implemented correctly since it gives correct results up to the cutoff frequency.

Numerically it hardly matters - leveling off at 1e-10 vs 1e-30 is irrelevent to the actual signals. I would just like to understand where the difference comes from.

  • $\begingroup$ That does look suspiciously like an intermediate value used in the C filtering is being stored in a lower-bit-count variable. This would only show up in the stop band of the filter, though. $\endgroup$ – Peter K. Apr 29 '16 at 17:10
  • $\begingroup$ You can see my implementation here, if it helps: github.com/shadowk29/CUSUM/blob/master/bessel.c. The actual filtering is done here: github.com/shadowk29/CUSUM/blob/master/detector.c in filter_signal(). Could you explain why numerical roundoff error would only show up in the stop band? $\endgroup$ – KBriggs Apr 30 '16 at 15:53
  • $\begingroup$ Is that a frequency response plot? It doesn't look like a Bessel filter $\endgroup$ – endolith May 5 '17 at 15:12
  • 1
    $\begingroup$ @endolith, no that's a PSD of the signal off one of my sensors. I was just using it to compare the two different implementations. I found the error - the least significant bit of the input to the second case was being flipped, so it was a bug unrelated to the filter implementation. $\endgroup$ – KBriggs Jun 9 '17 at 20:15

Your Answer

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

Browse other questions tagged or ask your own question.