## When to choose a ZDF

DSP, Plug-in and Host development discussion.
mystran
KVRAF
4981 posts since 12 Feb, 2006 from Helsinki, Finland
stratum wrote:
Z1202 wrote:Maybe you should use voltage law in a different way, like I proposed? This way the cycles will be implicitly covered by the absolute variables. And I did detail an example showing how this happens (your ABCDEF circuit).
Step 4 ('Write Kirchhoff's equations for the nodes') implies the usage of the current law. The voltage law always yields 3 equations (one being redundant) for that circuit because the corresponding graph has 3 cycles. The algorithm needs to proceed from that point and needs to figure out what to do about the third equation.
You get a redundant equation because your circuit is floating. Whether you use voltage or current laws has nothing to do with it. Pick a ground node and the problem goes away.
If you'd like Signaldust to return, please ask Katinka Tuisku to resign.

stratum
KVRAF
1841 posts since 29 May, 2012
Z1202 wrote:
stratum wrote:
Z1202 wrote:Maybe you should use voltage law in a different way, like I proposed? This way the cycles will be implicitly covered by the absolute variables. And I did detail an example showing how this happens (your ABCDEF circuit).
Step 4 implies the usage of the current law. The voltage law always yields 3 equations (one being redundant) for that circuit because the corresponding graph has 3 cycles. The algorithm needs to proceed from that point and needs to figure out what to do about the third equation.
Again. I propose to abandon the way how you use the voltage law, and address it differently, in the way I showed. The absolute voltages are the thing which does the trick here (lets you do without your version of voltage law, which becomes implicitly contained in your choice of variables).
To me this sounds like 'abandon the voltage law', which is OK, apparently it's not necessary at all for nodal analysis. The current law implicitly covers that case as your analysis shows.
~stratum~

stratum
KVRAF
1841 posts since 29 May, 2012
You get a redundant equation because your circuit is floating. Whether you use voltage or current laws has nothing to do with it. Pick a ground node and the problem goes away.
That's not the problem. I wasn't doing nodal analysis properly, you have told it, it's mesh analysis+nodal analysis mixed together.
~stratum~

Z1202
KVRian
958 posts since 12 Apr, 2002
mystran wrote:
Z1202 wrote:It's not that I'm trying to claim that MNA doesn't work. And as long as it's happening fully under the hood of a software piece, it doesn't matter. It's that I'm regularly hearing of people doing this by hand, and I can't figure out the reason why, except maybe inavailability of descriptions of other "more human-friendly" methods
I think people just have irrational fear of code-generators, even though there are problems (like this one) where anything else is just utter non-sense.
Often it helpful to follow the process of solving the equations as this can give you insights about how the system works.

As to my "time-varying" comment (again, maybe you didn't see the edit).

Imagine the capacitor equation: I=C * dV/dt.
We rarely think that this equation is valid only in the time-invariant case. My guess would be that nodal analysis will assume that C is constant. So the result will not be applicable to the case where it's not.

I'm not sure, but I think one can construct cases where even a varying resistor can have a similar effect.

OTOH, I have a feeling that that "natural choice of state variables" is not subject to this problem (or at least less likely to be).

Edit: or maybe it isn't. The truly natural choice of the variable would be the capacitor charge Whatever.

Edit2: still, I vaguely remember an article from simulanalog.org, where choosing currents as state variables in a capacitor-based system turned out to be not such a good idea. IIRC it even could lead to a zero division. Maybe I should refresh my memory. OTOH, I never had this problem with naturally-chosen state variables.
Last edited by Z1202 on Thu Nov 30, 2017 2:21 pm, edited 2 times in total.

Z1202
KVRian
958 posts since 12 Apr, 2002
stratum wrote:To me this sounds like 'abandon the voltage law', which is OK, apparently it's not necessary at all for nodal analysis. The current law implicitly covers that case as your analysis shows.
Actually it's not the current law that covers it, but a different choice of variables. If you think about it, working with absolute voltages leaves no other choice but to have the total voltage difference across a loop zero.

stratum
KVRAF
1841 posts since 29 May, 2012
Z1202 wrote:
stratum wrote:To me this sounds like 'abandon the voltage law', which is OK, apparently it's not necessary at all for nodal analysis. The current law implicitly covers that case as your analysis shows.
Actually it's not the current law that covers it, but a different choice of variables. If you think about it, working with absolute voltages leaves no other choice but to have the total voltage difference across a loop zero.
But there are 3 loops in that graph and the voltage rule says write equations for all of them. Regardless of whether we use absolute or relative choice of variables, usage of the rule requires that we deal with that fact: we have those 3 equations when the rule is followed.
~stratum~

Z1202
KVRian
958 posts since 12 Apr, 2002
stratum wrote:
Z1202 wrote:
stratum wrote:To me this sounds like 'abandon the voltage law', which is OK, apparently it's not necessary at all for nodal analysis. The current law implicitly covers that case as your analysis shows.
Actually it's not the current law that covers it, but a different choice of variables. If you think about it, working with absolute voltages leaves no other choice but to have the total voltage difference across a loop zero.
But there are 3 loops in that graph and the voltage rule says write equations for all of them. Regardless of whether we use absolute or relative choice of variables, usage of the rule requires that we deal with that fact: we have those 3 equations when the rule is followed.
Forget the voltage rule. It's implicitly contained in the variable choice and you don't need to write those equations at all.

stratum
KVRAF
1841 posts since 29 May, 2012
Forget the voltage rule. It's implicitly contained in the variable choice and you don't need to write those equations at all.
Yes, that seems to be the case. It's redundant and not needed at all for nodal analysis.
~stratum~

JCJR
KVRAF
2342 posts since 17 Apr, 2005 from S.E. TN
Hi Max, Vadim, whoever else might be interested-- I tested Max's first order code, named "BLT DF" filters, please suggest a better name if that isn't a good one. Compared against Vadim's Trapezoidal first-order filters.

Under my simple test if both do not have identical transfer functions, then at least they are VERY similar. I test compared LP, HP and AP at Fc = 100 Hz and also Fc = 5 kHz. Six tests. Test input was stereo 10 sec duration, -6 dB peak sine wave sweep, 20 to 20 kHz at 44.1 k samplerate.

All six tests null at all locations to -90.3 dB, probably because I forgot to disable dither in Reaper when rendering the test files. Am guessing it would null at or nearer to -inf if done without dither but am too lazy to repeat the tests.

At both 100 Hz Fc and 5 kHz Fc, the LP and HP -3 dB points are "very close to the same" but I don't measure exactly identical. They may be exact, with the error due to my testing somehow. At both 100 and 5 kHz the -3 dB points of LP and HP seemed "a tiny bit different from each other" but it could be pilot error. Close enough for rock'n'roll anyway. IOW, both Vadim's and Max's filters seem to perform identically, and both MAY have the identical VERY SLIGHT disagreement of LP -3 dB frequency versus HP -3 dB frequency. Or maybe I should have tested it some better way for more accurate result..

For the Fc=5 kHz LP and HP, I measured -3 dB in the ballpark of 4900 Hz but it is ambiguous to measure-- I just counted the number of samples in 8 wave periods surrounding -3 dB, but it might be closer than that to 5 kHz (or exactly 5 kHz for all I know). Anyway even if it is really "about 4900 Hz" then that is much closer than required by rock'n'roll.

I do think Max's first order code works better than some other "old fashioned" first order codes I tried in the past, begging the question why professors would put sucky code in their books. But because at least casually both kinds of filter seem to behave so similer, it comes full-circle back to, at least in the first-order case, is there a reason to choose one or the other?

Maybe one of the filters works better on certain kinds of transients? Or maybe one filter is better-behaved when swept? Someone else would be better qualified to test that if so interested. The code is in Reaper jsfx, using jsfx "name space objects" which are similar enough to C that it shouldn't take long to convert to C/C++ should someone desire. I stopped using C++ for awhile because simple things are less frustrating in jsfx. If I'm gonna get frustrated doing code it would be like having a job.

Code: Select all

``````desc: FirstOrdFiltCompare jcjr First Order Filter Test
//version: 0.1
// This effect Copyright (C) 2017 and later James Chandler Jr
//   http://www.errnum.com
//tags: First Order Filters
//author: James Chandler Jr

//A quick comparison test of first order LP, HP and AP filters
//One a trapezoid "zdf" filter adapted from Vadim Zavalishin
//The second a "BLT Direct Form" [z-1] filter adapted from Max Mikhailov suggestions
//
//Both filters process the left channel of the input file (identical input for both filters)
//Trapezoid filter writes to Left channel output
//BLT DF filter writes to Right channel output
//
//Test setup lines are above the line "@slider" near the bottom of code
//Uncomment one pair of lines to do a test, or make your own lines to test in another fashion
//
//I used a 10 second duration stereo sine wave sweep
//20 to 20 kHz at 44.1 k samplerate and -6 dB peak amplitude
//
//I tested LP, HP and AP at Fc = 100 Hz and Fc = 5 kHz
//Both filters tested "identical" but perhaps some other test could find a difference, dunno

@init
//generic utility functions

//should only call js_malloc within js @init section
JS_MALLOC_NO_CLEARMEM = 0;
JS_MALLOC_CLEARMEM = 1;
next_js_malloc = 0;
function js_malloc(a_size, a_clearmem)
local (l_result)
(
l_result = next_js_malloc;
(a_clearmem != 0) ?
memset(l_result, 0, a_size);
next_js_malloc += a_size;
l_result;
);

//First order filtertype named consts for FirstOrdTrapezoidFilter and FirstOrdBltDfFilter
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

//First order trapezoidal filter object
//Code adapted from Vadim Zavalishin's book "The Art of VA Filter Design"
//
//FiltType  : Use one of the filtertype named consts to set the return value of FirstOrdTrapezoidFilter_DoSamp()
//          : However, all values are simultaneously accessible after calling DoSamp() by reading properties--
//          : 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);
);

//only call after filter was created via _Init
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);
);

//only call after filter was created via _Init
//Returns the filtered sample
function FirstOrdTrapezoidFilter_DoSamp(a_InSamp)
local (l_v, l_result)
(
//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;
);

//First order BLT Direct Form filter object
//Code adapted from Max Mikhailov's suggestions:
//
//The BLT LP expresses in direct form (i.e. numerator/denominator, transfer-function form) as:
//b = [1-k, 1-k]/2; a = [1, k];
//Thus HP is:
//b = [1+k, -1-k]/2; a = [1, k];
//And AP is:
//b = [-k, 1]; a = [1, k];
//
//Where k = cos(w)/(1 + sin(w));
//and w = 2*pi*Fc/Fs;
//
//In practice they are usually implemented in opposite fashion:
//first your get the AP and then LP = (IN + AP)/2; HP = (IN - AP)/2;

//FiltType  : Use one of the filtertype named consts to set the return value of FirstOrdBltDfFilter_DoSamp()
//          : However, all values are simultaneously accessible after calling DoSamp() by reading properties--
//          : TheFilter.lp, TheFilter.hp, TheFilter.ap_A, TheFilter.ap_R
//a_FiltFC  : Filter center frequency in Hz
//a_SampRate: Samplerate of filter
function FirstOrdBltDfFilter_Init(a_FiltType, a_FiltFC, a_SampRate)
local (l_w)
(
this.FT = a_FiltType;
this.SR = a_SampRate;
this.Nyquist = floor(this.SR * 0.5);
this.FC = min(a_FiltFC, this.Nyquist - 1);

this.FBFF = 0.0; //abbreviation for the [z-1] feedback+feedforward calculated value
this.lp = 0;
this.hp = 0;
this.ap_A = 0;
this.ap_R = 0;

//calculate coefficient
l_w = (2 * \$pi * this.FC) / this.SR;
this.k = cos(l_w) / (1 + sin(l_w));
);

//only call after filter was created via _Init
function FirstOrdBltDfFilter_SetFC(a_FiltFC)
local (l_w)
(
this.FC = min(a_FiltFC, this.Nyquist - 1);
l_w = (2 * \$pi * this.FC) / this.SR;
this.k = cos(l_w) / (1 + sin(l_w));
);

//only call after filter was created via _Init
//Returns the filtered sample
function FirstOrdBltDfFilter_DoSamp(a_InSamp)
local (l_result)
(
//From Max Mikhailov suggestions
//this.ap_R = (a_InSamp * -1.0 * this.k) + this.FBFF;
this.ap_R = this.FBFF - (a_InSamp * this.k);
this.FBFF = a_InSamp + (this.ap_R * this.k);
this.lp = 0.5 * (a_InSamp + this.ap_R);
this.ap_A = -1.0 * this.ap_R;
this.hp = 0.5 * (a_InSamp + this.ap_A);

(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;
);

pdc_bot_ch = 0;
pdc_top_ch = 2;

//Instantiate two test filters

//Allpass test 100 Hz
//o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 100, srate);
//o_1stOrdBltDfFilt.FirstOrdBltDfFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 100, srate);

//Lopass test 100 Hz
//o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_LOPASS, 100, srate);
//o_1stOrdBltDfFilt.FirstOrdBltDfFilter_Init(FIRST_ORD_FILTTYPE_LOPASS, 100, srate);

//Hipass test 100 Hz
//o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_HIPASS, 100, srate);
//o_1stOrdBltDfFilt.FirstOrdBltDfFilter_Init(FIRST_ORD_FILTTYPE_HIPASS, 100, srate);

//Allpass test 5 kHz
//o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 5000, srate);
//o_1stOrdBltDfFilt.FirstOrdBltDfFilter_Init(FIRST_ORD_FILTTYPE_ALLPASS_RET, 5000, srate);

//Lopass test 5 kHz
//o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_LOPASS, 5000, srate);
//o_1stOrdBltDfFilt.FirstOrdBltDfFilter_Init(FIRST_ORD_FILTTYPE_LOPASS, 5000, srate);

//Hipass test 5 kHz
o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_Init(FIRST_ORD_FILTTYPE_HIPASS, 5000, srate);
o_1stOrdBltDfFilt.FirstOrdBltDfFilter_Init(FIRST_ORD_FILTTYPE_HIPASS, 5000, srate);

@slider
//no sliders declared

@block
//no block processing

@sample
//Apply test filters
//Use left channel input to feed both the trapezoid filter and the BLT DF filter
//Send trapezoid filter output to left channel and send BLT DF filter output to right channel
//So that the two transfer functions can be compared with a audio editor program or whatever
InSamp = spl0;
spl0 = o_1stOrdTrapezoidFilt.FirstOrdTrapezoidFilter_DoSamp(InSamp);
spl1 = o_1stOrdBltDfFilt.FirstOrdBltDfFilter_DoSamp(InSamp);

@gfx
//no graphics displayed``````

JCJR
KVRAF
2342 posts since 17 Apr, 2005 from S.E. TN
Re the "slight differences" between LP and HP Fc mentioned in earlier message-- DOH, it may just be a measurement oddity of the filter phase shifts, or some other dumb mistake.

Max M.
KVRist
262 posts since 20 Apr, 2005 from Moscow, Russian Federation
JCJR

I'll try to reply to most of your questions later, but meanwhile for the differences in tests: I suppose it's just a result of single-precision fp accuracy (I guess Reaper JS engine uses 32-bit floats for all this stuff). I think I'd be better to give "inverted" - more suitable for fp32 - coefficients (as in `k -> (1-k')`) right away (plus maybe a few tips on other minor numerical-stability tricks) . Well, never mind, I guess I'll put a more practical (adopted for a single-precision) pseudo-code in the later answer...

JCJR
KVRAF
2342 posts since 17 Apr, 2005 from S.E. TN
Max M. wrote:JCJR

I'll try to reply to most of your questions later, but meanwhile for the differences in tests: I suppose it's just a result of single-precision fp accuracy (I guess Reaper JS engine uses 32-bit floats for all this stuff). I think I'd be better to give "inverted" - more suitable for fp32 - coefficients (as in `k -> (1-k')`) right away (plus maybe a few tips on other minor numerical-stability tricks) . Well, never mind, I guess I'll put a more practical (adopted for a single-precision) pseudo-code in the later answer...
Thanks Max

Not to worry about the single-precision thang. Maybe sometimes Reaper can use double float audio data files, but I always use float 32 audio. But for the testing (generating the source sine sweep and rendering the result) the float32 (or float64 internal) was imported/exported via stereo 16 bit int files. Because I'm still using the last update of cool edit pro for stereo editing, which is fine for most uses but has a weird idea about the format of float audio files. Well, I also use Audacity a lot for encode and such, which is fine with 32 bit float files but isn't as good "analytic stereo editor" as old cool edit pro. I never updated it to adobe audition because couldn't find many folks who thought the audition versions brought much new to the table for the extra money.

Anyway, the jsfx is ALWAYS float64. It doesn't know any other variable type. You can do int or whatever math in a jsfx code, but all the variable memory is just a big array of double floats, and all the math is always done in double floats.

The "minor differences" I found between HP and LP cutoff and such, really are minor, only mentioned for sake of completeness and it could easily be some mistake I made in measurement. In the sine sweep test, your BLT DF and Vadim's Trapezoid code gives effectively the same results so far as I can see. Maybe differences would be noticed rapidly sweeping the filters or whatever, but I doubt if I could measure that accurately and probably won't be doing much filter sweeping, so it doesn't much matter to me. Its just interesting that the "gross behavior" of the two codes seems so similar.

matt42
KVRian
1057 posts since 9 Jan, 2006
Well I would look at it in terms of ZDF (with bilinear integrator) is just another topology of the bilinear transform. There isn't really such a thing as "BLT" filter in and of itself - typically direct forms were used as the whole system was transformed. With ZDF the integrators are transformed and the feedforward/feedback paths are resolved. But all these topologies (direct forms/zdf) are underpined by the BLT.

sorry if this gets TLDR We could take a Moog filter and take the BLT of the system, or take the BLT of the one pole filters and wire four of those together including zero delay feedback path, or even go one step deeper and take the BLT of the integrators wire those up to make one pole low pass filters and wire those together with the zero delay feedback path. All three versions would be identical (for LTI cases) because they are all essentially different typologies of BLT. I suppose it's just "luck" that the BLT preserves the frequency phase responses of the mapping otherwise I guess the above wouldn't hold true.

JCJR
KVRAF
2342 posts since 17 Apr, 2005 from S.E. TN
Thanks Matt. Ya I knew the trapezoidal filter is BLT-based. Was just trying to find a reasonable name for the variant which Max described. Am open to suggestions. Early in the thread Max had asserted that [z-1] filters can have exactly the same transfer function as the trapezoidal filters, but I had not seen "real good-behaving" examples of all three first order filters and asked him as curiosity question if he could suggest one which works "identical" to the trapezoid first order filters.

The allpass "core" of that one looks similar or identical to many allpass filters including the one JOS wrote about long ago for sample-delay interpolation. I used that one in my collection of "samplerate conversion objects" long ago. When looking for a good 1st order allpass a couple years ago for compressor code, found several versions of that one on the interwebz, with several "seemingly contradictory" tuning methods, causing wonder if in fact all of the examples could have even worked as claimed.

If I read zolzer about basing a filter family on allpass, was too dense to appreciate what he was getting at, or maybe the explanation was in abstract math that passed thru the eyes and out the back of the skull. Just saying, did not recall reading of such an easy way to get all the first-order filters, well-behaved, based on allpass, but it makes sense.

In comparing the two filters, didn't know what to expect. Would probably have considered response curves "pretty close" as quite good. Something that nulls at -24 or -48 dB not bad at all. It is just guessing, but I think if I bothered to put the sweep generator in my test program and did the nulling in the test program, odds are that the two filters would null near -inf dB, which is a lot better than "quite good". I think my -90.3 dB comparison null is because of accidentally dithering 16 bit output test files, but even if they really do differ by -90 dB, that is awfully close for government work! (In defense of old Cool Edit Pro, the reason its float files are odd, is because it was written before the Float WAV format was formalized.)

matt42
KVRian
1057 posts since 9 Jan, 2006
Hi JCJR, sorry if I was explaining the obvious