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
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.
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
