[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
ripple in Butterworth filter created by 'butter' (was "RE: 'long double'
From: |
Sergei Steshenko |
Subject: |
ripple in Butterworth filter created by 'butter' (was "RE: 'long double' support ?") |
Date: |
Tue, 6 May 2008 13:52:53 -0700 (PDT) |
--- "Schirmacher, Rolf" <address@hidden> 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.
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
#==================