I wanted to give your filter a try. I used ACD to generate the input values. These are:
analog coefficients:
a2 = 0
a1 = 0
a0 = 3947841760
b2 = 1
b1 = 125663.7061
b0 = 3947841760
IIR digital coefficients:
b0 = 0.06422539
b1 = 0.12845078
b2 = 0.06422539
a1 = -0.986290852
a2 = 0.243192411
Results:
OPTIMIZED IIR FILTER COEFFICIENTS:
b0 = 0.064543481
b1 = 0.173126494
b2 = -0.010531462
a1 = -1.049711599
a2 = 0.276921804
The above values were obtained when the optimization was carried out between 1kHz and 36kHz. I moved these up because of the higher sampling rate. I can only assume that your starting point was too far off for the SOLVER to get a stable result.
This particular filter shows why this is mainly (for audio) a problem with 44.1k and 48k SR. Here, when we have moved up to 96k SR the deviations in the "unoptimized" IIR digital filter are getting to 1 or 2 dB only by 20kHz or higher. It's been well demonstrated that even if you claim you can hear up there (not many can) there is no sensitivity to changes in the response of only 1 or 2 dB. With the lower sampling rates, the deviations move down into the range of hearing where they can start to be noticed. Hardware and software that can do 44.1 and 48k SR IIR filters is now getting more common and inexpensive but the same can not be said for 96kHz SR. Since I want to use a more affordable solution, I am left to tweak the IIR filter in order to get it to do what I want.
Pics attached:
Attachments
Last edited:
um, isn't the answer in the phrase "inside the unit circle"? - if you are talking biquads the root calculation is trivial, "inside the unit circle" another easy calc - only need something indirect Bistritz stability criterion - Wikipedia, the free encyclopedia "Routh-Hurwitz like" for higher order polynomials that may have no closed form root equations
Fair enough, I downloaded the Python optimization package maybe I'll play around with it. I tried some simulated annealing a while back but the initial guess was still too important. As you see from the .0001dB RIAA one can get carried away with obsessive accuracy.
EDIT - thanks again, yes you do need to give it a better hint.
EDIT - thanks again, yes you do need to give it a better hint.
Last edited:
Just a follow up. I dragged out my simulated annealing Python code and yes the answer in post 21 is a good local minima, certainly good enough for any reasonable purpose.
It took 7 or 8 tries for me to find the global minima which is better but sort of needless resolution. To jcx's point, I don't use the coefficients as input to the search but use the poles and zeros in both s and z domain. It then becomes trivial to search only inside the unit circle.
I might also add that that with finite poles and zeros there is no possible exact fit but (I can't prove this) the global minimum seems to be the only one with the most ripples in the residual error. See my pic of the 48k RIAA solution.
It took 7 or 8 tries for me to find the global minima which is better but sort of needless resolution. To jcx's point, I don't use the coefficients as input to the search but use the poles and zeros in both s and z domain. It then becomes trivial to search only inside the unit circle.
I might also add that that with finite poles and zeros there is no possible exact fit but (I can't prove this) the global minimum seems to be the only one with the most ripples in the residual error. See my pic of the 48k RIAA solution.
Last edited:
I agree that if I had the transfer function in the "poles and zeroes" factored form it would be very easy to determine if all poles were inside the unit circle. But don't have that. What I have is the coefficients, which I can generate for a wide variety of filters, basically all the filters that one would want to use for a loudspeaker crossover. It is POSSIBLE to factor the resulting polynomials into the factored form, but that is more coding effort than I want to do at this time. Since the optimization can be run on the coefficients directly that is the easiest to implement without adding a bunch more code that may or may not work well. I will keep it in mind, however, for later.
Also, I don't see a real need to find the global minimum. How do you know you have that anyway? Practically speaking one needs the filter to be "good enough" and there certainly is a lot of wiggle room in how you define that. For instance, if you simply change the maximum frequency of the interval over which the optimization is conducted the "optimized" coefficients will change because the problem has been reformulated.
As I mentioned a couple of posts ago, clear limits can be established for the magnitudes of the coefficients in the case of a stable system. One can use these to constrain their values during the optimization, and I might investigate this as a way to insure stability even in the case that there is gain.
Yes I agree. Bob Orban's technique using polynomial curve fitting reaches the best fit but uses math that most folks are unfamiliar with and since he only published a cursory description one can only compare the shape of the error curve to his. This is probably mostly of academic interest.
Has anyone tried optimization methods described in this newly published article written by Martin Vicanek - Matched Second Order Digital Filters ?
In my case, I would need some universal real time working optimization method for biquad RIAA and non-RIAA filters. I'm calculating the coefficients this way:
Resulting coefficients are not accurate compared to ideal at above say 5kHz frequencies until the sampling frequency is above lets say 100kHz. Sure I can fix the accuracy issue by using higher internal sampling frequency but as up-/downsampling audio data may lead to some side-effects I would like to try avoid going that path.
In my case, I would need some universal real time working optimization method for biquad RIAA and non-RIAA filters. I'm calculating the coefficients this way:
Code:
double a0, a1, a2, b0, b1, b2;
double fs = 96000.0;
//timeconstants (case RIAA):
// frequency -> time conversion 1/(2*pi*fc) (= R*C)
//poles
double p1 = 3180e-6; // 1/(2*pi*50.05Hz)
double p2 = 75e-6; // 2212Hz
//zeros
double z1 = 318e-6; // 500.5Hz
double z2 = 0.0; // use 3.18e-6 for Neumann pole (50kHz)
double pole1= exp(-1.0/(fs*p1));
double pole2 = exp(-1.0/(fs*p2));
double zero1= exp(-1.0/(fs*z1));
double zero2 = exp(-1.0/(fs*z2));
a0 = 1.0; // 1.0
a1 = -pole1 - pole2; // -0.967774
a2 = pole1 * pole2; // 0.0
b0 = 1.0; // 1.0
b1 = -zero1 - zero2; // -1.867057
b2 = zero1 * zero2; // 0.867481Resulting coefficients are not accurate compared to ideal at above say 5kHz frequencies until the sampling frequency is above lets say 100kHz. Sure I can fix the accuracy issue by using higher internal sampling frequency but as up-/downsampling audio data may lead to some side-effects I would like to try avoid going that path.
OK, Martin Vicanek kindly looked this case and tweaked his above mentioned method suitable for custom RIAA type filter. By info he mailed he got quite good results (< ±0.3dB magnitude differnce between analog and digital filter in range 20Hz - 20kHz for 44.1kHz sampling frequency). As said, that's fair result but I'm looking solution which can give max ±0.1dB error. So, I did abandon this method (and few other methods which I found because of it looks like those are designed to work with those common filter types (peak, lp, hp, notch, etc.)) --> too much work for me to tweak math suitable for custom filters like what the RIAA and non-RIAA filters usually are.
OP:
I couldn't get the optimizer spreadsheet working on OpenOffice.Calc 4.1.2. Is it compatible with Calc and do I need some addons (optimizer.py or scsolver.uno or something else). Pressing "Solve" -button gives some error saying "Incorreft format for script URI".
OP:
I couldn't get the optimizer spreadsheet working on OpenOffice.Calc 4.1.2. Is it compatible with Calc and do I need some addons (optimizer.py or scsolver.uno or something else). Pressing "Solve" -button gives some error saying "Incorreft format for script URI".
Unfortunately it's not Calc compatible. It uses VBA, a Microsoft scripting language, as well as Microsoft's Solver. I am not aware of a non-linear equation solver for Calc.
If you use it in Excel with the SOLVER add-on installed it should work for you. If you need some help drop me a PM.
There are non-linear sovers for Calc as an extention.
Looks like VBA source can't be listed in Calc so would it be possible for you to share the VBA code so I could make it work with Calc?
I looked back on how I did the optimization of the biquads. Not much magic... I just set up the SOLVER optimization and then let er rip. Excel SOLVER'a "GRG Nonlinear" is a powerful and fast minimization algorithm...
Here is all the code behind the button:
That's it! There is no other code in the spreadsheet.
.
Here is all the code behind the button:
Code:
Sub OptimizeButton_Click()
' Macro adapted from Jon Peltier web page:
' [url=http://peltiertech.com/Excel/SolverVBA.html]Using Solver in Excel VBA[/url]
' Solver is called using Application.Run to avoid reference problems
'
'Copy Initial Guess to cells that are SOLVER input
Range("D15:D19").Copy Destination:=Range("D33")
'Set Up Solver
Application.Run "SolverOptions", , , 0.001, , , , , , , False, 0.000001, False, , , False, False
Application.Run "SolverOk", "$D$41", 2, , "$D$33:$D$37", , "GRG Nonlinear"
'Execute Solver Run
Application.Run "SolverSolve", True
End Sub.
- Status
- This old topic is closed. If you want to reopen this topic, contact a moderator using the "Report Post" button.