LP/HP Filter

[SKyzZz]

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

[Vokbuz]

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]

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.

[earlevel]

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

[earlevel]

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.

A one-pole filter

[Vokbuz]

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]

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]

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]

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]

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.

[earlevel]

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.

[SKyzZz]

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]

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]

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.