1a2
> ## Copyright (C) 2015 Charles Praplan <charles.praplan@alumni.epfl.ch>
19c20,23
< ## @deftypefnx {Function File} {[@var{n}, @var{wc}] =} buttord (@dots{})
---
> ## @deftypefnx {Function File} {@var{n} =} buttord ([@var{wp1}, @var{wp2}], [@var{ws1}, @var{ws2}], @var{rp}, @var{rs}, 's')
> ## @deftypefnx {Function File} {[@var{n}, @var{wc_p}] =} buttord (@dots{})
> ## @deftypefnx {Function File} {[@var{n}, @var{wc_p}, @var{wc_s}] =} buttord (@dots{})
> ##
25,26c29,35
< ## attenuation in the stop band, also in decibels.  The output arguments @var{n}
< ## and @var{wc} can be given as inputs to @code{butter}.
---
> ## attenuation in the stop band, also in decibels.
> ##
> ## The output arguments @var{n} and @var{wc_p} (or @var{n} and @var{wc_n}) can be
> ## given as inputs to @code{butter}.
> ## Using @var{wc_p} makes the filter characteristic touch at least one pass band
> ## corner and using @var{wc_s} makes the characteristic touch at least one
> ## stop band corner.
40c49,52
< ## Theory: |H(W)|^2 = 1/[1+(W/Wc)^(2N)] = 10^(-R/10)
---
> ## For Laplace space filters, an additional input parameter ('s') must be entered.
> ## In that case W is in rad/s.
> ##
> ## Theory: For Low pass filters, |H(W)|^2 = 1/[1+(W/Wc)^(2N)] = 10^(-R/10).
42,46c54,61
< ## Ws,Rs and Wp,Rp.  For high pass filters, subtracting the band edges
< ## from Fs/2, performing the test, and swapping the resulting Wc back
< ## works beautifully.  For bandpass and bandstop filters this process
< ## significantly overdesigns.  Artificially dividing N by 2 in this case
< ## helps a lot, but it still overdesigns.
---
> ## Ws,Rs and Wp,Rp. Rounding N to the next greater integer, one can recalculate
> ## the allowable range for Wc (filter caracteristic touching the pass band edge or
> ## the stop band edge).
> ##
> ## For other types of filter, before making the above calculation, the requirements
> ## must be transformed to LP requirements. After calculation, Wc must be transformed
> ## back to original filter type.
> ##
51c66
< function [n, Wc] = buttord(Wp, Ws, Rp, Rs)
---
> function [n, Wc_p, Wc_s] = buttord(Wp, Ws, Rp, Rs, opt)
53c68
<   if (nargin != 4)
---
>   if (nargin != 4 && nargin != 5)
54a70,73
>   endif;
> 
>   if (nargin == 5 && strcmp(opt,'s'))
>     s_domain = 1;
55a75,78
>     s_domain = 0;
>   endif;
> 
>   if s_domain == 0
56a80,81
>   else
>     validate_filter_bandss ("buttord", Wp, Ws);
59,60c84,92
<   if (numel (Wp) == 2)
<     warning ("buttord: seems to overdesign bandpass and bandreject filters");
---
>   if s_domain == 1
>     Wpw = Wp; # No prewarp in case of analog filter
>     Wsw = Ws; # No prewarp in case of analog filter
>   else
>     ## sampling frequency of 2 Hz
>     T = 2;
> 
>     Wpw = (2/T)*tan(pi.*Wp./T); # prewarp
>     Wsw = (2/T)*tan(pi.*Ws./T); # prewarp
63c95,135
<   T = 2;
---
>   ## pass/stop band to low pass filter transform:
>   if (length(Wpw)==2 && length(Wsw)==2)
> 
>     if Wpw(1) > Wsw(1) % Band pass filter
> 
>       if Wpw(1)*Wpw(2) < Wsw(1)*Wsw(2)  %Modify band edges if not symmetrical
>         Wsw(2) = Wpw(1)*Wpw(2)/Wsw(1); % smaller stopband
>         %Wpw(2) = Wsw(1)*Wsw(2)/Wpw(1); % larger passband (less efficient for BP)
>       else
>         Wsw(1) = Wpw(1)*Wpw(2)/Wsw(2); % smaller stopband
>         %Wpw(1) = Wsw(1)*Wsw(2)/Wpw(2); % larger stopband (less efficient for BP)
>       end;
> 
>       w02 = Wpw(1)*Wpw(2);
>       wp=Wpw(2)-Wpw(1);
>       ws=Wsw(2)-Wsw(1);
> 
>     else               % Notch filter
> 
>       if Wpw(1)*Wpw(2) > Wsw(1)*Wsw(2)  %Modify band edges if not symmetrical
>         %Wsw(2) = Wpw(1)*Wpw(2)/Wsw(1); % larger stopband
>         Wpw(2) = Wsw(1)*Wsw(2)/Wpw(1); % smaller passband (less efficient for notch)
>       else
>         %Wsw(1) = Wpw(1)*Wpw(2)/Wsw(2); % larger stopband
>         Wpw(1) = Wsw(1)*Wsw(2)/Wpw(2); % smaller passband (less efficient for notch)
>       end;
> 
>       w02 = Wpw(1)*Wpw(2);
>       wp=w02/(Wpw(2)-Wpw(1));
>       ws=w02/(Wsw(2)-Wsw(1));
>     end
>     ws=ws/wp;
>     wp=1;
> 
>   elseif (Wpw > Wsw)   % High pass filter
>     wp = Wsw;
>     ws = Wpw;
>   else                 % Low pass filter
>     wp = Wpw;
>     ws = Wsw;
>   end
65,72d136
<   ## if high pass, reverse the sense of the test
<   stop = find(Wp > Ws);
<   Wp(stop) = 1-Wp(stop); # stop will be at most length 1, so no need to
<   Ws(stop) = 1-Ws(stop); # subtract from ones(1,length(stop))
< 
<   ## warp the target frequencies according to the bilinear transform
<   Ws = (2/T)*tan(pi*Ws./T);
<   Wp = (2/T)*tan(pi*Wp./T);
75,76d138
<   ## the factor 1/length(Wp) is an artificial correction for the
<   ## band pass/stop case, which otherwise significantly overdesigns.
79c141
<   n = ceil(max(0.5*(qs - qp)./log(Ws./Wp))/length(Wp));
---
>   n = ceil(max(0.5*(qs - qp)./log(ws./wp)));
82c144
<   Wc = exp(log(Wp) - qp/2/n);
---
>   if (length(Wpw)==2 && length(Wsw)==2)   % Band pass or notch filter
84,85c146,195
<   ## unwarp the returned frequency
<   Wc = atan(T/2*Wc)*T/pi;
---
>     if Wpw(1) > Wsw(1)                    % Band pass filter
>       w_prime_p = exp(log(Wpw) - qp/2/n); %   same formula as for LP
>       w_prime_s = exp(log(Wsw) - qs/2/n); %           "
>     else                                  % Notch filter
>       w_prime_p = exp(log(Wpw) + qp/2/n); %   same formula as for HP
>       w_prime_s = exp(log(Wsw) + qs/2/n); %           "
>     end
> 
>     % Applying LP to BP (respectively HP to notch) transformation to -3dB
>     % angular frequency :
>     %   s_prime/wc = Q(s/w0+w0/s)  or  w_prime/wc = Q(w/w0-w0/w)
>     % Here we need to inverse above equation:
>     %   w = abs(w_prime+-sqrt(w_prime^2+4*Q^2))/(2*Q/w0);
> 
>     % -3dB cutoff freq to match pass band
>     w0=sqrt(prod(Wpw));
>     Q=w0/diff(Wpw);      % BW at -Rp dB not at -3dB
>     wc=Wpw;
>     W_prime = w_prime_p(1)/wc(1); % same result with w_prime(2)/wc(2)
>     wa =abs(W_prime+sqrt(W_prime^2+4*Q^2))/(2*Q/w0);
>     wb =abs(W_prime-sqrt(W_prime^2+4*Q^2))/(2*Q/w0);
>     Wcw_p = [wb wa];
> 
>     % -3dB cutoff freq to match stop band
>     w0=sqrt(prod(Wsw));
>     Q=w0/diff(Wsw);      % BW at -Rs dB not at -3dB
>     wc=Wsw;
>     W_prime = w_prime_s(1)/wc(1); % same result with w_prime(2)/wc(2)
>     wa =abs(W_prime+sqrt(W_prime^2+4*Q^2))/(2*Q/w0);
>     wb =abs(W_prime-sqrt(W_prime^2+4*Q^2))/(2*Q/w0);
>     Wcw_s = [wb wa];
> 
>   elseif (Wpw > Wsw)   % High pass filter
>     Wcw_p = exp(log(Wpw) + qp/2/n); % -3dB cutoff freq to match pass band
>     Wcw_s = exp(log(Wsw) + qs/2/n); % -3dB cutoff freq to match stop band
>   else                 % Low pass filter
>     Wcw_p = exp(log(Wpw) - qp/2/n); % -3dB cutoff freq to match pass band
>     Wcw_s = exp(log(Wsw) - qs/2/n); % -3dB cutoff freq to match stop band
>   end
> 
> 
>   if s_domain == 1 # No prewarp in case of analog filter
>     Wc_p = Wcw_p;
>     Wc_s = Wcw_s;
>   else
>     ## sampling frequency of 2 Hz
>     T = 2;
>     Wc_p = atan(Wcw_p/(2/T)) * (T/pi); # inverse prewarp after evt. symmetrisation
>     Wc_s = atan(Wcw_s/(2/T)) * (T/pi); # inverse prewarp after evt. symmetrisation
>   end
87,88d196
<   ## if high pass, reverse the sense of the test
<   Wc(stop) = 1-Wc(stop);
