Band-pass filter help?

[clepsydrae]

Hi -- disclaimer: I know next to nothing about filters.

I have implemented what I understand to be a BiQuad filter in JSFX. It works fine with complex signals, AFAICT. (Meaning, it filters as it should).

But I have noticed that if I pass in a 100Hz sine and set the filter center frequency to 105Hz with a narrow width, the resulting outputs oscillates in amplitude, as if going in and out of phase.

If I pull up a ReaEQ instead and use a band pass in the same arrangement, there is no oscillation, just a lower overall signal, as I would expect. -- nope, see below

Is this due to the nature of a BiQuad filter, or does it point to a mistake in my implementation?

Thanks for any help!

[Time Waster]

I no filter expert either, but is article seems to suggest that biquads may not be the best at low frequencies: http://www.earlevel.com/main/2003/02/28/biquads/

[clepsydrae]

Thanks -- do you suspect the "instability" that article mentions is what is causing the oscillation I'm seeing? They mention the digital state variable filter as an alternative, which I'll check out.

I should mention that I'm using my BiQuad as a band pass, and I believe I implemented per this link (it was a while ago when I coded it up.) I think I used the "BW" method, and the second set of BPF coefficients.

[sai'ke]

I can't be sure from the description, but I suspect something is off in the implementation. By instability they usually mean instability in the sense of bounded input leads to bounded output. If it was really unstable, it would simply put the signal to infinity or -infinity (unless you are also clipping it somehow).

Do the other properties hold? Do you see the correct n dB/octave behavior? And the peak? Try testing it with white noise, you can nicely see the curve then.

In case you're looking for those state variable filters, there's this document that has nice zero feedback delay versions of an SVF:
https://cytomic.com/files/dsp/SvfLin...Optimised2.pdf

IHTH.

[clepsydrae]

Thanks Sai'ke!

I just discovered that I was wrong about ReaEQ -- it does the same thing. So maybe it's normal after all... Is this just unavoidable filter "ringing" perhaps?

If you place that demo code (or ReaEQ) on a track, put a 100Hz tone in front of it, set the band pass filter center to 105Hz and the octave width to .005 (in ReaEQ the bandwidth to .01), the output looks like this (this image from the JSFX, but similar pattern for ReaEQ):



...which represents a 5Hz oscillation, i.e. the difference between 105 and 100. It stabilizes out over time, which made me suspect ringing.

I just went over my implementation, and I can't find any mistakes in my code... in case anyone wants to check it out, this is a demo:

slider1:300<0,22050,1>filt cent
slider2:.8<0,1,.01>octave width

@init
recidx=0;
recsize=srate;
filtbuf=recidx+recsize;
outbuf=filtbuf+recsize;

@sample

function sinh(myval)
( ($e^myval - $e^(-myval))/2; );

function calcFilterVars(myfreq)
(
  // 0 to 1, 1 is wide, 0 is narrow
  octavewidth=slider2;
  
  w0=2*$pi*myfreq/srate;
  alpha=sin(w0)*sinh( log(2)/2 * octavewidth * w0/sin(w0) );
  b0=alpha;
  //b1=0;
  b2=-alpha;
  a0=1+alpha;
  a1=-2*cos(w0);
  a2=1-alpha;
);

function applyBiQuadFilter(centfreq, mybuf, mynumsamples, mydest)
(
  calcFilterVars(centfreq);

  mydest[0]=mydest[1]=0;
  sampidx=2;
  loop(mynumsamples-2,
    mydest[sampidx]=(b0/a0)*mybuf[sampidx] + (b2/a0)*mybuf[sampidx-2] -
      (a1/a0)*mydest[sampidx-1] - (a2/a0)*mydest[sampidx-2]; // b1 term is zero for BPF
    sampidx+=1;
  );
);

recbuf[recidx]=spl0;
spl1=spl0=outbuf[recidx];
recidx=(recidx+1)%recsize;

recidx==0?
(
  applyBiQuadFilter(slider1, recbuf, recsize, filtbuf);
  memcpy(outbuf,filtbuf,recsize);
);

(NOTE: since that implementation just crudely processes blocks at a time, and since the filter as written skips the first two terms, there is a hiccup every srate samples causing a jump in a spectral meter, especially if you're looking at a pure tone.)

From the cookbook:

y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]
                        - (a1/a0)*y[n-1] - (a2/a0)*y[n-2] 

alpha = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) )
b0 =   alpha
b1 =   0
b2 =  -alpha
a0 =   1 + alpha
a1 =  -2*cos(w0)
a2 =   1 - alpha

[bezusheist]

What do you expect when you crank the resonance ?

[sai'ke]

Ah, you mean the transient response? For a moment I thought you meant the wobble was persistent.

Having some transient response is normal for any filter. The response of a linear filter to a sine that starts is always a transient response plus a steady state response. The steady state response is the one it eventually evolves into. In this case the transient response is quite long because you select such a narrow band. For wider bands the transient response becomes shorter.

Your code doesn't seem to work on my end though. If I find some time, I'll have a closer look tomorrow.

[clepsydrae]

Originally Posted by sai'ke

Having some transient response is normal for any filter.

Ah, ok, thanks. So now my quest is to find a very narrow bandpass filter with minimal transient ringing. I don't (think I) care about response ripple in the either the stop band or pass band, really, I just want to roughly isolate the content in the pass band without temporal ringing. Any tips on filter choice?

I poked around with some other EQs trying to find something that doesn't ring as much... ReaFIR was better than the BiQuad, but still has some wobble, even when the pass band is wider and the stop bands aren't as drastically low. This makes me think an fft/tweak/ifft may not be a great choice.

I had my best success using iZotope Alloy 2 and setting up a brickwall HPF and LPF around the target frequency. This exhibited a ring but it stopped very quickly (~200-300ms). The Alloy 2 filters are extremely steep; apparently "elliptic" filters, according to their documentation. Those seem a little more intense to code than I want to get in to, though, yikes.

I was eyeballing Chebyshev type I. I could try to implement that but I wonder if I should expect the same lack of ring as the Alloy 2 elliptic brickwall?

Your code doesn't seem to work on my end though. If I find some time, I'll have a closer look tomorrow.

Ooops, sorry, it was dependent on pin routing. Fixed now, in case anyone is still curious. (Now takes L as input, outputs left to L, R).

Originally Posted by bezusheist

What do you expect when you crank the resonance ?

Knowing almost nothing about filters, I had no idea what to expect, hence my question. :-)

