1bce563cdb1dde53bfccb7ab59aeddbf3f15d420
[scilab.git] / scilab / modules / signal_processing / macros / zpell.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - 1989 - F.Delebecque
3 // Copyright (C) INRIA - 1996 - C. Bunks
4 //
5 // Copyright (C) 2012 - 2016 - Scilab Enterprises
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13
14 function [ze,po,gain]=zpell(epsilon,A,omegac,omegar)
15     //[ze,po,gain]=zpell(epsilon,A,omegac,omegar)
16     //Poles and zeros of prototype lowpass elliptic filter
17     //gain is the gain of the filter
18     //  epsilon :Ripple of filter in pass band (0<epsilon<1)
19     //  A       :Attenuation of filter in stop band (A>1)
20     //  omegac  :Pass band cut-off frequency in rd/s
21     //  omegar  :Stop band cut-off frequency in rd/s
22     //  ze      :Resulting zeros of filter
23     //  po      :Resulting poles of filter
24     //  gain    :Resulting gain of filter
25     //
26     //!
27
28     m1=(epsilon*epsilon)/(A*A-1);
29     K1=delip(1,sqrt(m1));
30     K1t=imag(delip(1/sqrt(m1),sqrt(m1)));
31     m=(omegac/omegar)^2;
32     K=delip(1,sqrt(m));
33     Kt=imag(delip(1/sqrt(m),sqrt(m)));
34     n=(K1t*K)/(K1*Kt);
35     order=round(n);
36     u0=-(Kt/K1t)*delip(sqrt(1/(1+epsilon*epsilon)),sqrt(1-m1));
37     even=2*int(order/2);
38     if order<>even then,
39         vmin=2*K/n;
40     else,
41         vmin=K/n;
42     end,
43     v=vmin:(2*K/n):K;
44     un=ones(1:max(size(v)));
45     zlambda=-un*Kt+%i*v;
46     plambda= u0*un+%i*v;
47     ze=%i*imag(%i*omegac*%sn(-%i*zlambda,m));
48     ze=[ze,conj(ze)];
49     po=%i*omegac*%sn(-%i*plambda,m);
50     po=[po,conj(po)];
51     if order<>even then,
52         po=[po,%i*omegac*%sn(-%i*u0,m)];
53     end,
54     gain=abs(real(prod(po))/real(prod(ze)));
55     if order==even then,
56         gain=gain/sqrt(1+epsilon^2);
57     end;
58 endfunction