Hi,
I am using the butter command as follows : [b,a] = butter(12, 0.95); This works. If i increase the order to above 64, or want a cut off at the PI frequency (fs/2), Wc=1, with a very high order, then the filter returned seems to be in error. Example of failure : [b,a] = butter(12, 0.99); [b,a] = butter(128, 0.95); I assume i have passed the limits of the Octave capability. Currently using Fedora 26, Octave v4.2.1. I had a quick look at sourceforge - nothing about limits of the order or the cut off frequency. Any guidance gratefully received. Regards, Shadders. -- Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
________________________________ From: Shadders <[hidden email]> To: [hidden email] Sent: Friday, January 19, 2018 12:34 AM Subject: Butter Command - Maximum Order and Cut Off Frequency - Issues Hi, I am using the butter command as follows : [b,a] = butter(12, 0.95); This works. If i increase the order to above 64, or want a cut off at the PI frequency (fs/2), Wc=1, with a very high order, then the filter returned seems to be in error. Example of failure : [b,a] = butter(12, 0.99); [b,a] = butter(128, 0.95); I assume i have passed the limits of the Octave capability. Currently using Fedora 26, Octave v4.2.1. I had a quick look at sourceforge - nothing about limits of the order or the cut off frequency. Any guidance gratefully received. Regards, Shadders. -- Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html _______________________________________________ "Example of failure"- what are the exact error messages ? --Sergei _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
________________________________ From: Shadders <[hidden email]> To: [hidden email] 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. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
In reply to this post by Shadders
On 01/18/2018 06:56 PM, Shadders wrote:
> 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) > ================================= for d=0.96+[.010:.0001:.015]; [b,a] = butter(12,d ); y = filter(b, a, x); plot(n(500:3600), y(500:3600)),title(d); d,pause(.5); end I have no idea what is going on but it's kind of pretty :) _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
In reply to this post by Sergei Steshenko
Thanks Sergei,
Much appreciated, Regards, Shadders. -- Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
Free forum by Nabble | Edit this page |