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
