0
$\begingroup$

Can anyone explain to me how to use this equation to find the Order? Let's say

$$\Omega_p/\Omega_s =0.509$$ and $$e/\delta = 0.00349$$

How can I find the Order? I'm so lost i don't know how to use this equation.

If there is another equation for determining the order to elliptic filters please provide it.

$$N = \frac{K(\Omega_p/\Omega_s)K\left(\sqrt{1-(\epsilon^2/\delta^2)}\right)}{K(\epsilon/\delta)K\left(\sqrt{1-(\Omega_p/\Omega_s)^2}\right)}$$

$$K(x) = \int_{0}^{\pi/2} \frac{d\theta}{\sqrt{1-x^2\sin^2\theta}}$$

$\endgroup$
1
  • $\begingroup$ You'll need a numerical solver to calculate the complete elliptic integral K(x). Octave (Matlab?) has it as ellipke(), wxMaxima as elliptic_kc(), for example. You'll need to check with the software you're using, or make your own (Gaussian quadrature works very well). $\endgroup$ Dec 18, 2022 at 17:23

2 Answers 2

3
$\begingroup$

If we apply your formula, we get:

# Python
from math import  ceil
from numpy import sqrt
from scipy.special import ellipk
def K(x):
    return ellipk(x*x)

omega_ratio = 0.509
e_delta_ratio = 0.00349

N = K(omega_ratio)*K(sqrt(1-e_delta_ratio*e_delta_ratio))/(K(e_delta_ratio)*K(sqrt(1-omega_ratio*omega_ratio)))
print(N)   # 3.5420113113759575
print(ceil(N))  # 4
% matlab
k1 = 0.00349;
k = 0.509;
capk = ellipke([k^2 1-k^2]);
capk1 = ellipke([(k1^2) 1-(k1^2)]);
tmp=(capk(1)*capk1(2)/(capk(2)*capk1(1)))  % 3.5420
order=ceil(tmp)   % 4

I don't know where the formula you give comes from, nor do I know whether that formula is correct.

I think the algorithm in matlab ellipord function may be more reliable.

Why don't you use matlab ellipord function?

By the way, it looks like

function [order,wn] = ellipord(wp,ws,rp,rs,opt)
%ELLIPORD Elliptic filter order selection.
%   [N, Wp] = ELLIPORD(Wp, Ws, Rp, Rs) returns the order N of the lowest
%   order digital elliptic filter which has a passband ripple of no more
%   than Rp dB and a stopband attenuation of at least Rs dB. Wp and Ws are
%   the passband and stopband edge frequencies, normalized from 0 to 1
%   (where 1 corresponds to pi radians/sample). For example,
%       Lowpass:    Wp = .1,      Ws = .2
%       Highpass:   Wp = .2,      Ws = .1
%       Bandpass:   Wp = [.2 .7], Ws = [.1 .8]
%       Bandstop:   Wp = [.1 .8], Ws = [.2 .7]
%   ELLIPORD also returns Wp, the elliptic natural frequency to use with
%   ELLIP to achieve the specifications.
%
%   [N, Wp] = ELLIPORD(Wp, Ws, Rp, Rs, 's') does the computation for an
%   analog filter, in which case Wp and Ws are in radians/second.
%
%   NOTE: If Rs is much much greater than Rp, or Wp and Ws are very close,
%   the estimated order can be infinite due to limitations of numerical
%   precision.
%
%   % Example 1:
%   %   For 1000 Hz data, design a lowpass filter with less than 3 dB of
%   %   ripple in the passband defined from 0 to 40 Hz and at least 60 dB
%   %   of ripple in the stopband defined from 150 Hz to the Nyquist
%   %   frequency (500 Hz):
%
%   Wp = 40/500; Ws = 150/500;
%   Rp = 3; Rs = 60;
%   [n,Wp] = ellipord(Wp,Ws,Rp,Rs)      % Gives mimimum order of filter
%   [b,a] = ellip(n,Rp,Rs,Wp);          % Elliptic filter design
%   freqz(b,a,512,1000);                % Plots the frequency response
%   title('n=4 Elliptic Lowpass Filter')
%
%   % Example 2:
%   %   Design a bandpass filter with a passband of 60 Hz to 200 Hz, with
%   %   less than 3 dB of ripple in the passband, and 40 dB attenuation in
%   %   the stopbands that are 50 Hz wide on both sides of the passband.
%
%   Wp = [60 200]/500; Ws = [50 250]/500;
%   Rp = 3; Rs = 40;
%   [n,Wp] = ellipord(Wp,Ws,Rp,Rs)      % Gives mimimum order of filter
%   [b,a] = ellip(n,Rp,Rs,Wp);          % Elliptic filter design
%   freqz(b,a,512,1000)                 % Plots the frequency response
%
%   See also ELLIP, BUTTORD, CHEB1ORD, CHEB2ORD, DESIGNFILT.

%   Author(s): L. Shure, 6-9-88
%              T. Krauss, 11-18-92, updated
%   Copyright 1988-2018 The MathWorks, Inc.

%   Reference(s):
%       [1] Rabiner and Gold, p 241.

