Re: AW: Bessel-filter problems

[Sergei Steshenko]

--- On Tue, 4/26/11, Bernd Heinze <address@hidden> wrote:

> From: Bernd Heinze <address@hidden>
> Subject: AW: Bessel-filter problems
> To: address@hidden
> Date: Tuesday, April 26, 2011, 1:10 PM
> Hello Christian,
> 
> thanks for your reply. Unfortunately I'm not an filter
> expert, I played around with Diadem and decided that a
> Bessel filter with said parameters gives me the best
> results. I assumed if it works in Diadem it will also work
> in Octave.
> 
> 
> > But are you sure that you really want to construct
> such a filter? Your ratio
> > between corner frequencies and the Nyquist frequency
> (2.5 GHz ???) is quite
> > high. Throw in a filter order of 10 and you will get a
> filter with
> > coefficients that are numerically ill-behaved.
> 
> The data is taken from an oscilloscope which has a sampling
> rate of 5 GHz. As my signal contains many high frequencies
> this sampling rate is necessary. (Although it might be
> possible to reduce it by a factor of 2 or 4)
> 
> > Constructing a butterworth filter from your data shows
> what I mean: The
> > resulting coefficients range from 1 to 1e5, i.e. they
> differ in scale by
> > about 17 bits.
> 
> I assume this means that by using the filter I might get
> larger "errors" or imprecissions in my signal?
> 
> Regards,
> Bernd

One has to be quite careful dealing with numeric accuracy. Please find
below an Email from 2008 of mine with a practical demo/recipe attached
(butterworth_numeric_accuracy.m).

Regards,
  Sergei.

The Email from 2008:

--- On Tue, 5/6/08, Sergei Steshenko <address@hidden> wrote:

> From: Sergei Steshenko <address@hidden>
> Subject: ripple in Butterworth filter created by 'butter' (was "RE: 'long 
> double' support ?")
> To: "Schirmacher, Rolf" <address@hidden>, address@hidden, address@hidden, 
> address@hidden
> Cc: address@hidden
> Date: Tuesday, May 6, 2008, 1:52 PM
> 
> --- "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
> -----Inline Attachment Follows-----
> 
> 
> -----Inline Attachment Follows-----
> 
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www.cae.wisc.edu/mailman/listinfo/help-octave
>

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