LP/HP Filter

DSP, Plug-in and Host development discussion.
SKyzZz
KVRist

Topic Starter

35 posts since 17 Jan, 2021

Post Mon Oct 11, 2021 3:49 am

Vokbuz wrote:
Mon Oct 11, 2021 2:21 am
Please show your "process" method.

For http://www.earlevel.com/main/2016/12/08 ... e-grapher/ try this:

a coefficients (zeros):
0.7921204236492381
1.0

b coefficients (poles):
1.0
0.20787957635076193

Code: Select all

      // Process methods
      // IIR processing
      T Process(T in, int channel) {
        assert(channel < MAXNUMCHANNELS);
        switch (mType)
        {
        case OnepoleHighpass:
          z1[channel] = in * a0 + z1[channel] * b1;
          in -= z1[channel];
          return in;
          break;
        case OnepoleLowpass:
          z1[channel] = in * a0 + z1[channel] * b1;
          return z1[channel];
          break;
        default: // for all two poles
          double out = in * a0 + z1[channel];
          z1[channel] = in * a1 + z2[channel] - b1 * out;
          z2[channel] = in * a2 - b2 * out;
          return out;
          break;
        }
      } // Process one sample

      void ProcessBlock(T** inputs, T** outputs, int nChans, int nFrames) {
        for (auto c = 0; c < nChans; c++)
        {
          for (auto s = 0; s < nFrames; s++)
          {
            Process(inputs[c][s], c);
          }
        }
      } // Process entire buffer

User avatar
Vokbuz
KVRist
242 posts since 24 Aug, 2014 from Moscow

Post Mon Oct 11, 2021 4:11 am

Your highpass works because you set lowpass coefficients and then subtracts those from input. X - LP = HP.
With this process method you need to set a0 and b1 for lowpass and highpass identically:
b1 = exp(-2.0 * M_PI * mFreq);
a0 = 1.0 - b1;

SKyzZz
KVRist

Topic Starter

35 posts since 17 Jan, 2021

Post Mon Oct 11, 2021 4:27 am

Vokbuz wrote:
Mon Oct 11, 2021 4:11 am
Your highpass works because you set lowpass coefficients and then subtracts those from input. X - LP = HP.
With this process method you need to set a0 and b1 for lowpass and highpass identically:
b1 = exp(-2.0 * M_PI * mFreq);
a0 = 1.0 - b1;
I set:

Code: Select all

 
       case OnepoleHighpass:
          b1 = exp(-2.0 * M_PI * mFreq);
          a0 = 1.0 - b1;
          break;
                              
        case OnepoleLowpass:
          b1 = exp(-2.0 * M_PI * mFreq);
          a0 = 1.0 - b1;
          break;
current: 3kHz - 6db/oct.
Screenshot 2021-10-11 at 15.23.43.png
etalon: 3kHz - 6db/oct.
Screenshot 2021-10-11 at 15.24.45.png
The only problem is that this is absolutely wrong.
You do not have the required permissions to view the files attached to this post.

earlevel
KVRian
595 posts since 4 Apr, 2010

Post Mon Oct 11, 2021 7:54 pm

Sorry, I didn't take time to read everyone's post, but...

The implementation of the one-pole is slightly different from that of the biquads. Sorry, but I make the website as I go, and at the moment of doing the one-pole filters, I chose to keep the simple one-pole evaluation from looking awkward in dealing with the signs. In other words, the b1 coefficient is the opposite sign that that implemented in the biquad. So it won't work right with the frequency response grapher either—the frequency response grapher is of the same orientation as the biquad calculator (which also can generae one-pole filters). So, flip the sign and all is well.

In other words, if you use the one-pole code to calculate the coefficiens and use its code to also run the filter, it works. If you use the biquad code to calculate coefficients and to run the filter, it works. If you use the biquad code to calculate the coefficients and run the file in the one-pole code...nope.

The latest calculator adds more first order types, phase display, and doubles as a frequency response grapher if you paste in the coefficients:

