berndporr / iirj

An efficient IIR filter library written in JAVA
Apache License 2.0
134 stars 36 forks source link

Question on Butterworth highPass, cutoff==Nyquist, coefficients are NaN? #16

Closed mickvf closed 3 years ago

mickvf commented 4 years ago

Hi and thanks for this wonderful library - I am learning a lot!

I wanted to run this by you to see if you had any thoughts. Perhaps I am doing something wrong - would not be the first time in history =). I am trying to filter seismic time series data so we generally work in the sub-100 Hz realm. Much of the time, 20 samples per sec.

For a low-pass Butterworth filter order 4, sampling rate 20/sec and cutoff frequencies 1 through 10, all is fine.

For the high-pass case however, I run into trouble right at the Nyquist barrier with cutoff = 10.

Stepping thru the call to filter(), with a breakpoint in the Cascade class, ~line 70, I can see that the m_biquads array (index 0) has NaN for the 3 "b" coefficients (m_b0, m_b1, m_b2). This leads to all NaNs for my data after the filter call.

So then I had a go at stepping into the initialization of the filter:

butterworth.highPass(4, rate, freqCutoff); // 4, 20.0, 10.0

and to me things seem OK in highPass() and setupHighPass(). I did manage to step a little further, into the Biquad initialization. But in the setCoefficients() method, I don't detect any NaNs or div by zero or anything really obvious a mortal like me might spot.

So finally my question to you: Are there particular cases where I might need to make additional checks of my web service user's inputs in which I can predict that things are going to get unstable, or catch the exception and return no data?

Thanks for your patience and have a wonderful day! Mick Van Fossen Seattle

mickvf commented 4 years ago

Hi again. With some further thought and development, I've decided that the onus is on me to sense when NaNs are coming back: I think it's mainly an artifact of us working mainly with such low sampling rates (20), and in any case, our users aren't expected to be working so close to the Nyquist frequency. I'm returning a flat line timeseries in these edge cases plus the necessary information for them to figure out what's happening.

Once again thank for this project!