Book: The Art of VA Filter Design 2.1.2

[Z1202]

I was thinking of how the naive step in the filter I'm using works, the b1 turns the naive integrator to a leaky integrator, and the a0 is the DC gain adjustment. I thought it might be possible to apply a0 and b1 in the same way to the trapezoidal integrator, and it is, just not in the way I expected. I'm hoping to make a 1p lowpass with improved accuracy in both frequency and phase near nyquist, but I do realize the tradeoffs involved. I'm not adept with figuring out the transfer functions like you did, but I'll study what you said and see what I can do!

The problem is that b1 is embedded into the naive integrator. It can be expressed as feedback path instead:
u[n]=a0*x[n]
y[n]=u[n]+b1*y[n-1]=u[n]+( 1 - (1-b1) )*y[n-1]= ( u[n] - (1-b1)*y[n-1] ) + y[n-1]
Thus you get a naive integrator whose input is u[n] - (1-b1)*y[n-1]. Then in principle you can replace it with the trapezoidal one, but (depending on whether you feed back y[n], which is delayless feedback, or y[n-1]) you can get exactly the normal BLT TPT filter, most probably (too lazy to think about the details now) with offset cutoff and gain.

However, I'm not sure what you're aiming at. Maybe you first should decide what kind of frequency response you want to obtain. The frequency response is defined solely by the transfer function. Then you can decide on the topology. Also, there are not that many different transfer functions for a 1-pole filter. They all have the familiar form:
H(z)=(b0+b1*z^-1)/(1+a1*z^-1) (I've switched the usage of a and b coefficients to a more common one). The cutoff is determined by the a1 coefficient, whereas the b0 and b1 coefficients define the type of the filter response (lowpass, highpass etc). Note: to keep the response consistent for different cutoffs, b0 and b1 also generally should depend on a1.

Edit: roughly speaking, b0 and b1 define the mixture between the lowpass and highpass components of the output (this is more intuitively visible in the continuous-time domain, but the idea is totally similar) and that is all there to it. You can only choose the amount of mixing. In digital world there is some variation in regards to what you can choose as the "basis" lowpass and highpass signals, but in the end you get the same mixtures. So e.g. naive lowpass is the same as the BLT lowpass with some BLT highpass signal mixed in (the mixture coefficients may depend on the cutoff, I don't want to think about it too much now ).

Edit2: what you could do is a simple experiment. Choose a1 corresponding to some typical cutoff value, e.g. a1=-0.9. Under this limitation, the frequency response will be defined solely by the balance between b0 and b1 (because multiplying both b0 and b1 by the same value solely changes the total gain of the filter). Thus let b0=cos(alpha), b1=sin(alpha). By simply varying alpha you can obtain all possible 1-pole response shapes (for a given cutoff). This will include various versions of lowpass, highpass, low-shelving, high-shelving and allpass. The available shapes will vary with the cutoff choice though, as a1 is changed from values close to -1, corresponding to low cutoffs to 0 and even further to 1, corresponding to cutoffs close to Nyquist, where note that you need |a1|<1 to keep the filter stable.

[camsr]

Do you have a graphic of the topology for the basic 1-pole you described? Or can any of them be used?

I've made some progress on my filter, I am using a two-sample integrator (I can't call it trapezoidal) which is adjustable in tandem with the cutoff. So far, the FR is very promising, but I will have to check phase as well. It attempts to match the analog FR via adjustment of the integrator.

[Z1202]

Do you have a graphic of the topology for the basic 1-pole you described? Or can any of them be used?

As long as we are talking purely of transfer functions, the topology is irrelevant, it doesn't make a difference for the static response (with the exception of precision loss issues which can vary from topology to topology, but are usually not critical for 1-poles).

I've made some progress on my filter, I am using a two-sample integrator (I can't call it trapezoidal) which is adjustable in tandem with the cutoff. So far, the FR is very promising, but I will have to check phase as well. It attempts to match the analog FR via adjustment of the integrator.

Would be interesting to see your results (I guess it would be more appropriate to share them in your own thread).

[camsr]