Biquad calculator v3
My audio DSP blog: earlevel.com

earlevel
KVRian
595 posts since 4 Apr, 2010

Post Mon Oct 11, 2021 8:59 pm

The key is in the second paragraph of the one-pole article:
Specifically, b1 = -e-2πFc and a0 = 1 – |b1|. (In the implementation, we’ll roll the minus sign in the summation into b1 and make it positive, for convenience.)
That is, a0 is 1 minus the magnitude of b1. In my biquad, my filter computation mimicked the natural order of the derivation of the equations, which subtracted the b (feedback) contributions. To do the same for the one-pole, my "process" evaluation would have read "return z1 = in * a0 - z1 * b1;" Fine so far, but that would put a minus sign on b1 and the a0 computation would either look like "a0 = 1.0 - -b1", or "a0 = 1.0 + b1", the former not pretty, the latter making it not intuitive that it matches the explanation of "1 minus the magnitude of b1". And the b1 calculation would have looked a little silly too, tacking on a minus sign that would be negated in the process calculation. I think either way would have confused someone, somewhere, I punted to go with the way that made the coefficient calculations look less stupid. :wink:

A one-pole filter
My audio DSP blog: earlevel.com

User avatar
Vokbuz
KVRist
242 posts since 24 Aug, 2014 from Moscow

Post Tue Oct 12, 2021 1:09 am

SKyzZz wrote:
Mon Oct 11, 2021 4:27 am
The only problem is that this is absolutely wrong.
If this is absolutely wrong, then HP should be wrong too. But you wrote before that it works fine. Are you sure that etalon is correct?

juha_p
KVRian
699 posts since 21 Feb, 2006 from FI

Post Tue Oct 12, 2021 2:18 am

Vokbuz wrote:
Tue Oct 12, 2021 1:09 am
SKyzZz wrote:
Mon Oct 11, 2021 4:27 am
The only problem is that this is absolutely wrong.
If this is absolutely wrong, then HP should be wrong too. But you wrote before that it works fine. Are you sure that etalon is correct?
Upper image has IIM/MZT type filter response (this LP filter is implemented in calculator v3 earlevel linked in his post (second last in filter type menu) ) and the lower image BLT type filter response.

SKyzZz
KVRist

Topic Starter

35 posts since 17 Jan, 2021

Post Tue Oct 12, 2021 3:11 am

juha_p wrote:
Tue Oct 12, 2021 2:18 am
Vokbuz wrote:
Tue Oct 12, 2021 1:09 am
SKyzZz wrote:
Mon Oct 11, 2021 4:27 am
The only problem is that this is absolutely wrong.
If this is absolutely wrong, then HP should be wrong too. But you wrote before that it works fine. Are you sure that etalon is correct?
Upper image has IIM/MZT type filter response (this LP filter is implemented in calculator v3 earlevel linked in his post (second last in filter type menu) ) and the lower image BLT type filter response.
Thanks for the answer. Do you have docs for BLT type filter or code examples?

juha_p
KVRian
699 posts since 21 Feb, 2006 from FI

Post Tue Oct 12, 2021 3:21 am

SKyzZz wrote:
Tue Oct 12, 2021 3:11 am
Thanks for the answer. Do you have docs for BLT type filter or code examples?
I did already give link to project in github and RBJ's Cookbook earlier and example code I gave earlier for LPF/HPF are BLT implementations.

earlevel
KVRian
595 posts since 4 Apr, 2010

Post Tue Oct 12, 2021 1:46 pm

SKyzZz wrote:
Mon Oct 11, 2021 4:27 am
The only problem is that this is absolutely wrong.
Too dicey to try to figure out exactly where you are at this point, so I'll just get to what I see:

I have to make some assumptions because there are no numbers on your charts, but I see you have two charts, one is apparently the code you are using, a one-pole lowpass. It looks pretty much as I would expect, within the limitation of not knowing the relative magnitude from your chart.

