2a3
> ## Copyright (C) 2015 Charles Praplan <charles.praplan@alumni.epfl.ch>
19a21
> ## @deftypefnx {Function File} {@var{n} =} cheb1ord ([@var{wp1}, @var{wp2}], [@var{ws1}, @var{ws2}], @var{rp}, @var{rs}, 's')
20a23
> ## @deftypefnx {Function File} {[@var{n}, @var{wc_p}, @var{wc_s}] =} cheb1ord (@dots{})
26,27c29,35
< ## is the minimum attenuation in the stop band, also in decibels.  The output
< ## arguments @var{n} and @var{wc} can be given as inputs to @code{cheby1}.
---
> ## is the minimum attenuation in the stop band, also in decibels.
> ##
> ## The output arguments @var{n} and @var{wc_p} (or @var{n} and @var{wc_s}) can be
> ## given as inputs to @code{cheby1}.
> ## 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.
40a49,51
> ## For Laplace space filters, an additional input parameter ('s') must be entered.
> ## In that case W is in rad/s.
> ##
44c55,59
< function [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
---
> function [n, Wc_p, Wc_s] = cheb1ord(Wp, Ws, Rp, Rs, opt)
> 
>   if (nargin != 4 && nargin != 5)
>     print_usage ();
>   endif;
46,47c61,62
<   if nargin != 4
<     print_usage;
---
>   if (nargin == 5 && strcmp(opt,'s'))
>     s_domain = 1;
48a64,67
>     s_domain = 0;
>   endif;
> 
>   if s_domain == 0
49a69,70
>   else
>     validate_filter_bandss ("cheb1ord", Wp, Ws);
52c73,78
<   T = 2;
---
>   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;
54,55c80,82
<   ## returned frequency is the same as the input frequency
<   Wc = Wp;
---
>     Wpw = (2/T)*tan(pi.*Wp./T); # prewarp
>     Wsw = (2/T)*tan(pi.*Ws./T); # prewarp
>   endif
57,78c84,126
<   ## warp the target frequencies according to the bilinear transform
<   Ws = (2/T)*tan(pi*Ws./T);
<   Wp = (2/T)*tan(pi*Wp./T);
< 
<   if (Wp(1) < Ws(1))
<     ## low pass
<     if (length(Wp) == 1)
<       Wa = Ws/Wp;
<     else
<       ## FIXME: Implement band reject filter type
<       error ("cheb1ord: band reject is not yet implemented");
<     endif;
<   else
<     ## if high pass, reverse the sense of the test
<     if (length(Wp) == 1)
<       Wa = Wp/Ws;
<     else
<       ## band pass
<       Wa=(Ws.^2 - Wp(1)*Wp(2))./(Ws*(Wp(1)-Wp(2)));
<     endif;
<   endif;
<   Wa = min(abs(Wa));
---
>   ## 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
> 
>   Wa = ws./wp;
83a132,192
> 
> 
>   ## compute stopband frequency limits to make the the filter characteristic
>   ## touch either at least one stop band corner or one pass band corner.
> 
>   epsilon=1/sqrt(10^(.1*abs(Rs)) - 1);
>   k=cosh(1/n*acosh(sqrt(1/(10^(.1*abs(Rp)) - 1))/epsilon)); % k=fstop/fpass
> 
>   if (length(Wpw)==2 && length(Wsw)==2)   % Band pass or notch filter
> 
> 
>     if Wpw(1) > Wsw(1)                    % Band pass filter
>       w_prime_p = Wpw;                    %   same formula as for LP
>       w_prime_s = Wsw/k;                  %           "
>     else                                  % Notch filter
>       w_prime_p = Wpw;                    %   same formula as for HP
>       w_prime_s = k*Wsw;                  %           "
>     end
> 
>     % Applying LP to BP (respectively HP to notch) transformation to
>     % angular frequency to be returned:
>     %   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);
> 
>     % freq to be returned 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];
> 
>     % freq to be returned 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 = Wpw;         % to match pass band
>     Wcw_s = Wsw*k;       % to match stop band
>   else                 % Low pass filter
>     Wcw_p = Wpw;         % to match pass band
>     Wcw_s = Wsw/k;       % 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