I would like to compare to the analog prototype, but I'm not sure how. I have Octave (like matlab) but I'm not familiar with the signal functions. Maybe I can post some screenshots from SPAN later.

[1eqinfinity]

Hello everyone!
I'm having fun with linear multimode lader (chapter 4.4 of The Art). I'm stuck at the notch, and either I got some idea wrong, or my algebra is wrong (however I double checked it).
So what I think I should be doing:
On the left there's some desired transfer function H that is a ratio of polynomials.
On the right I have

Code: Select all

a0*(1 + s)^4 + a1*(1 + s)^3 + a2*(1 + s)^2 + a3*(1 + s) + a4
----------------------------------------------------------------------
k + (1 + s)^4

that is a linear combination of several equations 4.5 for different n (page 66).
Then
- Proceed as if k = 0.
- Substract 1 from all s, so that s+1 on the right becomes s, and e.g. s^2 on the left becomes (s-1)^2;
- Bring H on the left to the form that has the same denominator k + (1 + s)^4 as on the right and write it in the decreasing order of powers of s.
- The coefficients of corresponding powers of s on the left are those that I'm after.

LP, HP and BP are working fine (except that if k > 3.1 my implementation occasionally explodes somewhere in the 3kHz-5kHz range).
On page 68 Vadim proposes (s`^2 + w)^2 / (s + 1)^4 as a transfer function for the notch. Solving that I get really weird stuff (see the code please). Using analog frequency for the value of w yields a lowpass (\O_o/) with a somewhat flatter slope. Using angular frequency (which I think what I should do, so w = analog_freq * pi / sample_rate) creates extreme overdrive because of the form of the coefficients I got after solving.

What am I doing wrong?

Here's my code:

Code: Select all

double Process(double in) {
		auto gDen = g / (1.0 + g); // g = analog_freq * pi / sample_rate
		auto g2 = g * g;
		auto g3 = g2 * g;
		auto G = g2 * g2;

		auto S = g3 * s1 + g2 * s2 + g * s3 + s4;

		// feedback shaping
		auto hp = S - fbsS;
		fbsS = fbsS + 2.0 * hp * fbsF; // fbsF = 40.0 * pi / sample_rate;
		S = hp;
		
		assert(isfinite(S));

		auto u = (in - k * S) / (1 + k * G);

		double v; // temp variable
		v = (u - s1) * gDen;
		auto y1 = v + s1;
		s1 = y1 + v;

		v = (y1 - s2) * gDen;
		auto y2 = v + s2;
		s2 = y2 + v;

		v = (y2 - s3) * gDen;
		auto y3 = v + s3;
		s3 = y3 + v;

		v = (y3 - s4) * gDen;
		auto y4 = v + s4;
		s4 = y4 + v;

		double r = 0.0;

		switch (passType) {
				default:
				case FLT_LP:
					r = y4;
					break;
				case FLT_HP:
					r = 1.0 * u - 4.0 * y1 + 6.0 * y2 + -4.0 * y3 + y4;
					break;
				case FLT_BP:
					r = y2 - 2.0 * y3 + y4;
					break;
				case FLT_BR:
					v = 1.0 / (f * f); // f = analog frequency
					r = v * u + 4.0 * v * y1 + (2.0 + 6.0 * v) * y2 - 4.0 *(1.0 + v) * y3 + (v + 2.0 / f + 1.0)*y4;
					break;
			}

		assert(isfinite(r));
		return r;
		}

[Z1202]

On page 68 Vadim proposes (s`^2 + w)^2 / (s + 1)^4 as a transfer function for the notch.

Not sure what was I smoking back then Looking at it now, I think it only works at w=1. Need to fix that in the book, thanks!

Edit: although maybe I have been smoking this: "Notably, these modes are not necessarily perfectly matching their descriptions, the parameters may have some weird side effects.". At any rate they work properly at w0=1

IIRC I didn't try out all the modes, just wrote them out of my head, apparently making some mistakes. My apologies.

