I haven't analyzed this filter, but generally speaking the poles
theoretically being in the unit circle isn't a guarantee of stability
when it comes time to implement the filter as a difference equation.
When doing so, it is possible that round-off errors nudge those poles,
in actuality, outside the unit circle if the theoretical poles are real
close to the unit circle, of which your example is true.
You may need to go to a higher order system, or you may need to
investigate doing the filter not as a single difference equation, but as
multiple stages of lower order filters that are each stable, i.e.,
"second order sections". That is why there is a set of routines with
the filter code for convert poles and zeros to second order
coefficients, etc.
Dan
On 11/12/2011 10:17 AM, Alberto Cuadrado wrote:
> I think the transfer function you provide is right (see help on "filter"
> function). In the vector "a" returned by butter.m, a(1) corresponds to 1
> and a(N+1) corresponds to a_N*z^-N, being "N" the filter order. However,
> the poles of the filter are roots(a), not 1./roots(a), because if you
> multiply by z^N both numerator and denominator of the transfer function,
> then a(1) corresponds to z^N and a(N+1) corresponds to a_N. If you want to
> use the "roots" function to obtain the z^-1 roots, you should flip
> (left-right) the vector "a" first. Answering your questions (I'm not sure,
> just guessing):
>
> (1) I don't know why I'm getting different abs(roots(a)). ¿Maybe different
> versions of Octave (3.2.4) and the signal package (1.0.11)?
>
> (2) That's a good question because a Butterworth filter shouldn't be
> unstable. I don't know the internals of butter.m very well but maybe it's
> related to the transform from s-plane to z-plane and/or numerical precision
> errors. But anyway you're right: if the filter is less narrow, the current
> implementation of the butter function returns a stable filter:
>
> [b,a]=butter(4,4e-5);abs(roots(a))
> ans =
>
> 0.99996
> 0.99996
> 0.99988
> 0.99988
>
>
> (3) It has to do with "poly" and "roots" functions not being exactly
> inverses of each other. I don't know if it's a bug or just unavoidable
> round-off errors.
>
>
> Alberto
>
> 2011/11/12 Jerry<lancebo...@qwest.net>
>
>> Hi Alberto,
>>
>> I hope you'll bear with me because I'm only an occasional user of Octave.
>>
>> I am experiencing frustration with understanding butter because I can't
>> find any documentation of what the results a and b actually are—there are
>> at least four different interpretations in use, including which is the
>> numerator and which is the denominator, and how the signs of the
>> coefficients are, which can be written in two different ways which derive
>> from how one originally wrote the difference equation. I am assuming that
>> the transfer function in question is this:
>>
>>
>>
>> using normal indexing notation, not Octave matrix-based notation.
>>
>> Today I updated my signals package after noting that there have been two
>> updates in the few days since I originally installed. I am now running
>> 1.1.1 but I don't see any changes for this problem
>>
>> On Nov 11, 2011, at 7:26 AM, Alberto Cuadrado wrote:
>>
>> Hello Jerry,
>>
>> The filter is actually unstable, so I think it's not a bug in filter, but
>> in butter:
>>
>> abs(roots(a))
>> ans =
>>
>> 1.000062437074018
>> 1.000062437074018
>> 0.999855481546637
>> 0.999855481546637
>>
>>
>> My result from running the code below (original posting) is this:
>>
>> abs(roots(a))
>> roots =
>>
>> 1.000131563207420
>> 0.999959711216214
>> 0.999959711216214
>> 0.999784860238433
>>
>> which is obviously different than your results.
>>
>> The actual roots that I get are
>>
>> octave-3.4.0:1> r = roots(a)
>> r =
>>
>> 1.000131563207420 + 0.000000000000000i
>> 0.999959694504841 + 0.000182815205375i
>> 0.999959694504841 - 0.000182815205375i
>> 0.999784860238433 + 0.000000000000000i
>>
>> which check out, giving 10^-15 type of answers.
>>
>> Since I have assumed (as stated above) that the polynomial is in z^-1,
>> then the actual poles in the z-plane are
>>
>> octave-3.4.0:3> rz = 1./r
>> rz =
>>
>> 0.999868454099180 + 0.000000000000000i
>> 1.000040273694318 - 0.000182829937069i
>> 1.000040273694318 + 0.000182829937069i
>> 1.000215186056644 + 0.000000000000000i
>>
>> and
>>
>> octave-3.4.0:6> abs(rz)
>> ans =
>>
>> 0.999868454099180
>> 1.000040290407038
>> 1.000040290407038
>> 1.000215186056644
>>
>> showing three unstable poles. No bueno.
>>
>> So: (1) Why did you get different abs(roots(a)), (2) why is butter
>> returning unstable results? and (3) why does butter return inconsistent
>> results when computing coefficients versus poles?
>>
>> It's possible (I suppose) that (2) is caused by my extremely narrow LPF
>> request. If so, I wonder if there is a way to change the computation so
>> that things don't get lost in the round-off.
>>
>> I wonder if the poles returned by butter are reliable and can be converted
>> into coefficients.
>>
>> Jerry
>>
>>
>>
>>
>> By the way, I think you have a mistake in the code. To obtain the
>> normalized frequency for the second input parameter of butter, it should be
>> f_cutoff/f_Nyquist, if I'm not wrong.
>>
>>
>> No, you're right.
>>
>>
>> Regards
>>
>> Alberto
>>
>>
>> 2011/11/11 Jerry<lancebo...@qwest.net>
>>
>>> Hi list,
>>>
>>> There is a bug in filter which is part of the signals add-on package
>>> which, when given stable poles, can generate an unstable impulse response.
>>> I probably have an extreme case (a very narrow passband low pass) but
>>> nonetheless, the above statement is true. I didn't actually factor the "a"
>>> polynomial from the first call to butter to see if it is consistent with
>>> the poles in the second call to butter. The poles from this fourth-order
>>> Butterworth are _very_ close to the unit circle but from the numerical
>>> results, there shouldn't be an unstable result--there are plenty of non-9's
>>> at the RHS of the numbers.
>>>
>>> Jerry
>>>
>>> signal *| 1.0.11 |
>>> OS X 10.6.8, using the pre-built Octave binary from the HPC site, 3.4.0.
>>>
>>> This code shows the problem.
>>>
>>> N = 500000;
>>> f_c = 1.0e9;
>>> T_counter = 1.0 / f_c;
>>> f_Nyquist = f_c / 2.0;
>>> f_cutoff = 2.0e4; % For more fun, make this 3.0e4.
>>> Filter_Order = 4;
>>>
>>> grid("on");
>>> format long;
>>>
>>> % Compute a Butterworth filter.
>>> [b, a] = butter(Filter_Order, f_cutoff / f_c);
>>> disp("b = "); disp(b);
>>> disp("a = "); disp(a);
>>>
>>> % Do zero-pole-gain form.
>>> [z, p, g] = butter(Filter_Order, f_cutoff / f_c);
>>> disp("z = "); disp(z);
>>> disp("p = "); disp(p);
>>> disp("g = "); disp(g);
>>> disp("Pole magnitudes = "); disp(abs(p));
>>>
>>> % Make time vector, impulse; filter the impulse.
>>> t = [0 : T_counter : (N - 1) * T_counter];
>>> Impulse = [1, zeros(1, N - 1)];
>>> h = filter(b, a, Impulse);
>>>
>>> % Plot the impulse response.
>>> plot(t, h); grid("on");
>>> pause;
>>>
>>> ------------------------------------------------------------------------------
>>> RSA(R) Conference 2012
>>> Save $700 by Nov 18
>>> Register now
>>> http://p.sf.net/sfu/rsa-sfdev2dev1
>>> _______________________________________________
>>> Octave-dev mailing list
>>> Octave-dev@lists.sourceforge.net
>>> https://lists.sourceforge.net/lists/listinfo/octave-dev
>>>
>>
>>
>>
>>
>
>
>
> ------------------------------------------------------------------------------
> RSA(R) Conference 2012
> Save $700 by Nov 18
> Register now
> http://p.sf.net/sfu/rsa-sfdev2dev1
>
>
>
> _______________________________________________
> Octave-dev mailing list
> Octave-dev@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/octave-dev
--
Dan Sebald
email: daniel(DOT)sebald(AT)ieee(DOT)org
URL: http://www(DOT)dansebald(DOT)com
------------------------------------------------------------------------------
RSA(R) Conference 2012
Save $700 by Nov 18
Register now
http://p.sf.net/sfu/rsa-sfdev2dev1
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