In the second you show is from "etalon"—I don't know what that is exactly, but it's certainly appears to be a one-pole, one-zero lowpass—the zero at Nyquist.

So it looks like you are comparing different filters, not sure if you expect them to look the same, but they won't.
My audio DSP blog: earlevel.com

earlevel
KVRian
595 posts since 4 Apr, 2010

Post Tue Oct 12, 2021 2:18 pm

OK, I gave the basic answer to you on my reply to you comment on my site, I'll give a quickie here for a first-order, BLT-derived filter (zero at Nyquist):

Code: Select all

    // calculate coefficients
    double norm, a0, a1, b1;
    double K = tan(M_PI * mFc);
    double norm = 1 / (1 + K * oneOverQ + K2);
    a0 = a1 = norm;
    b1 = (1 - 1 / K) * norm;
...
// process a sample
inline double Process(double in) {
    double out = in * a0 + z1;
    z1 = in * a1 - b1 * out;
    return out;
}
Note that the handler is the same as my biquad code, with the second order removed. You could get a little more efficient by observing that a0 is always the same as a1, if process will always be used only for lowpass. Of course, the same is true of biquad evaluation, depending on filter type.
My audio DSP blog: earlevel.com

SKyzZz
KVRist

Topic Starter

35 posts since 17 Jan, 2021

Post Wed Oct 13, 2021 6:41 am

earlevel wrote:
Tue Oct 12, 2021 2:18 pm
OK, I gave the basic answer to you on my reply to you comment on my site, I'll give a quickie here for a first-order, BLT-derived filter (zero at Nyquist):

Code: Select all

    // calculate coefficients
    double norm, a0, a1, b1;
    double K = tan(M_PI * mFc);
    double norm = 1 / (1 + K * oneOverQ + K2);
    a0 = a1 = norm;
    b1 = (1 - 1 / K) * norm;
...
// process a sample
inline double Process(double in) {
    double out = in * a0 + z1;
    z1 = in * a1 - b1 * out;
    return out;
}
Note that the handler is the same as my biquad code, with the second order removed. You could get a little more efficient by observing that a0 is always the same as a1, if process will always be used only for lowpass. Of course, the same is true of biquad evaluation, depending on filter type.
Yesterday, with the help of substitution, I got the desired slice, but for now, I'm thinking how to write it correctly in this case. As you can see, I took the code from BiquidLowPassFilter:



Code: Select all

        case BqLowpass:
          K = tan(M_PI * mFreq);  // 3000
          norm = 1 / (1 + K / mQ + K * K);          
          a0 = 0.17832624470902195; //K * K * norm; ?? 
          a1 = a0;
          a2 = 0;
          b1 = -0.6433475105819562; //2 * (K * K - 1) * norm; ?? 
          b2 = 0;
          break;
        }

Code: Select all

        default: // for all two poles
          double out = in * a0 + z1[channel];
          z1[channel] = in * a1 + z2[channel] - b1 * out;
          z2[channel] = in * a2 - b2 * out;
          return out;
          break;
        }
This is the working code that displays the filter that I need. Your code is essentially not working at all. The signal just disappears.

SKyzZz
KVRist

Topic Starter

35 posts since 17 Jan, 2021

Post Wed Oct 13, 2021 10:46 am

Solution:

var K = Math.tan(Math.PI * Fc / Fs);
...
case "lowpass 1p1z":
norm = 1 / (1 / K + 1);
a0 = a1 = norm;
b1 = (1 - 1 / K) * norm;
a2 = b2 = 0;
break;

earlevel
KVRian
595 posts since 4 Apr, 2010

Post Wed Oct 13, 2021 12:31 pm

SKyzZz wrote:
Wed Oct 13, 2021 6:41 am
This is the working code that displays the filter that I need. Your code is essentially not working at all. The signal just disappears.
From your solution that followed, it appears that you missed that I gave a process routine that omits the second pole/zero, hence there is no a2/b2 to zero. :wink:
My audio DSP blog: earlevel.com

Return to “DSP and Plug-in Development”