narginchk(4,5);  
nargoutchk(0,2); 

if nargin > 4
    opt = convertStringsToChars(opt);
end

% Cast to enforce precision rules
wp = signal.internal.sigcasttofloat(wp,'double','ellipord','Wp','allownumeric');
ws = signal.internal.sigcasttofloat(ws,'double','ellipord','Ws','allownumeric');
rp = signal.internal.sigcasttofloat(rp,'double','ellipord','Rp','allownumeric');
rs = signal.internal.sigcasttofloat(rs,'double','ellipord','Rs','allownumeric');

if nargin == 4
    opt = 'z';
elseif nargin == 5
    if ~strcmp(opt,'z') && ~strcmp(opt,'s')
        error(message('signal:ellipord:InvalidParam'));
    end
end

[msg,msgobj]=freqchk(wp,ws,opt);
if ~isempty(msg), error(msgobj); end

ftype = 2*(length(wp) - 1);
if wp(1) < ws(1)
    ftype = ftype + 1;  % low (1) or reject (3)
else
    ftype = ftype + 2;  % high (2) or pass (4)
end

% first, prewarp frequencies from digital (unit circle) to analog (imag. axis):
if strcmp(opt,'z')  % digital
    WP=tan(pi*wp/2);
    WS=tan(pi*ws/2);
else  % don't have to if analog already
    WP=wp;
    WS=ws;
end

% next, transform to low pass prototype with passband edge of 1 and stopband
% edges determined by the following: (see Rabiner and Gold, p.258)
if ftype == 1   % low
    WA=WS/WP;
    order = findelliporder(WA,rp,rs);
elseif ftype == 2   % high
    WA=WP/WS;
    order = findelliporder(WA,rp,rs);
elseif ftype == 3   % stop
    if strcmp(opt,'s')
        % For analog bandstop, convert back to digital
        wp = 2*atan(wp)/pi;
        ws = 2*atan(ws)/pi;
    end
    Fpass1 = wp(1);
    Fpass2 = wp(2);
    Fstop1 = ws(1);
    Fstop2 = ws(2);
    c = sin(pi*(Fpass1+Fpass2))/(sin(pi*Fpass1)+sin(pi*Fpass2));
    wpa = abs(sin(pi*Fpass2)/(cos(pi*Fpass2)-c));
    ws1 = sin(pi*Fstop1)/(cos(pi*Fstop1)-c);
    ws2 = sin(pi*Fstop2)/(cos(pi*Fstop2)-c);
    wsa = min(abs([ws1,ws2]));
    order = ellipord(wpa,wsa,rp,rs,'s');
elseif ftype == 4   % pass
    WA=(WS.^2 - WP(1)*WP(2))./(WS*(WP(1)-WP(2)));
    order = findelliporder(WA,rp,rs);
end


% natural frequencies are simply the passband edges (WP).
% finally, transform frequencies from analog to digital if necessary:
if strcmp(opt,'z')  % digital
    wn = wp;
else
    wn = WP;
end

%--------------------------------------------------------------------------
function order = findelliporder(WA,rp,rs)
% find the minimum order elliptic filter to meet the more demanding spec:
WA = min(abs(WA));
epsilon = sqrt(10^(0.1*rp)-1);
k1 = epsilon/sqrt(10^(0.1*rs)-1);
k = 1/WA;
capk = ellipke([k^2 1-k^2]);
capk1 = ellipke([(k1^2) 1-(k1^2)]);
order = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1)));

% if both warnings are in effect, only print the first one
if (1-k1^2) == 1
    warning(message('signal:ellipord:MustBeFiniteAttn', 'ellipord'))
elseif k^2 == 1
    warning(message('signal:ellipord:MustBeFiniteBandEdges', 'ellipord'))
end


For science calculation, you can find many APIs in different programming languages.

# python
from  scipy.special import ellipk
ellipk(0.509)

# matlab
ellipke(0.509)

# mathematica
EllipticK[0.509]


# All of them output the same value:  1.8617756243745627
# All of them are  Complete elliptic integral of the first kind 

They all imply the same integral: $$ \int_0^{\frac{\pi }{2}} \frac{1}{\sqrt{1-0.509 \sin ^2(\theta)}} \, \mathrm{d} \theta \approx 1.86178 $$

% however, for matlab there are some weird syntax
[a,b]=ellipke(0.509)     %[a,b]=[1.8618,1.3461]  1st kind complete and 2nd complete
a=ellipke([0.509,0.609])   %a=[1.8618,1.9594]    looks like a map operation on the list with func=1st-kind-complete
$\endgroup$
1
$\begingroup$

I'd just use Python:

from math import  ceil
from numpy import sqrt
from scipy.special import ellipk
def K(x):
    return ellipk(x*x)

omega_ratio = 0.509
e_delta_ratio = 0.00349

N = K(omega_ratio)*K(sqrt(1-e_delta_ratio*e_delta_ratio))/(K(e_delta_ratio)*K(sqrt(1-omega_ratio*omega_ratio)))
print(N)   # 3.5420113113759575
print(ceil(N))  # 4

which yields

3.5420113113759575

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.