Edit: as for the explosion, I wouldn't expect it to happen at 3.1 maybe at 3.99 but not 3.1. Although feedback shaping is somewhat inconvenient to analyse. Does it still explode if you remove the feedback shaper? Also there could be a bug in the shaper (which I didn't check for in your code).
Edit2: it seems that feedback hp doesn't affect G in any way, which is strange

[1eqinfinity]

Thanks for the reply and for your book too! Nothing to apologize for.
So the workflow was correct?
Because I've found this thread http://www.kvraudio.com/forum/viewtopic.php?t=207647 and mystran proposes to combine 2-pole LP and HP into a single transfer function. So I could put (1 + s^2) / (1 +s)^2 on the left, right?

Edit: yes, it explodes without the feedback shaper too.

Edit 2: maybe I'm setting g wrong. I do

Code: Select all

g = analog_freq * pi / sample_rate

[Z1202]

Because I've found this thread http://www.kvraudio.com/forum/viewtopic.php?t=207647 and mystran proposes to combine 2-pole LP and HP into a single transfer function. So I could put (1 + s^2) / (1 +s)^2 on the left, right?

Could you point to the exact message? I'm not sure what specifically you're referring to and don't have time now to read through the entire thread.

Edit: yes, it explodes without the feedback shaper too.

Edit 2: maybe I'm setting g wrong. I do

Code: Select all

g = analog_freq * pi / sample_rate

Hmmm, strange. Otherwise upon the first look the code seems correct. The value of g shouldn't affect the explosion, unless it's very large or negative.

[1eqinfinity]

With w = 1 I get a proper notch, thanks!

I'm refering to this mystran's sentence in his very first message:

For notch filter you need to do something more clever, and add hipass and lowpass responses together (at least 2-pole for each such that they happen to cancel each other near the cutoff).

I derived the H of notch as Hlp^2 + Hhp^2

[Z1202]

Actually you should be able to construct the modes yourself. E.g. the idea behind the notch is that s^2+1 has a zero at s=j. The square (s^2+1)^2 ensures that in the limits at s=0 and s=infty we get 1. You could try to follow the similar logic for all other modes and even construct your own. Maybe I should have given more explanation in the book.

[Z1202]

With w = 1 I get a proper notch, thanks!

I'm refering to this mystran's sentence in his very first message:

For notch filter you need to do something more clever, and add hipass and lowpass responses together (at least 2-pole for each such that they happen to cancel each other near the cutoff).

I derived the H of notch as Hlp^2 + Hhp^2

Seems to me exactly the same. 2-pole highpass has s^2 and 2-pole lowpass has 1 in the numerator. However you still need to take care of the limiting values at 0 and infinity somehow (if you just use s^2+1 you'll get a zero at infinity).

[Z1202]

Here's another notch s^4-1. Maybe should add it to the book. Again the idea is that at s=j we have a zero and at 0 and infty we have the same as the denominator.

[1eqinfinity]

The square (s^2+1)^2 ensures that in the limits at s=0 and s=infty we get 1.

Sorry, I don't quite get it. So we map the [0, +inf) freq range to the unity circle. If we let the H for notch be (1 + s^2) / (1 + s)^2 it will also have 1 at the inf. If I'm not mistaken, (inf^2 + 1) / (inf + 1)^2 is 1.

[1eqinfinity]

Funny thing, If I boldly solve for (Hlp^2 + Hhp^2)^2 I get a peak mode with some beareable freq response distortion to the right of the f, but only with nonzero feedback.
i.e H = (1+s^2)^2 / (1+s)^4

Edit: with zero fb it's a gentle notch.

[Z1202]

The square (s^2+1)^2 ensures that in the limits at s=0 and s=infty we get 1.

Sorry, I don't quite get it. So we map the [0, +inf) freq range to the unity circle. If we let the H for notch be (1 + s^2) / (1 + s)^2 it will also have 1 at the inf. If I'm not mistaken, (inf^2 + 1) / (inf + 1)^2 is 1.

Yes, but for the ladder the denominator is (1+s)^4, not (1+s)^2. Unless you cancel it with an additional (1+s)^2 factor in the numerator (resulting in yet another notch filter).

Also I'm not sure why you referred to the unit circle...