help-octave
[Top][All Lists]
Advanced

[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
#==================

reply via email to

[Prev in Thread] Current Thread [Next in Thread]