[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Butter Command - Maximum Order and Cut Off Frequency - Issues
From: |
Sergei Steshenko |
Subject: |
Re: Butter Command - Maximum Order and Cut Off Frequency - Issues |
Date: |
Fri, 19 Jan 2018 00:21:55 +0000 (UTC) |
________________________________
From: Shadders <address@hidden>
To: address@hidden
Sent: Friday, January 19, 2018 1:57 AM
Subject: Re: Butter Command - Maximum Order and Cut Off Frequency - Issues
Hi,
No error messages, but the plots have bad values.
The following works :
[b,a] = butter(12, 0.95) - plots are ok (multiple sine waves, different
frequencies, distinct)
The following plots a single exponentially increasing sine wave which
increases to values from 3e305 to -3e305, and terminates plotting at approx
1/8 of the expected values for the x axis :
[b,a] = butter(12, 0.98)
The output of the filter has from the 1/8 distance, NaN in the remainder of
the vector.
Code that replicates failure - simple
=================================
nc = 5; % Number of cycles of the waveform
fs = 192000;
ts = 1/fs;
n = 0:192000; % Variable used for indexing
x = zeros(1, length(n)); % Input signal initialisation
[b,a] = butter(12, 0.98); % Create the butterworth response
k1 = fs/1000;
k1n=0:(nc*k1);
x1k = sin(2*pi*k1n/k1);
x(1, 1000:1000+length(x1k)-1) = x1k;
y = filter(b, a, x); % Use the filter function
figure(1);
plot(n, y)
=================================
Regards,
Shadders.
--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
_______________________________________________
"the plots have bad values" - this is not a limitation of Octave, it's a
limitation of floating point arithmetic. I had the same problem with SciLab and
with Octave.
The workaround is not to use the direct polynomial form (B(z)/A(z)), but rather
to use roots of B (zeros) and roots of A (poles) and to calculate the value
through roots - it decreases rounding errors.
Read the function documentation to see how the needed roots can be obtained and
you can use this (my "production" function) to calculate values:
###########################################################################################################
function freq_response = eval_filter_through_zpg(zeros_vector, poles_vector,
gain_scalar, zfrequency_range)
numerator = ones(1, length(zfrequency_range));
for zero_number = 1:length(zeros_vector)
numerator = numerator .* (zfrequency_range - zeros_vector(zero_number));
endfor
denominator = ones(1, length(zfrequency_range));
for pole_number = 1:length(poles_vector)
denominator = denominator .* (zfrequency_range -
poles_vector(pole_number));
endfor
freq_response = gain_scalar * numerator ./ denominator;
endfunction
#==========
zfrequency_range can be understood from
two_pi = 2 * pi;
sample_rate = 44100; # change as needed
number_of_points_in_fft = 2 ^ 16; # change as needed
half_number_of_points_in_fft = number_of_points_in_fft / 2;
relative_all_freqs = (0:half_number_of_points_in_fft) / number_of_points_in_fft;
relative_omegas = two_pi * relative_all_freqs;
relative_iomegas = i * relative_omegas;
zfrequency_range = exp(relative_iomegas);
...
Even with the proposed workaround you can still be hit with floating point
rounding issues, but in more extreme cases (i.e. higher order filter, closer to
Nyquist cutoff frequency).
This all should be in the documentation for 'butter' function, but I'm not an
Octave developer.
--Sergei.