ripple in Butterworth filter created by 'butter' (was "RE: 'long double' support ?")

[Sergei Steshenko-2]

--- "Schirmacher, Rolf" <Rolf.Schirmacher@...> wrote:

>
> Perhaps you might give more details on the package / function you encouter
> the issue with and/or your actual problem as the root cause and/or best
> solution / workaround might be different to a "brute force" approach of
> increasing precision?
>
> Rolf
>

Hello All,

the problem I had was with Butterworth filter created by 'butter'.

The problem was/is ripple in the passband - there can be no ripple anywhere because
of the definition of Butterworth filter.

The problem is not specific to 'octave' - a couple of years ago I found the same problem in
'scilab'.

Attached is a script which demonstrates the problem.

The script also demonstrates how to avoid the problem using the

"
[z,p,g] = butter(...)
   return filter as zero-pole-gain rather than coefficients of the
   numerator and denominator polynomials.
"

form of 'butter' - I think this should be THE recommended way to use the function.

With

zorder = 4

the ripple is about 0.006; with

zorder = 5

or higher the filter becomes completely incoherent - just try to run the script.


I suggest to copy-paste the script into documentation - this will hopefully help others
to avoid the trap. And it's a usage example which isn't in the documentation yet

Thanks,
  Sergei.

Applications From Scratch: http://appsfromscratch.berlios.de/


      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile.  Try it now.  http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ
# written by Sergei Steshenko.

max_frequency = 10000;


frequency_range = 1:max_frequency;
zfrequency_range = exp(i * pi * frequency_range / max_frequency);


lower_zcutoff = 0.01;
higher_zcutoff = 0.02;

zorder = 4; # try to change this to 5 or more to see drastic influence of numeric accuracy issues
            # even at 4 ripple is visible in passband

figure_number = 1;

#-------------------
# bad accuracy BEGIN

[bad_accuracy_nominator, bad_accuracy_denominator] = butter(zorder, [lower_zcutoff, higher_zcutoff]);


bad_accuracy_frequency_response = polyval(bad_accuracy_nominator, zfrequency_range) ./ polyval(bad_accuracy_denominator, zfrequency_range);
figure(figure_number++);
semilogx(frequency_range, abs(bad_accuracy_frequency_response));
title("Bad accuracy frequency response - ripple in passband");

figure(figure_number++);
semilogx(frequency_range, unwrap(arg(bad_accuracy_frequency_response), pi), frequency_range, arg(bad_accuracy_frequency_response));
title("Bad accuracy phase response");

# bad accuracy END
#=================


#--------------------
# good accuracy BEGIN

[good_accuracy_zeros, good_accuracy_poles, good_accuracy_gain] = butter(zorder, [lower_zcutoff, higher_zcutoff]);

znom = ones(1, length(frequency_range));
for zero_number = 1:length(good_accuracy_zeros)
  znom = znom .* (good_accuracy_zeros(zero_number) - zfrequency_range);
endfor

zdenom = ones(1, length(frequency_range));
for pole_number = 1:length(good_accuracy_poles)
  zdenom = zdenom .* (good_accuracy_poles(pole_number) - zfrequency_range);
endfor

good_accuracy_frequncy_response = good_accuracy_gain * znom ./ zdenom;

figure(figure_number++);
semilogx(frequency_range, abs(good_accuracy_frequncy_response));
title("Good accuracy frequency response - NO ripple in passband");

figure(figure_number++);
semilogx(frequency_range, unwrap(arg(good_accuracy_frequncy_response), pi), frequency_range, arg(good_accuracy_frequncy_response));
title("Good accuracy phase response");

# good accuracy END
#==================

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave