Oversampling for an EQ

[JCJR]

Apologies if this is thread-drift but got curious enough about the "jos allpass delay" idea of oversampling to put together a simple test, to look for aliasing and freq response anomalies.

The attached jsfx code should run ok if pasted into a new reaper FX instance.

Because the upsampling is several one-pole IIR allpass in parallel, one can't simply upsample then pick off 1 out of N samples and measure for frequency response and aliasing.

For instance, if 4 X oversampling, the 0.75 sample delayed branch would have flat freq response and no aliasing. And the 0.5 sample delayed branch would have flat freq response and no aliasing. And the same for the 0.25 sample delayed branch and the non-delayed branch.

Therefore to get it to "show bad warts" the branches need to be interleaved together and fed serially thru an IIR filter network of some kind, and then pick out 1 out of 4 of the IIR filter network outputs, without any antialias filtering, to measure what this silly abuse would do to the frequency response and aliasing.

The test project code can be easily modified. As written, it has a quick'n'dirty sine sweep generator that runs from 20 Hz to 20 kHz, about 1 second per octave, and holds about 20 kHz for 1 second at the end before shutting off. So if you put the plugin on an empty reaper track, hilite about an 11 second region from the beginning and render to disk, you get a stereo test output file.

I am running Reaper at 44.1 k samplerate so it is running the test filter network at 4X oversample in my testing so far.

The left stereo output channel is the raw output from the quick'n'dirty sine sweep generator, peaking at -6 dB FS. TURN DOWN YER SPEAKERS IF YER EARS OR TWEETERS ARE FRAGILE!! YOU HAVE BEEN WARNED!!

The right stereo output channel oversamples 4X by interleaving the sine sweep samples with three JOS allpass IIR delays set for 0.75, 0.5, and 0.25 sample delays. For the "test filter network" I feed the 4X stream thru 4 instances of Vadim's first-order ZDF allpass filter, all of them tuned to 1 kHz. Then 1 out of 4 sample outputs from this four-filter cascade is written to the right stereo output, to accomplish a crude nasty downsample with no antialias filtering at all. Just discarding 3 out of 4 filter output samples.

I picked four poles of allpass filter for the "EQ filter stand-in" because it ought to have flat frequency response though of course with all that phase shifting the right channel won't null against the un-processed left channel.

Another way to test might be something like put a non-oversampled lowpass, shelf or peaking filter on the left not-oversampled output and put the same-tuned oversampled lowpass, shelf or peaking filter on the right oversampled output. However that would not null either, even if working perfectly, because of BLT frequency warping. Regardless whether the test filters would be tuned high or low, the oversampled filter network would have "less warped" response than the non-oversampled filter network.

So anyway I'm not trying to "prove a point" just honestly wondering how to realistically evaluate this silly way to oversample. I could try tuning the cascade of allpass filters to something like 100 Hz or 10000 Hz to see if changing the time constants changes the behavior.

However, with the four-pole allpass cascade tuned to 1 kHz, analyzing in Cool Edit Pro looking at the spectral window display, I'm not seeing any obvious increase of aliasing on the oversampled and crudely decimated right channel. The left unprocessed and right processed low level harmonics are about the same level at all frequencies 20 to 20 kHz. Even at 20 kHz the loudest low-level "noise" harmonics are no more than -115 dB below the frequency peak. Maybe some other tool would show different results.

However, the crude up/downsampled processed file shows a little bit of gain starting around 10 kHz and rising up to about +2.5 dB gain at 20 kHz. Which is surprising because I would have expected a little bit of FIR lowpass filtering, maybe a few dB of loss in the top octave. I had originally suspected maybe the gain is some kind of symptom of aliasing, though the spectral display doesn't show aliasing splatter. Though maybe in this test the aliasing "just happens to reflect back to the same frequency" or something? Dunno.

I can post screen shots or the test WAV if somebody wants to see them. A two-tone intermodulation test might be something useful. Maybe this silly oversampling method would exaggerate intermodulation, dunno.

Maybe there are better ways to show the warts. Maybe this particular test is somehow too unfairly sympathetic.

Code: Select all

desc:TestAllpassUpsample
// James Chandler Jr
//   http://www.errnum.com
// License: GPL - http://www.gnu.org/licenses/gpl.html

