I recently translated the code from
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:
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.