[sai'ke]

With such a narrow bandwidth it sounds more like an engineering problem than a music problem. It would probably help if you told us what you are designing this for and what the constraints are, so that we would not miss a solution that might not be a linear filter.

Elliptic filters ring a lot, so I wouldn't use that. What you want though is a bit of a contradiction. Typically, if you restrict yourself to real-time causal linear filters (cannot look into the future) then narrow band = longer transient.

One option to reduce the transient could also be to find better initial conditions for the filter, but that would require knowing a lot about what problem you're trying to solve.

Are you going to be changing the input / cutoff frequency? Is the bandwidth fixed? Does it have to be real-time and causal? Is a group delay allowed?

[clepsydrae]

Thanks for the consideration -- I'm toying with a simple pitch-detection algorithm. Imagine a guitar tuner that knows it is tuning the note "E". I want to filter everything out except a band around that frequency, then schmitt trigger on the result.

It works great when the biquad is relatively wide, but when I narrow it works much better, except that the resulting oscillating amplitude causes the pitch readings from the schmitt trigger to wobble a little (which makes sense, given how schmitt works.) It's clear from my tests that if I had a non-wobbling result, the schmitt would be even more accurate. (And if the frequency coming in perfectly matches the band pass center, the tight filter resonance causes the result to grow asymptotically rather than wobble, which also confuses the schmitt.)

(I'm aware that I could use autocorrelation as well, which I'm exploring.)

Originally Posted by sai'ke

Elliptic filters ring a lot

Huh, that's interesting. Does that imply that iZotope is doing something special behind the scenes besides a generic elliptic? (...asking because when I used the HP/LP combination there was so little ringing apparent.)

Until you said that, my latest scheme was to use SciPy to generate the coefficients of an n-pole elliptic, and then implement that in JSFX via Direct-Form I (I've been reading up a little, so I may be off in the weeds here.) I might do this anyway, just for the educational sake of it. :-)

Are you going to be changing the input / cutoff frequency?

Yes -- I'll need to move the filter band around all the time. But whatever filter used should have no effect no the resulting frequency that is passed through.

Is the bandwidth fixed?

Ideally I could scale the bandwidth at will. Even if it was fixed, I would need to tune it at time of filter design.

Does it have to be real-time and causal? Is a group delay allowed?

Real-time/causal: no, as long as processing takes less than, say, 1 second.

Group delay: I think it's allowed as well, though I profess a weak understanding.

Basically, I don't care what happens to the signal or when it happens, as long as the pitch isn't messed with and the band is tight and the ringing is minimized. (That's why I also said that ripple in the pass band is fine.)

Thanks for any ideas!

[sai'ke]

I don't know what it does. It's possible that it does something clever with the initials to reduce the ringing. I don't know, without documentation it's guesswork. Are you sure you had the same bandwidth? Any filter with a tiny bandwidth will be slow to respond.

If you're set on the filter approach and you really just want to detect one band, you could consider a Goertzel filter to do the complete tone detection (this is what they use to detect tones from phone keypads). It's kind of like an FFT for a single bin written as a filter. It has a very controllable width, but you'd need a lot of samples for a high resolution (same problem really).
https://www.embedded.com/design/conf...tzel-Algorithm

Issue is though, for any real tuning, it'd probably not be a perfect pitch and it might be a lot easier to tune if the pitch is shown, rather than trying to hit some pitch.

Could also have a look at the cepstrum. The inverse FFT of the log of the magnitude of an FFT of a short window. It's probably a bit more robust than single bin picking since it also uses information from the harmonics.

Something like this would also seem worth playing with (time domain route):
http://miracle.otago.ac.nz/tartini/p...Find_Pitch.pdf

[ErBird]

I'm finding your code/implementation very peculiar. Have you taken a look at the RBJ Highpass/Lowpass that comes with Reaper? Certainly you don't need all this complexity for a simple biquad bandpass fiter.

[clepsydrae]

EDIT: you can probably skip this post and see below -- I think I found a workable solution...

------------------

Thanks for the ideas!

So I've been playing with SciPy's iirfilter and I'm able to design filters and put the coefficients in to JSFX and they work generally as expected.

But sometimes the results are unstable (output going to infinity). Maybe this is because I have no way to get the precise coefficients from python to JSFX? I'm doing this:

from scipy import signal
b, a = signal.iirfilter(8, 100, rp=5, rs=60, btype='highpass', ftype='ellip', output='ba', fs=44100)
for i in range(len(a)): print('a[' + str(i) + '] = {:.100f}'.format(a[i]) + ';')
for i in range(len(b)): print('b[' + str(i) + '] = {:.100f}'.format(b[i]) + ';')

Which generates e.g.:

a[0] = 1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000;
a[1] = -7.9643984992593175320507725700736045837402343750000000000000000000000000000000000000000000000000000000;
a[2] = etc...

...which I copy/paste into my testing JSFX. For lower-order filters, or for mid-range frequencies, it seems to work fine, but it gets unstable pretty easily at e.g. 100Hz or orders above 5, give or take.

I know that JSFX isn't using most of those digits when it's interpreting that python output. Is that the likely cause of the instability? SciPy outputs float64's, so I had my fingers crossed that the decimal representation might come into JSFX intact.

If I go to the trouble of computing those coefficients in JSFX, should I expect better (meaning stable) results? Or are such filters inherently unstable somehow? Or maybe the config parameters I'm providing are unreasonable?

For the record (or in case anyone else wants to use it), the JSFX test code is this:

@init

filtlen=9;
filtinbuf=0;
filtoutbuf=filtinbuf+filtlen;
a=filtoutbuf+filtlen;
b=a+filtlen;

filtidx=0;

// 100Hz highpass elliptic, 8th order, 60dB stopband rejection, 5dB passband ripple

a[0] = 1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000;
a[1] = -7.9643984992593175320507725700736045837402343750000000000000000000000000000000000000000000000000000000;
a[2] = 27.7534825146954347019345732405781745910644531250000000000000000000000000000000000000000000000000000000;
a[3] = -55.2684920460229136551788542419672012329101562500000000000000000000000000000000000000000000000000000000;
a[4] = 68.7941716398878355676060891710221767425537109375000000000000000000000000000000000000000000000000000000;
a[5] = -54.8074689089812565612191974651068449020385742187500000000000000000000000000000000000000000000000000000;
a[6] = 27.2924281122364398299851018236950039863586425781250000000000000000000000000000000000000000000000000000;
a[7] = -7.7667814203223066726877732435241341590881347656250000000000000000000000000000000000000000000000000000;
a[8] = 0.9670586077661216251044606906361877918243408203125000000000000000000000000000000000000000000000000000;

b[0] = 0.5527712070143768752217283690697513520717620849609375000000000000000000000000000000000000000000000000;
b[1] = -4.4219146405610425532017870864365249872207641601562500000000000000000000000000000000000000000000000000;
b[2] = 15.4760637427930820564370151259936392307281494140625000000000000000000000000000000000000000000000000000;
b[3] = -30.9513625183507130600446544121950864791870117187500000000000000000000000000000000000000000000000000000;
b[4] = 38.6888844182085875900156679563224315643310546875000000000000000000000000000000000000000000000000000000;
b[5] = -30.9513625183507130600446544121950864791870117187500000000000000000000000000000000000000000000000000000;
b[6] = 15.4760637427930838327938545262441039085388183593750000000000000000000000000000000000000000000000000000;
b[7] = -4.4219146405610425532017870864365249872207641601562500000000000000000000000000000000000000000000000000;
b[8] = 0.5527712070143768752217283690697513520717620849609375000000000000000000000000000000000000000000000000;

@sample

// see https://stackoverflow.com/questions/16879642/how-to-implement-chebyshev-type-2-lpf-in-java
newx=(spl0+spl1)/2;
newy=newx*b[0];
loopidx=0;
loop(filtlen-1,
  newy +=  filtinbuf[(filtidx+loopidx)%filtlen]*b[(loopidx+1)%filtlen] -
          filtoutbuf[(filtidx+loopidx)%filtlen]*a[(loopidx+1)%filtlen];
  loopidx+=1;
);

filtidx-=1;
filtidx<0?filtidx=filtlen-1;
filtinbuf[filtidx]=newx;
filtoutbuf[filtidx]=newy;
spl0=newy;
spl1=newy;

[clepsydrae]

Originally Posted by sai'ke

Are you sure you had the same bandwidth? Any filter with a tiny bandwidth will be slow to respond.

When I used the elliptic HP/LP combo I may have had a wider bandwidth to the equivalent biquad possibly because the slopes were so much steeper and I may have been able to get away with wider separation as a result.

[clepsydrae]

Hey, so I think I found something that will work... a good old fashioned fft/ifft with ~10Hz on either side of the estimated frequency as the pass band, and the rest zeroed out. (I had passed on that because of what I saw in ReaFIR in testing, but reconsidered in light of my last post.)

In some initial tests using a Hamming window it's pretty close to ideal for my purposes. It starts to freak out when the estimated frequency is too far from the edges of the fft bin frequencies, but when it's close(ish) to on-center the amplitude doesn't wobble much, which is great. And I have a separate idea for improving the estimated frequency position to prevent the off-center case, so I think I'm good!

Anyhow I won't bore you with the details. I really appreciate all the help and education!

[sai'ke]

If that's what you've decided on, it's probably worthwhile to have a look at the Goertzel filter, which allows you to calculate just that one bin more cheaply.

[flicker177]

But I have noticed that if I pass in a 100Hz sine and set the filter center frequency to 105Hz with a narrow width, the resulting outputs oscillates in amplitude, as if going in and out of phase.

If you're looking at your low frequency filter with an FFT based analyzer you'll see bouncing and instability. FFT's don't work well at low frequencies because of the very limited resolution down there.

I made a sweep generator/plotter as jsfx effects for exactly that reason when I was developing filters for a dbx decoder, also written in jsfx. You'll find that in the Reaper Stash in the "js effects" section. Search for "sweeper" and you'll find it. It comes with a Reaper project that has the routing all set up so you can just plug your filter in as an effect and plot it with good resolution. If its *very* sharp you need to run the sweep slowly.

[edit]
Out of curiosity I ran your filter here (the chebychev). I must be missing something, I just used your code as a jsfx and when I try to plot it it totally overloads, even with the input signal set at -120dB. Dunno why...
[/edit]
--Bill

[jcjr]

Hi clepsydrae. Maybe already mentioned somewhere in the thread, but the ringing related to high-selectivity filtering is an aspect of the Heisenberg Uncertainty Principle.

In quantum mechanics it means as we make increasingly accurate measurement of perhaps particle velocity, our knowledge of the particle's position becomes increasingly less accurate. And vice-versa.

In audio analysis, the more accurately we measure frequency, the less accurately we can know the time at which the frequency occurs, and vice-versa.

For instance a dirac delta function is a very brief "spike" for instance in digital audio a signal that is all zero-value samples except one lonely sample with a value of 1.0. In time domain, we can "exactly" locate the timestamp of this single value of 1.0 in the sea of 0.0 values.

But if we take the fourier transform of the dirac delta, we measure "white noise" all frequency harmonics of equal amplitude.

OK, if you do the equivalent of a "dirac delta" in frequency domain-- Set one harmonic to 1.0 and all other harmonics to 0.0. When we do the inverse fourier transform to get time domain, we get a pure sine wave of a precisely-known frequency, but it stretches in duration forever.

So we can precisely know the frequency but nothing about the event time, or we can precisely know the event time but nothing about the frequency. Or we can know any "mix" of the two. If we want better time-discrimination (less ringing) we have to put up with less frequency-resolution. And vice-versa. The engineering is finding a tradeoff between the two that works best to solve whatever problem is at hand.

[clepsydrae]

Thanks @flicker177 for the ideas and the plotter tool -- I will check that out. I was examining the wobble in the recorded waveform, not the fft (see third post). The second JSFX I posted was to demonstrate that the filter runs away, so it's working as intended. :-) I was trying to understand if it was just inherent to the filter or because I was copy/pasting coefficients from python.

And thanks @jcjr for an intuitive explanation; I still only half-grasp the math going on here, but I think I understand the limitations of the technique. The FFT/iFFT method is convenient for me because I can adjust bandwidth and so forth without needing to design filter coefficients, even if I now understand there is no free lunch in terms of the ringing, etc.

Have you taken a look at the RBJ Highpass/Lowpass that comes with Reaper? Certainly you don't need all this complexity for a simple biquad bandpass fiter.

Thanks -- those are good to know about, but far too gentle a rolloff for what I was looking for. However the LPF may be just the ticket for my next question, which I will start another thread for...

[clepsydrae]

New thread of craziness: https://forum.cockos.com/showthread....11#post2074211

[flicker177]

Here's a great website that explains a lot of this stuff in a way that even I can begin to understand it:

http://www.earlevel.com/

Check it out...

--Bill

[jcjr]

Originally Posted by clepsydrae

And thanks @jcjr for an intuitive explanation; I still only half-grasp the math going on here, but I think I understand the limitations of the technique. The FFT/iFFT method is convenient for me because I can adjust bandwidth and so forth without needing to design filter coefficients, even if I now understand there is no free lunch in terms of the ringing, etc.

I am a dunce wit ze math. Only wished to mention that it is a "fundamental problem" not generally solvable simply by finding some "ideal filter".

Would like to second the earlier mention of Goertzel filter if you wish to do "fft-like" frequency analysis but only at a single frequency of interest. Also, as best I recall you can tune your Goertzel to any arbitrary center frequency, not just the bin frequencies mathematically related to sample rate and FFT size.

Long ago when I used a Goertzel for "chromatic filter bank" analysis, I was analyzing limited-length blocks of audio, for instance maybe an eighth-note or quarter-note duration. In that usage, I windowed the audio data before the Goertzel, Blackman or some other window can't recall. I don't know if the windowing would be necessary to suppress side-bands, but if it behaves similar to FFT then I'd guess you would want to window the data in similar fashion if analyzing fixed-size blocks, not just running continuous audio thru the filter.

The nearest I recall reading for "tweaking Heisenberg's nose" involved heavily-overlapped windowed FFT. Details should be looked up before trying, exactly how to window and such. But the basic idea-- Maybe take 64 samples and put them in a zero-padded big FFT block. Maybe 4096 or even 65536 FFT block size with only 64 audio samples and the rest zeros, properly windowed.

That way you have fine-grained harmonic analysis of a very short piece of audio, pretty good time-resolution and pretty good frequency-resolution. It would just be very slow, going crab-wise thru audio only a few samples at a time, each FFT analyzing big blocks of mostly zeros. But it is still an "uncertainty trade-off".

[clepsydrae]

Thanks -- I'm not trying to detect a known frequency, I'm trying to precisely measure the pitch of a sound. Unless I misunderstand, it seems Goertzel isn't the best for that.

I've moved on to autocorrelation as a method now (hence the other thread) and it looks promising so far... just have some remaining filter-related kinks to work out.