@init
  //should only call js_malloc within js @init section
  JS_MALLOC_NO_CLEARMEM = 0;
  JS_MALLOC_CLEARMEM = 1;
  g_next_js_malloc = 0;  
  function js_malloc(a_size, a_clearmem)
  local (l_result)
  (
    l_result = g_next_js_malloc;
    (a_clearmem != 0) ?
      memset(l_result, 0, a_size);
    g_next_js_malloc += a_size;
    l_result;
  );
  
  //initialize sin lookup table
  SIN_TBL_SIZE = 4096;
  SIN_TBL = js_malloc(SIN_TBL_SIZE + 1, JS_MALLOC_CLEARMEM);
  TWO_PI = 2 * $pi;
  _phaseinc = TWO_PI / 4096;
  _phase = _phaseinc;
  SIN_TBL[0] = 0;
  SIN_TBL[4096] = 0;
  //SIN_TBL[0] = 0 and also SIN_TBL[4096] = 0 to ease phase wrapping interpolation
  _i = 1;
  while (_i < 4096) //fill out elements [1] thru [4095]
  (
    SIN_TBL[_i] = sin(_phase);
    _i += 1;
    _phase += _phaseinc;
  );
  
  //function SinSweep_Init(a_BeginFreq_Hz, a_EndFreq_Hz, a_DurationSecs, a_PeakAmplitude)
  function SinSweep_Init() //20 to 20 kHz 1 octave per second
  (
    this.AmpMul = 0.5;
    this.BeginFreq_Hz = 20;
    this.BeginFreq_SinTblStepsPerSec = this.BeginFreq_Hz * SIN_TBL_SIZE;
    this.BeginFreq_SinTblIncPerSample = this.BeginFreq_SinTblStepsPerSec / srate;
    this.EndFreq_Hz = 20000;
    this.EndFreq_SinTblStepsPerSec = this.EndFreq_Hz * SIN_TBL_SIZE;
    this.EndFreq_SinTblIncPerSample = this.EndFreq_SinTblStepsPerSec / srate;
    this.OneOctavePerSecIncMul = 2 ^ (1 / srate);
    this.SinTblPhaseIdx = 0;
    this.SinTblIncPerSample = this.BeginFreq_SinTblIncPerSample;
    this.EndFreqPlaySampleCount = 0; //play the top 20 kHz for 1 second before outputting silence
  );
  
  function SinSweep_DoSamp()
  local (l_result, l_int, l_frac)
  (
    l_int = floor(this.SinTblPhaseIdx);
    l_frac = this.SinTblPhaseIdx - l_int;
    l_result = SIN_TBL[l_int] * (1 - l_frac) + SIN_TBL[l_int + 1] * l_frac; //linear interpolate
    this.SinTblPhaseIdx += this.SinTblIncPerSample;
    (this.SinTblPhaseIdx >= SIN_TBL_SIZE) ? //phase wrap
      this.SinTblPhaseIdx -= SIN_TBL_SIZE;
    (this.SinTblIncPerSample < this.EndFreq_SinTblIncPerSample) ? //gradually inc freq up to 20 kHz
      this.SinTblIncPerSample *= this.OneOctavePerSecIncMul
    :
    (
      (this.EndFreqPlaySampleCount < srate) ?
        this.EndFreqPlaySampleCount += 1
      :
        l_result = 0; //silent after 1 second of 20 kHz, don't risk frying deaf people's tweeters
    );
    l_result *= this.AmpMul;
  );  
  
  
  //https://ccrma.stanford.edu/~jos/pasp/First_Order_Allpass_Interpolation.html
  //https://www.dsprelated.com/freebooks/pasp/Delay_Line_Interpolation.html
  function JOS_AP_FracDly_Init(a_FracSampDly)
  (
    this.FracSampDly = max(min(a_FracSampDly, 0.99), 0.01);
    this.k = (1 - this.FracSampDly) / (1 + this.FracSampDly);
    this.p1 = 0.0;
    this.last_p1 = 0.0;
  );

  //only call after filter was created via _Init
  function JOS_AP_FracDly_SetFracSampDly(a_FracSampDly)
  (
    this.FracSampDly = max(min(a_FracSampDly, 0.99), 0.01);
    this.k = (1 - this.FracSampDly) / (1 + this.FracSampDly);
  );

  //only call after filter was created via _Init
  //Returns the delayed sample
  function JOS_AP_FracDly_DoSamp(a_InSamp)
  local (l_result)
  (
    this.p1 = a_InSamp - (this.k * this.last_p1);
    l_result = this.last_p1 + (this.k * this.p1);
    this.last_p1 = this.p1;
    l_result;
  );
  
  //First order trapezoidal filter object
  //Code adapted from Vadim Zavalishin's book "The Art of VA Filter Design"
  FIRST_ORD_FILTTYPE_LOPASS = 0;
  FIRST_ORD_FILTTYPE_HIPASS = 1;
  FIRST_ORD_FILTTYPE_ALLPASS_ADV = 2; //+180 degrees phase shift at DC, descending to 0 degrees phase shift at nyquist
  FIRST_ORD_FILTTYPE_ALLPASS_RET = 3; //0 degrees phase shift at DC, descending to -180 degrees phase shift at nyquist
  
  //FiltTyps: Use one of the above to set the return value of FirstOrdTrapezoidFilter_DoSamp()
  //        : However, all values are simultaneously accessible after calling DoSamp() by reading
  //        : TheFilter.lp, TheFilter.hp, TheFilter.ap_A, TheFilter.ap_R
  //a_FiltFC: Filter center frequency in Hz
  //a_SampRate: Samplerate of filter
  function FirstOrdTrapezoidFilter_Init(a_FiltType, a_FiltFC, a_SampRate)
  (
    this.FT = a_FiltType;
    this.SR = a_SampRate;
    this.Nyquist = floor(this.SR * 0.5);
    this.FC = min(a_FiltFC, this.Nyquist - 1);

    this.s = 0.0;
    this.lp = 0;
    this.hp = 0;
    this.ap_A = 0;
    this.ap_R = 0;
    
    //calculate coefficient 
    this.g = tan($pi * this.FC / this.SR);
    this.g /= (1 + this.g);
  );
  
  function FirstOrdTrapezoidFilter_SetFC(a_FiltFC)
  (
    this.FC = min(a_FiltFC, this.Nyquist - 1);
    this.g = tan($pi * this.FC / this.SR);
    this.g /= (1 + this.g);
  );
  
  //Returns the filtered sample
  function FirstOrdTrapezoidFilter_DoSamp(a_InSamp)
  local (l_v, l_result)
  (
    //Vadim Zavalishin code
    //v = (x-z1_state)*g/(1+g);
    //y = v + z1_state;
    //z1_state = y + v;
    l_v = (a_InSamp - this.s) * this.g;
    this.lp = l_v + this.s;
    this.s = this.lp + l_v;
    this.hp = a_InSamp - this.lp;
    this.ap_A = this.hp - this.lp;
    this.ap_R = this.lp - this.hp;
    
    (this.FT == FIRST_ORD_FILTTYPE_LOPASS) ?
      l_result = this.lp
    :
      (this.FT == FIRST_ORD_FILTTYPE_HIPASS) ?
        l_result = this.hp
      :
        (this.FT == FIRST_ORD_FILTTYPE_ALLPASS_ADV) ?
          l_result = this.ap_A
        :
          l_result = this.ap_R;
    l_result;
  );
   
  (srate >= 176400) ?
    g_OvrSmpMul = 1.0
  :
  (
    (srate >= 88200) ?
      g_OvrSmpMul = 2.0
    :
      g_OvrSmpMul = 4.0;
  );  

  o_SineSweeper.SinSweep_Init();
  
  //init JOS fractional delay filters
  o_JosApDly_075_L.JOS_AP_FracDly_Init(0.75);
  o_JosApDly_050_L.JOS_AP_FracDly_Init(0.50);
  o_JosApDly_025_L.JOS_AP_FracDly_Init(0.25);
  o_JosApDly_075_R.JOS_AP_FracDly_Init(0.75);
  o_JosApDly_050_R.JOS_AP_FracDly_Init(0.50);
  o_JosApDly_025_R.JOS_AP_FracDly_Init(0.25);
 
  //init test first order allpass filters
  o_TestAllpassFilt_1_R.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 1000, srate * g_OvrSmpMul);
  o_TestAllpassFilt_2_R.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 1000, srate * g_OvrSmpMul);
  o_TestAllpassFilt_3_R.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 1000, srate * g_OvrSmpMul);
  o_TestAllpassFilt_4_R.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 1000, srate * g_OvrSmpMul);


@sample
  (play_state == 1) ? //if playing
  (
    s_Tmp_Sine = o_SineSweeper.SinSweep_DoSamp();
    
    (g_OvrSmpMul >= 4) ? //4X oversample
    (
      //s_Tmp_L = o_JosApDly_075_L.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      //s_Tmp_L = o_JosApDly_050_L.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      //s_Tmp_L = o_JosApDly_025_L.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
  
      s_Tmp_R = o_JosApDly_075_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      s_Tmp_R = o_TestAllpassFilt_1_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_2_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_3_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_4_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
  
      s_Tmp_R = o_JosApDly_050_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      s_Tmp_R = o_TestAllpassFilt_1_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_2_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_3_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_4_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
  
      s_Tmp_R = o_JosApDly_025_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      s_Tmp_R = o_TestAllpassFilt_1_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_2_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_3_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      s_Tmp_R = o_TestAllpassFilt_4_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
    )
    :
    (
      (g_OvrSmpMul >= 2) ? //2X oversample
      (
        //s_Tmp_L = o_JosApDly_050_L.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
  
        s_Tmp_R = o_JosApDly_050_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
        s_Tmp_R = o_TestAllpassFilt_1_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
        s_Tmp_R = o_TestAllpassFilt_2_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
        s_Tmp_R = o_TestAllpassFilt_3_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
        s_Tmp_R = o_TestAllpassFilt_4_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
      ); 
    );
    
    //left side output is original sine sweep
    spl0 = s_Tmp_Sine; 
  
    //right side output is "crude decimated" oversampled filter by taking the output
    //  of the last pass thru the allpass filter cascade
    s_Tmp_R = o_TestAllpassFilt_1_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_Sine);
    s_Tmp_R = o_TestAllpassFilt_2_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
    s_Tmp_R = o_TestAllpassFilt_3_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
    spl1 = o_TestAllpassFilt_4_R.FirstOrdTrapezoidFilter_DoSamp(s_Tmp_R);
  )
  : //else if not playing, output silence
  (
    spl0 = 0;
    spl1 = 0;
  );

