[Scilab-users] Unexpected result using horner

[Samuel Gougeon]

Le 25/07/2022 à 02:44, Federico Miyara a écrit :
>
> Samuel,
>
>> Sorry, but nothing is clear to me about your statements:
>>
>> 1) what shows you that iir() is correct while the analpf + horner way 
>> is not?
>
> There are several reasons:
>
> a) The result of using iir() is consistent with the expected frequency 
> response, since it is equal to the analog filter response except at 
> high frequency, an expected artifact for IIR filters. The other 
> solution is completely different.
>
> b) The degree of the denominator should be 6 but when using horner it 
> reduces to 3.
>
> My own informal and not very deep analysis suggests that as all the 
> poles are very close to unity, may be horner() performs some 
> simplification and simplifies things incorrectly, which might 
> dramatically change the frequency response.
>
>> 2) With the analpf + horner method, assuming that it is not correct, 
>> what shows you that horner is not correct, while analpf is correct, 
>> instead of the opposite or both incorrect?
>
> analpf() is a very simple algorithm, at least for Butterworth, since 
> explicit formulas for the poles exist and in all the cases I have 
> tested the result is the expected one, particularly in this case.
>
>> 3) do you have a reference about the equivalence?
>
> The substitution of the bilinear transformation is the usual method to 
> get an IIR digital filter from an analog prototype. They aren't 
> completely equivalent, but very similar up to about half Nyquist 
> frequency. Any book on digital signal processing includes that 
> transformation. Also https://en.wikipedia.org/wiki/Bilinear_transform
>
>> 4) have you tried after simp_mode(%f)?
>
> No, but this might be the answer considering my reply to 2). I'll try 
> it later.


To get equivalent transfer functions, you need
a) indeed to turn simp_mode(%f),
b) to run
hBPz1= iir(3, "bp", "butt", fo/Fs*[2^(-1/6), 2^(1/6)], [0 0])
instead of
hBPz1= iir(3, "bp", "butt", fo/Fs*[2^(1/6),2^(-1/6)], [0 0])

Then the curves of both transfer functions are superimposed.

clf,  plot(linspace(0.5,1.5,100),  [hBPsz,  hBPz1]) When working with polynomials or rationals, one must keep in 
mind that every time that their roots or/and poles have close values, 
then results become very sensitive to numerical round-off errors or 
uncertainties. The roots of a polynomial of degree N can't be computed 
from its coefficients with a numerical relative accuracy better than 
%eps^(1/N). If this uncertainty becomes bigger than the distance between 
roots values, then things become really bad.
Samuel

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20220725/90881584/attachment.htm>