[matt42]

Hi JCJR,

I have time on my hands today, so I wrote 3 versions of the same EQ plugin for testing.

Version 1: Standard BLT peak filter.
Version 2: BLT peak filter with x4 upsampling using a 20th order elliptical filter
Version 3: BLT peak filter with x4 upsampling using allpass interpolation

All the versions use the same parameter set for the peak filter:
Gain = +10dB. Frequency = 15 kHz. Q = 0.1. Base sample rate 44100.

Here are the frequency response curves:

Standard BLT peak filter:


Obvious cramping at nyquist but also left side of slope is narrower than the analog prototype.

Elliptical x4:


Better matched response. Can see the sharp elliptical filter cut off at nyquist.

Allpass Interpolation x4:


Phase shifted aliases are messing up the high frequency response.

[matt42]

JCJR,

I think this is also relevant to the discussion and might also explain your test results - Below is a plot of the allpass upsampled EQ's harmonic response to a 10 kHz sine:


No distortion! The out of phase aliases that mess up the frequency response are mirrors of the original signal, they reflect back down perfectly to the original partials and so don't cause distortion. They will however mess with the frequency response depending on their phase and will cause issues if we perform non-linear processing at the higher sample rate. I guess the implication is we can use the lowest order upsampling filter that produces acceptable ripple in the high end response as we don't have to worry about distortion.

EDIT: Sorry last graph. I also graphed the full band response of the all pass interpolator (I'm very bored with nothing to do today ). The original signal was at 11025, so nyquist was 5512.5 kHz then upsampled to 44100 with the allpass interpolation:


So, there's filtering in the pass band and not much rejection of the mirrors above nyquist (these will fold down and shift the passband response after downsampling).

(Also for full disclosure the test is a little fast and dirty. I decimated by a factor of 4 using the elliptical filter and reconstructed the signal with the interpolator. There's a tight notch at 5512.5 from the elliptical filter decimation which is mirrored at 16537.5. Otherwise the response is accurate.)

[Max M.]

Yes, the "no distortion" thing is expected since just upsampling/downsampling (with or without any additional linear filtering in between) does not introduce any "new" harmonics (unlike non-linear things similar to for example a waveshaping would do). The mirrored harmonics simply return to their original places in the end.

... we can use the lowest order upsampling filter that produces acceptable ripple

And that would be a problem because this will require both amplitude and phase in the given region to be matched for "real" and "mirror" parts in a specific way. Now count that there they (amplitude and phase) also depend on your main EQ filter, thus, in a best case: any change in the main EQ will also require a corresponding change in the allpass, and in the worst case: no allpass setting will be able to produce the required phase shift (unless it's a high-order allpass able of an arbitrary phase curve).

I'm not saying the idea can't work at all, it's just about complexity it brings in. E.g. try to test even more minimalist version of the scheme:

• - x2 upsampling (zero-padding)
  - first order allpass
  - first order BLT lowpass
  - x2 decimation

By tweaking the allpass frequency we will be able to get all sorts of filters in a range [from BLT lowpass to subtly cutting highshelf] with "matching lowpass" curves also present among them.
With two issues:
a. changes in the main EQ (the BLT lowpass in this case) require changes in the allpass (unlike standard antialias filtering where the filter is constant).
b. the resulting phase curve is similar to what you would get with a regular "matching designs".

I guess we may fight "b" by increasing the oversampling factor and/or increasing the allpass order. But "a" is pretty discouraging since for the second-order EQ and higher oversampling ratios we'll need quite transcendental trigonometry to compute the required allpass phase-shift.

At first glance I suspect that even a second-order "quasi-linear phase" antialiasing filter (with comparable amount of computations per sample and not depending on the main EQ) would be a better alternative (but it's difficult to say for sure w/o further analysis).

[JCJR]

Thanks very much Matt for showing how to demonstrate a problem. Yes I was just wondering about the allpass interpolation only for linear processing such as nX oversampled EQ. I wouldn't have expected it to be of much use, at least without additional antialias filtering, for resampling or nonlinear processing.

Interesting that the aliases reflect down to the original frequencies and don't show distortion harmonics (as was "puzzling" in my test). I'd guessed maybe that was a possibility. Would be interesting to see if that would still hold true on a 2-tone intermodulation test.

Maybe I'll port some EQ filters into jsfx and play with it some more. Or not. Long ago I was real impressed with the performance of those polyphase halfband filters considering the "relatively low number" of math ops necessary, but OTOH the op count is fairly high compared to stringing a few 2nd order IIRs together.

And also many thanks to Max M for additional excellent explanations.

[matt42]

Would be interesting to see if that would still hold true on a 2-tone intermodulation test.

Yes, it would hold. There can be no intermodulations with zero padding (even the allpass interpolator can be viewed in these terms) . We can intuit why. Firstly we can add signals together and be sure that the result only contains what we put in, no new intermodulations can occur. So if we add some sines together we just have the harmonic content we put in. If we zero pad, then the non zero samples of this signal exactly equal the sum of the original sines. We'd get the same result if we zero padded separately and then sum, or if we combined first then zero pad. So logically no intermodulation is occurring.

[matt42]

Hi Max M,

I might have misunderstood your reply. Anyway regarding the "lowest order... acceptable" comment - I just meant, literally, whatever is good enough. Perhaps with a not quite so steep slope we get a bit of ripple near nyquist at least we don't have to be worried about distortion. For me this would still be a reasonably high order, and I'm definitely not suggesting the allpass scheme. But perhaps your comment is more general and I've misunderstood.

[JustinJ]

After a bit more work on the FIR filtering and a lot more digging about, a few revelations.

My problem isn't that there's some high end leakage. I've been using Voxengo's SPAN which has a default tilt in the spectrum which more readily showed high end. Getting rid of the tilt shows what I think is a noise reduction issue. So, I'm thinking my filter isn't reducing the audio enough. It's getting to around -60db which is still audible. I think it's the aliasing energy bouncing back.

What I need is a FIR filter that will bounce the aliasing back at least -110db which won't be audible. According to another forum post by Andy Simper it's possible to get this with 143 taps with a 71 sample latency.

[Max M.]

If I'm not mistaken that -60dB supression of the "mirror" would mean about 0.0087dB error (-60db -> 0.001; (1 - 0.001) -> -0.0087dB) in your final EQ curve (roughly, assuming it's flying somewhere around unity). And I suspect it is (in general) a lower error than the error you get by having only x2 oversampling. I.e. dont't forget that while aliasing is the only source of the error in your null test, it will be not the only source of your EQ vs. "ideal EQ" curve mismatch.

[JustinJ]

Mission accomplished!

I used the following FIR designer: http://arc.id.au/FilterDesign.html

To produce coefficients for a 65 tap Kaiser FIR filter using the following definitions:

Fa = 0
Fb = 22050
length = 65
Fs = 88200
Att = 96db

It's a steep low pass that does an extremely good job of getting rid of aliasing beyond Nyquist/2 (44.1kHz in my case). It filters any reflective signal into the -96db and lower range making it inaudible.

It has a latency of 16 samples.

Here's the unrolled optimized JSFX code:

Code: Select all

function init_FIR_filter()
  instance(
    b0, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11, b12, b13, b14, b15, b16, b17, b18, b19,
    b20, b21, b22, b23, b24, b25, b26, b27, b28, b29, b30, b31, b32, b33, b34, b35, b36, b37, b38, b39,
    b40, b41, b42, b43, b44, b45, b46, b47, b48, b49, b50, b51, b52, b53, b54, b55, b56, b57, b58, b59,
    b60, b61, b62, b63, b64, b65,
    w0, w1, w2, w3, w4, w5, w6, w7, w8, w9, w10, w11, w12, w13, w14, w15, w16, w17, w18, w19,
    w20, w21, w22, w23, w24, w25, w26, w27, w28, w29, w30, w31, w32, w33, w34, w35, w36, w37, w38, w39,
    w40, w41, w42, w43, w44, w45, w46, w47, w48, w49, w50, w51, w52, w53, w54, w55, w56, w57, w58, w59,
    w60, w61, w62, w63, w64, w65
  )
(
  // Set the coefficients taking into account they're mirrored around the centre tap.
  w0 = w64 = -0.000000;
  w1 = w63 = -0.000016; 
  w2 = w62 = 0.000000;
  w3 = w61 = 0.000067; 
  w4 = w60 = -0.000000; 
  w5 = w59 = -0.000189; 
  w6 = w58 = -0.000000; 
  w7 = w57 = 0.000438;
  w8 = w56 = -0.000000; 
  w9 = w55 = -0.000891; 
  w10 = w54 = 0.000000; 
  w11 = w53 = 0.001653; 
  w12 = w52 = -0.000000; 
  w13 = w51 = -0.002862; 
  w14 = w50 = 0.000000; 
  w15 = w49 = 0.004690; 
  w16 = w48 = -0.000000; 
  w17 = w47 = -0.007364; 
  w18 = w46 = 0.000000; 
  w19 = w45 = 0.011189; 
  w20 = w44 = -0.000000; 
  w21 = w43 = -0.016631; 
  w22 = w42 = 0.000000; 
  w23 = w41 = 0.024500; 
  w24 = w40 = -0.000000; 
  w25 = w39 = -0.036479; 
  w26 = w38 = 0.000000; 
  w27 = w37 = 0.056928; 
  w28 = w36 = -0.000000; 
  w29 = w35 = -0.101933; 
  w30 = w34 = 0.000000; 
  w31 = w33 = 0.316897; 
  w32 = 0.500000; 
);

function do_FIR_filter(input)
  instance(
    b0, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11, b12, b13, b14, b15, b16, b17, b18, b19,
    b20, b21, b22, b23, b24, b25, b26, b27, b28, b29, b30, b31, b32, b33, b34, b35, b36, b37, b38, b39,
    b40, b41, b42, b43, b44, b45, b46, b47, b48, b49, b50, b51, b52, b53, b54, b55, b56, b57, b58, b59,
    b60, b61, b62, b63, b64, b65,
    w0, w1, w2, w3, w4, w5, w6, w7, w8, w9, w10, w11, w12, w13, w14, w15, w16, w17, w18, w19,
    w20, w21, w22, w23, w24, w25, w26, w27, w28, w29, w30, w31, w32, w33, w34, w35, w36, w37, w38, w39,
    w40, w41, w42, w43, w44, w45, w46, w47, w48, w49, w50, w51, w52, w53, w54, w55, w56, w57, w58, w59,
    w60, w61, w62, w63, w64, w65
  )  
(
  // Scroll the window
  b65 = b64; b64 = b63; b63 = b62; b62 = b61; b61 = b60;
  b60 = b59; b59 = b58; b58 = b57; b57 = b56; b56 = b55; b55 = b54; b54 = b53; b53 = b52; b52 = b51; b51 = b50;
  b50 = b49; b49 = b48; b48 = b47; b47 = b46; b46 = b45; b45 = b44; b44 = b43; b43 = b42; b42 = b41; b41 = b40;
  b40 = b39; b39 = b38; b38 = b37; b37 = b36; b36 = b35; b35 = b34; b34 = b33; b33 = b32; b32 = b31; b31 = b30;
  b30 = b29; b29 = b28; b28 = b27; b27 = b26; b26 = b25; b25 = b24; b24 = b23; b23 = b22; b22 = b21; b21 = b20;
  b20 = b19; b19 = b18; b18 = b17; b17 = b16; b16 = b15; b15 = b14; b14 = b13; b13 = b12; b12 = b11; b11 = b10;
  b10 = b9; b9 = b8; b8 = b7; b7 = b6; b6 = b5; b5 = b4; b4 = b3; b3 = b2; b2 = b1; b1 = b0;

  // Introduce the new sample into the window
  b0 = input;

  // Compute and return
  b0*w0 + b1 * w1 + b2 * w2 + b3 * w3 + b4 * w4 + b5 * w5 + b6 * w6 + b7 * w7 + b8 * w8 + b9 * w9 +
  b10*w10 + b11 * w11 + b12 * w12 + b13 * w13 + b14 * w14 + b15 * w15 + b16 * w16 + b17 * w17 + b18 * w18 + b19 * w19 +
  b20*w20 + b21 * w21 + b22 * w22 + b23 * w23 + b24 * w24 + b25 * w25 + b26 * w26 + b27 * w27 + b28 * w28 + b29 * w29 +
  b30*w30 + b31 * w31 + b32 * w32 + b33 * w33 + b34 * w34 + b35 * w35 + b36 * w36 + b37 * w37 + b38 * w38 + b39 * w39 +
  b40*w40 + b41 * w41 + b42 * w42 + b43 * w43 + b44 * w44 + b45 * w45 + b46 * w46 + b47 * w47 + b48 * w48 + b49 * w49 +
  b50*w50 + b51 * w51 + b52 * w52 + b53 * w53 + b54 * w54 + b55 * w55 + b56 * w56 + b57 * w57 + b58 * w58 + b59 * w59 +
  b60*w60 + b61 * w61 + b62 * w62 + b63 * w63 + b64 * w64 + b65 * w65;
);

Thanks for everyone for helping on here and the fascinating conversation.

[JCJR]

Mission accomplished!

Hi Justin, I'm glad you got something you like out of it.

[JCJR]

Max I am ignorant but suspect that only the most golden of pristine ears could hear naked oversampling high frequency phase oddities unless such oddities get massively big. My possibly wrongheaded concern that BLT phase should be warped exactly the same amount as the BLT amplitude was already mentioned by Vadim-- Hoping that multiple bands might mix together "similar to analog expectations" even if the actual center frequencies might not agree with analog in all cases. If an entire filter network wrapped inside an up/downsample block "works about as expected" then it is an open question whether some minor HF phase anomalies from the up/downsample process could be heard by any except the most golden pristine ears. OTOH maybe minor up/downsample HF phase oddities would be easily audible. Dunno.

For high quality resampling or oversampled mastering processing or other critical tasks, it would seem worthy to use only linear phase and high fidelity up/downsampling. Do as little accidental damage as is possible just in case somebody somewhere might be able to hear it. But IMO a causal EQ screws up the phase so drastically that slight up/downsampling phase oddities might not make much difference? (At least so long as the up/downsampling phase weirdness doesn't mess up the enclosed EQ filter network's amplitude response.)

So anyway I'm retired and don't know why I recently got re-interested in allpass resampling. Was always interested in low-parts-count electronic circuits. Maybe it is about the same thing "low parts count DSP". Interesting to consider even if it ultimately turns out "not very useful". There is a certain charm in the one or two-transistor circuit or whatever.

I don't dispute your idea that allpass upsampling for EQ is probably a dead end useless but remain curious about it. Did some further testing following Matt42's example. The possibility that MAYBE some very lightweight antialias post-filtering applied to an allpass upsample could clean the signal enough to un-do most of the aliased frequency response interaction with the EQ filters.

The first step was to closer-test the behavior without any antialias filtering, as described in the next message. A possible further step might be to test whether a cheap anti-alias filter could "mostly fix the major problems".

[JCJR]

Hi Matt42. I imported my old RBJ filters object into jsfx and did some testing on the allpass oversampling. With the "no antialias filtering" I saw problems but not as severe as you found. Three possible guesses about the discrepancies--

_1 I used bigger peaking filter Q than you did.
_2 Maybe there are procedural differences in what got tested (jsfx code attached at end of message).
_3 Maybe the code is so buggy that it just accidentally shows fewer problems but otherwise a meaningless test.

Am inclined to guess _1 most probable. You stated a peaking filter Q of 0.2. Long ago when I originally wrote the RBJ filter object, it was set to clamp Q in the range of 0.5 to 100. I realize that Q can sometimes be lower than 0.5 but so far as I knew 0.5 was the lowest Q in "conventional designed 2nd order active filters".

So I used Q = 0.5, gain = +12 dB. As best I recall RBJ accomplishes peaking filter constant bandwidth by increasing nominal Q along with increasing boost and decreasing Q along with increasing cut. So a nominal peaking filter Q of 0.5 would be a bigger "actual Q" with 12 dB boost, and would use less than 0.5 for any amounts of filter cut.

I modified the sine sweep to hold 20 Hz for 1 second at the beginning, then sweep at 1 octave per second, then hold for 1 second at 20 kHz at the end. Which was intended to maybe make it easier on FFT based spectral analysis on the file heads and tails. The following screen shots, Voxengo Span still shows a glitchy artifact around 20 Hz that isn't in the actual audio. The Cooledit Pro waveform screenshots are zoomed-in to ignore the extra 1 second steady freq hold on the heads and tails.

The filter input sine sweep peaks at -15 dB FS, so the max peaking filter output should kiss -3 dB in the CoolEdit Pro screenshots. As best I recall RBJ defined bandwidth at the middle-gain positions (in log space) so the bandwidth of each filter test is the distance between -9 dB intersections. Because the sweep is 1 octave per second, it can be "kinda sorta accurately estimated" in CoolEdit Pro by hiliting that region and reading off the region size in seconds. Ferinstance if it happens to be 2.6 seconds then it is 2.6 octaves wide. But this is not a super-accurate way to do it. Just kinda-sorta a purt close estimate.

Apologies there are so many screenshots. It was just an interesting diversion for the day. Knowing what warts to look for, I might see how effective some lighhtweight IIR antialias filtering might "mostly fix up" the artifacts. Perhaps it will turn out impossible to do it in a lightweight fashion.









You can change the RBJ test filter frequency on lines 393 and 394:
o_PeakingFilt_L.TRBJFilt_Init(RBJ_FILTTYPE_PEAKING, 15000, 0.5, 4.0, srate);
o_PeakingFilt_R.TRBJFilt_Init(RBJ_FILTTYPE_PEAKING, 15000, 0.5, 4.0, srate * g_OvrSmpMul);

Code: Select all

desc:TestAllpassUpsample
// James Chandler Jr
//   http://www.errnum.com
// License: GPL - http://www.gnu.org/licenses/gpl.html

@init
  //should only call js_malloc within js @init section
  JS_MALLOC_NO_CLEARMEM = 0;
  JS_MALLOC_CLEARMEM = 1;
  g_next_js_malloc = 0;  
  function js_malloc(a_size, a_clearmem)
  local (l_result)
  (
    l_result = g_next_js_malloc;
    (a_clearmem != 0) ?
      memset(l_result, 0, a_size);
    g_next_js_malloc += a_size;
    l_result;
  );
  
  //initialize sin lookup table
  SIN_TBL_SIZE = 4096;
  SIN_TBL = js_malloc(SIN_TBL_SIZE + 1, JS_MALLOC_CLEARMEM);
  TWO_PI = 2 * $pi;
  _phaseinc = TWO_PI / 4096;
  _phase = _phaseinc;
  SIN_TBL[0] = 0;
  SIN_TBL[4096] = 0;
  //SIN_TBL[0] = 0 and also SIN_TBL[4096] = 0 to ease phase wrapping interpolation
  _i = 1;
  while (_i < 4096) //fill out elements [1] thru [4095]
  (
    SIN_TBL[_i] = sin(_phase);
    _i += 1;
    _phase += _phaseinc;
  );
  
  function SinSweep_Init() //20 to 20 kHz 1 octave per second
  (
    this.AmpMul = 0.1778; //peaking at -15 dB FS
    this.BeginFreq_Hz = 20;
    this.BeginFreq_SinTblStepsPerSec = this.BeginFreq_Hz * SIN_TBL_SIZE;
    this.BeginFreq_SinTblIncPerSample = this.BeginFreq_SinTblStepsPerSec / srate;
    this.EndFreq_Hz = 20000;
    this.EndFreq_SinTblStepsPerSec = this.EndFreq_Hz * SIN_TBL_SIZE;
    this.EndFreq_SinTblIncPerSample = this.EndFreq_SinTblStepsPerSec / srate;
    this.OneOctavePerSecIncMul = 2 ^ (1 / srate);
    this.SinTblPhaseIdx = 0;
    this.SinTblIncPerSample = this.BeginFreq_SinTblIncPerSample;
    this.BeginFreqPlaySampleCount = 0; //play the bottom 20 Hz for 1 secong before beginning sweep
    this.EndFreqPlaySampleCount = 0; //play the top 20 kHz for 1 second before outputting silence
  );
  
  function SinSweep_DoSamp()
  local (l_result, l_int, l_frac)
  (
    l_int = floor(this.SinTblPhaseIdx);
    l_frac = this.SinTblPhaseIdx - l_int;
    l_result = SIN_TBL[l_int] * (1 - l_frac) + SIN_TBL[l_int + 1] * l_frac; //linear interpolate
    this.SinTblPhaseIdx += this.SinTblIncPerSample;
    (this.SinTblPhaseIdx >= SIN_TBL_SIZE) ? //phase wrap
      this.SinTblPhaseIdx -= SIN_TBL_SIZE;
    (this.SinTblIncPerSample < this.EndFreq_SinTblIncPerSample) ? //gradually inc freq up to 20 kHz
    (
      (this.BeginFreqPlaySampleCount < srate) ?
        this.BeginFreqPlaySampleCount += 1
      :
        this.SinTblIncPerSample *= this.OneOctavePerSecIncMul;
    )
    :
    (
      (this.EndFreqPlaySampleCount < srate) ?
        this.EndFreqPlaySampleCount += 1
      :
        l_result = 0; //silent after 1 second of 20 kHz, don't risk frying deaf people's tweeters
    );
    l_result *= this.AmpMul;
  );  
  
  
  //https://ccrma.stanford.edu/~jos/pasp/First_Order_Allpass_Interpolation.html
  //https://www.dsprelated.com/freebooks/pasp/Delay_Line_Interpolation.html
  function JOS_AP_FracDly_Init(a_FracSampDly)
  (
    this.FracSampDly = max(min(a_FracSampDly, 0.99), 0.01);
    this.k = (1 - this.FracSampDly) / (1 + this.FracSampDly);
    this.p1 = 0.0;
    this.last_p1 = 0.0;
  );

  //only call after filter was created via _Init
  function JOS_AP_FracDly_SetFracSampDly(a_FracSampDly)
  (
    this.FracSampDly = max(min(a_FracSampDly, 0.99), 0.01);
    this.k = (1 - this.FracSampDly) / (1 + this.FracSampDly);
  );

  //only call after filter was created via _Init
  //Returns the delayed sample
  function JOS_AP_FracDly_DoSamp(a_InSamp)
  local (l_result)
  (
    this.p1 = a_InSamp - (this.k * this.last_p1);
    l_result = this.last_p1 + (this.k * this.p1);
    this.last_p1 = this.p1;
    l_result;
  );
  
  //First order trapezoidal filter object
  //Code adapted from Vadim Zavalishin's book "The Art of VA Filter Design"
  FIRST_ORD_FILTTYPE_LOPASS = 0;
  FIRST_ORD_FILTTYPE_HIPASS = 1;
  FIRST_ORD_FILTTYPE_ALLPASS_ADV = 2; //+180 degrees phase shift at DC, descending to 0 degrees phase shift at nyquist
  FIRST_ORD_FILTTYPE_ALLPASS_RET = 3; //0 degrees phase shift at DC, descending to -180 degrees phase shift at nyquist
  
  //FiltTypes: Use one of the above to set the return value of FirstOrdTrapezoidFilter_DoSamp()
  //        : However, all values are simultaneously accessible after calling DoSamp() by reading
  //        : TheFilter.lp, TheFilter.hp, TheFilter.ap_A, TheFilter.ap_R
  //a_FiltFC: Filter center frequency in Hz
  //a_SampRate: Samplerate of filter
  function FirstOrdTrapezoidFilter_Init(a_FiltType, a_FiltFC, a_SampRate)
  (
    this.FT = a_FiltType;
    this.SR = a_SampRate;
    this.Nyquist = floor(this.SR * 0.5);
    this.FC = min(a_FiltFC, this.Nyquist - 1);

    this.s = 0.0;
    this.lp = 0;
    this.hp = 0;
    this.ap_A = 0;
    this.ap_R = 0;
    
    //calculate coefficient 
    this.g = tan($pi * this.FC / this.SR);
    this.g /= (1 + this.g);
  );
  
  function FirstOrdTrapezoidFilter_SetFC(a_FiltFC)
  (
    this.FC = min(a_FiltFC, this.Nyquist - 1);
    this.g = tan($pi * this.FC / this.SR);
    this.g /= (1 + this.g);
  );
  
  //Returns the filtered sample
  function FirstOrdTrapezoidFilter_DoSamp(a_InSamp)
  local (l_v, l_result)
  (
    //Vadim Zavalishin code
    //v = (x-z1_state)*g/(1+g);
    //y = v + z1_state;
    //z1_state = y + v;
    l_v = (a_InSamp - this.s) * this.g;
    this.lp = l_v + this.s;
    this.s = this.lp + l_v;
    this.hp = a_InSamp - this.lp;
    this.ap_A = this.hp - this.lp;
    this.ap_R = this.lp - this.hp;
    
    (this.FT == FIRST_ORD_FILTTYPE_LOPASS) ?
      l_result = this.lp
    :
      (this.FT == FIRST_ORD_FILTTYPE_HIPASS) ?
        l_result = this.hp
      :
        (this.FT == FIRST_ORD_FILTTYPE_ALLPASS_ADV) ?
          l_result = this.ap_A
        :
          l_result = this.ap_R;
    l_result;
  );
   
  //RBJ Cookbook Filter Object
  RBJ_FILTTYPE_LOPASS = 1;
  RBJ_FILTTYPE_BANDPASS = 2;
  RBJ_FILTTYPE_HIPASS = 3;
  RBJ_FILTTYPE_NOTCH = 4;
  RBJ_FILTTYPE_ALLPASS = 5;
  RBJ_FILTTYPE_PEAKING = 6;
  RBJ_FILTTYPE_LOSHELF = 7;
  RBJ_FILTTYPE_HISHELF = 8;

  //A     = sqrt[ 10^(dBgain/20) ]   = 10^(dBgain/40)    (for for peaking and shelving EQ filters only)
  //omega = 2*PI*frequency/sample_rate
  //sn    = sin(omega)
  //cs    = cos(omega)
  //alpha = sn/(2*Q)
  //      = sn*sinh[ ln(2)/2 * bandwidth * omega/sn ]     (if bandwidth is specified instead of Q)
  //beta  = sqrt[ (A^2 + 1)/S - (A-1)^2 ]   (for shelving EQ filters only)
  //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]

  function TRBJFilt_CalcInternalVars()
  local (l_A, l_omega, l_sn, l_cs, l_alpha, l_beta)
  (
    (this.FT > RBJ_FILTTYPE_HISHELF) ? //trap out of range filter type error
      this.FT = RBJ_FILTTYPE_LOPASS;
    (this.FT >= RBJ_FILTTYPE_LOSHELF) ?
      this.Q = 0.5 ^ 0.5; //set to 0.7071 for shelving filters
    l_A = pow(this.GainMul, 0.5); //for peaking and shelving EQ filters only
    l_omega = (2.0 * $pi * this.FC) / this.SR;
    l_sn = sin(l_omega);
    l_cs = cos(l_omega);
    l_alpha = l_sn / (2.0 * this.Q);
    l_beta = pow(((l_A * l_A) + 1) - ((l_A - 1) * (l_A - 1)), 0.5); //presume S var is steep as possible, 1.0

    (this.FT == RBJ_FILTTYPE_LOPASS) ?
    (
      this.B0 = (1.0 - l_cs) / 2.0;
      this.B1 = 1.0 - l_cs;
      this.B2 = this.B0;
      this.A0 = 1.0 + l_alpha;
      this.A1 = -2.0 * l_cs;
      this.A2 = 1.0 - l_alpha;
    );
    (this.FT == RBJ_FILTTYPE_BANDPASS) ?
    (
      this.B0 = l_alpha;
      this.B1 = 0.0;
      this.B2 = -l_alpha;
      this.A0 = 1.0 + l_alpha;
      this.A1 = -2.0 * l_cs;
      this.A2 = 1.0 - l_alpha;
    );
    (this.FT == RBJ_FILTTYPE_HIPASS) ?
    (
      this.B0 = (1.0 + l_cs) / 2.0;
      this.B1 = -(1.0 + l_cs);
      this.B2 = this.B0;
      this.A0 = 1.0 + l_alpha;
      this.A1 = -2.0 * l_cs;
      this.A2 = 1.0 - l_alpha;
    );
    (this.FT == RBJ_FILTTYPE_NOTCH) ?
    (
      this.B0 = 1.0;
      this.B1 = -2.0 * l_cs;
      this.B2 = 1.0;
      this.A0 = 1.0 + l_alpha;
      this.A1 = this.B1;
      this.A2 = 1.0 - l_alpha;
    );
    (this.FT == RBJ_FILTTYPE_ALLPASS) ?
    (
      this.B0 = 1.0 - l_alpha;
      this.B1 = -2.0 * l_cs;
      this.B2 = 1.0 + l_alpha;
      this.A0 = this.B2;
      this.A1 = this.B1;
      this.A2 = this.B0;
    );
    (this.FT == RBJ_FILTTYPE_PEAKING) ?
    (
      this.B0 = 1.0 + l_alpha * l_A;
      this.B1 = -2.0 * l_cs;
      this.B2 = 1.0 - l_alpha * l_A;
      this.A0 = 1.0 + l_alpha / l_A;
      this.A1 = this.B1;
      this.A2 = 1.0 - l_alpha / l_A;
    );
    (this.FT == RBJ_FILTTYPE_LOSHELF) ?
    (
      this.B0 = l_A * ((l_A + 1.0) - ((l_A - 1.0) * l_cs) + (l_beta * l_sn));
      this.B1 = 2.0 * l_A * ((l_A - 1.0) - (l_A + 1.0) * l_cs);
      this.B2 = l_A * ((l_A + 1.0) - ((l_A - 1.0) * l_cs) - (l_beta * l_sn));
      this.A0 = (l_A + 1.0) + ((l_A - 1.0) * l_cs) + (l_beta * l_sn);
      this.A1 = -2.0 * ((l_A - 1.0) + ((l_A + 1.0) * l_cs));
      this.A2 = (l_A + 1.0) + ((l_A - 1.0) * l_cs) - (l_beta * l_sn);
    );
    (this.FT == RBJ_FILTTYPE_HISHELF) ?
    (
      this.B0 = l_A * ((l_A + 1.0) + ((l_A - 1.0) * l_cs) + (l_beta * l_sn));
      this.B1 = -2.0 * l_A * ((l_A - 1.0) + (l_A + 1.0) * l_cs);
      this.B2 = l_A * ((l_A + 1.0) + ((l_A - 1.0) * l_cs) - (l_beta * l_sn));
      this.A0 = (l_A + 1.0) - ((l_A - 1.0) * l_cs) + (l_beta * l_sn);
      this.A1 = 2.0 * ((l_A - 1.0) - ((l_A + 1.0) * l_cs));
      this.A2 = (l_A + 1.0) - ((l_A - 1.0) * l_cs) - (l_beta * l_sn);
    );
    this.B0_Div_A0 = this.B0 / this.A0;
    this.B1_Div_A0 = this.B1 / this.A0;
    this.B2_Div_A0 = this.B2 / this.A0;
    this.A1_Div_A0 = this.A1 / this.A0;
    this.A2_Div_A0 = this.A2 / this.A0;
  );
  
  function TRBJFilt_ResetStateVars()
  (
    this.X_1 = 0.0;
    this.X_2 = 0.0;
    this.Y_1 = 0.0;
    this.Y_2 = 0.0;
  );
  
  function TRBJFilt_Init(a_FiltType, a_FiltFC, a_FiltQ, a_FiltGain, a_SampRate)
  //pass this.Q = 0.70710678 or pow(0.5, 0.5) for shelving filters, auto done in TRBJFilt_CalcInternalVars()
  (
    this.SR = a_SampRate;
    this.FC = min(a_FiltFC, this.SR * 0.47619); //max of 21 KHz at 44.1 K samplerate
    this.Q = min(max(a_FiltQ, 0.5), 100);
    this.GainMul = a_FiltGain;
    this.FT = a_FiltType;
    this.TRBJFilt_CalcInternalVars();
    this.TRBJFilt_ResetStateVars();
  );

  function TRBJFilt_DoSamp(a_InSamp)
  local (l_result)
  (
    l_result = (this.B0_Div_A0 * a_InSamp) + (this.B1_Div_A0 * this.X_1) + (this.B2_Div_A0 * this.X_2) -
    (this.A1_Div_A0 * this.Y_1) - (this.A2_Div_A0 * this.Y_2);
    this.Y_2 = this.Y_1;
    this.Y_1 = l_result;
    this.X_2  = this.X_1;
    this.X_1  = a_InSamp;
    (this.FT < RBJ_FILTTYPE_PEAKING) ?
      l_result *= this.GainMul;
    l_result;
  );

  //In the following setter functions, set a_ReCalc > zero on the last setter function called
  //For example, if you only change one setting at a time, maybe you only repeatedly
  //  change the filter Fc as a tracking filter or to follow user's mouse EQ adjustment, then every time
  //  you change the Fc you should pass a_ReCalc = 1 or whatever to make the filter always uses current settings.
  //However if you need to change two or more settings for tracking automation or whatever, then it
  //  would be cpu wasteful to recalc the filter everytime you only change one setting.
  //For example if you are simultaneously changing BOTH the Fc AND the Q, then the filter only needs
  //  recalc after you have set both the Fc and Q to new values.
  //This saves redundantly calling redundant wasteful recalcs.
  //For instance if you happen to be automating Fc, Q, and Gain all as a group,
  //  then you only have to recalc on the LAST SETTING you change in that group:
  //  MyFilt.TRBJFilt_SetFc(NewFc, 0); //change Fc but don't recalc yet
  //  MyFilt.TRBJFilt_SetQ(NewQ, 0); //change the Q but don't recalc yet
  //  MyFilt.TRBJFilt_SetOutGain(NewGain, 1);
  //  //finally after making the last change in your group of changes, force that one final recalc
  //  //to update ALL THREE recently changed settings with just one pass thru TRBJFilt_CalcInternalVars()
  function TRBJFilt_SetSampleRate(a_SampRate, a_ReCalc)
  (
    this.SR = a_SampRate;
    this.FC = min(this.FC, this.SR * 0.47619); //max of 21 KHz at 44.1 K samplerate
    (a_ReCalc > 0) ?
      this.TRBJFilt_CalcInternalVars();
  );

  function TRBJFilt_SetFc(a_FiltFC, a_ReCalc)
  (
    this.FC = min(a_FiltFC, this.SR * 0.47619); //max of 21 KHz at 44.1 K samplerate
    (a_ReCalc > 0) ?
      this.TRBJFilt_CalcInternalVars();
  );

  function TRBJFilt_SetQ(a_FiltQ, a_ReCalc)
  (
    this.Q = min(max(a_FiltQ, 0.5), 100);
    (a_ReCalc > 0) ?
      this.TRBJFilt_CalcInternalVars();
  );

  function TRBJFilt_SetOutGain(a_FiltGain, a_ReCalc)
  (
    this.GainMul = a_FiltGain;
    (a_ReCalc > 0) ?
      this.TRBJFilt_CalcInternalVars();
  );

  function TRBJFilt_SetFilterType(a_FiltType, a_ReCalc)
  (
    this.FT = a_FiltType;
    (a_ReCalc > 0) ?
      this.TRBJFilt_CalcInternalVars();
  );

  (srate >= 176400) ?
    g_OvrSmpMul = 1.0
  :
  (
    (srate >= 88200) ?
      g_OvrSmpMul = 2.0
    :
      g_OvrSmpMul = 4.0;
  );  

  o_SineSweeper.SinSweep_Init();
  
  //init JOS fractional delay filters
  //o_JosApDly_075_L.JOS_AP_FracDly_Init(0.75);
  //o_JosApDly_050_L.JOS_AP_FracDly_Init(0.50);
  //o_JosApDly_025_L.JOS_AP_FracDly_Init(0.25);
  o_JosApDly_075_R.JOS_AP_FracDly_Init(0.75);
  o_JosApDly_050_R.JOS_AP_FracDly_Init(0.50);
  o_JosApDly_025_R.JOS_AP_FracDly_Init(0.25);
 
  //init test peaking filters
  o_PeakingFilt_L.TRBJFilt_Init(RBJ_FILTTYPE_PEAKING, 15000, 0.5, 4.0, srate);
  o_PeakingFilt_R.TRBJFilt_Init(RBJ_FILTTYPE_PEAKING, 15000, 0.5, 4.0, srate * g_OvrSmpMul);


@sample
  (play_state == 1) ? //if playing
  (
    s_Tmp_Sine = o_SineSweeper.SinSweep_DoSamp();
    
    (g_OvrSmpMul >= 4) ? //4X oversample
    (  
      s_Tmp_R = o_JosApDly_075_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      s_Tmp_R = o_PeakingFilt_R.TRBJFilt_DoSamp(s_Tmp_R);
    
      s_Tmp_R = o_JosApDly_050_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      s_Tmp_R = o_PeakingFilt_R.TRBJFilt_DoSamp(s_Tmp_R);
    
      s_Tmp_R = o_JosApDly_025_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
      s_Tmp_R = o_PeakingFilt_R.TRBJFilt_DoSamp(s_Tmp_R);
    )
    :
    (
      (g_OvrSmpMul >= 2) ? //2X oversample
      (
        s_Tmp_R = o_JosApDly_050_R.JOS_AP_FracDly_DoSamp(s_Tmp_Sine);
        s_Tmp_R = o_PeakingFilt_R.TRBJFilt_DoSamp(s_Tmp_R);
      ); 
    );
    
    //left side output is original srate
    spl0 = o_PeakingFilt_L.TRBJFilt_DoSamp(s_Tmp_Sine);
  
    //right side output is "crude decimated" oversampled filter by taking the output
    //  of the last pass thru the peaking filter
    spl1 = o_PeakingFilt_R.TRBJFilt_DoSamp(s_Tmp_Sine);
  )
  : //else if not playing, output silence
  (
    spl0 = 0;
    spl1 = 0;
  );

[matt42]

Hi JCJR,

I found the difference in our results. As far as I can tell the amount of delays to interpolate the signal should go 0 -> 0.25 -> 0.5 -> 0.75. However, I just noticed, you are doing the opposite, which as far as I can tell is "wrong" . However the results are much better for oversampling an EQ. Here's the full band response:



So, it's not perfect, but interestingly all the aliases (the crudely coloured in areas) center around the two green lines. These are the points that reflect back to nyquist on down sampling. So in effect all the errors centre around nyquist, which is actually pretty handy (EDIT: at least in the case of an EQ, not so much for a non linear effect).

Here's my peak filter test with that order of allpass delays:



Not bad at all, just increasing overshoot the closer to nyquist - from visual inspection of the graph, the error at nyquist is well under 2dB.

[JCJR]

I found the difference in our results. As far as I can tell the amount of delays to interpolate the signal should go 0 -> 0.25 -> 0.5 -> 0.75. However, I just noticed, you are doing the opposite, which as far as I can tell is "wrong" . However the results are much better for oversampling an EQ. Here's the full band response:

Thanks matt42
You may be correct that is the wrong way to do it. The simplest things can be points of confusion. It was a few months ago I tested it for quick'n'dirty allpass intersample peak detection and I think I did the tests right but haven't recently repeated the tests on the "order of allpass delays" and maybe had an error in the past. I think I first "mistakenly" thought the order would be as you say, 0 -> 0.25 -> 0.5 -> 0.75 and then fixed the error by reversing the direction.

Here is the reasoning for that order-- Let us name the last input sample we pumped thru the three allpass delays x[n - 1]. And the very last "4X upsample" input we fed into the EQ network is also x[n - 1]. The state of the 3 allpass delays and the state of the EQ network all reflect a last input of x[n - 1], before we are fed the next sample x[n].

Now we get called to process the now-current sample of x[n]. We first want interpolation 0.25 sample past the last x[n - 1] sample we previously fed the EQ network. The quarter-sample past x[n - 1] ought to be x[n - 1 + 0.25] = x[n - 0.75]? And then the half-sample position, x[n - 1 + 0.5] = x[n - 0.5], followed by the three-quarter sample interpolation point of x[n - 1 + 0.75] = x[n - 0.25], and then finally we feed the fourth oversample of this current cycle, which is x[n]. Doesn't that seem the correct order?