function[o1,o2,o3,o4]=ellconv(i1,i2,i3,i4,i5)% ELLCONV  Converts between time-varying ellipse representations.%%   [Y1,Y2,Y3,Y4]=ELLCONV(X1,X2,X3,X4,STR) converts a time-varying ellipse%   from one representation to another, as specified by the string STR.  %%   STR is of the form 'in2out' where the names 'in' and 'out' may be any %   of the following%  %       Name    Arguments            Description%       +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%         'xy'  [X,Y,PHIX,PHIY]      X and Y component form   %         'pn'  [P,N,PHIP,PHIN]      Positively / negatively rotating form%         'ab'  [A,B,THETA,PHI]      Semi-major / semi-minor axis form  %  %   thus [X,Y,PHIX,PHIY]=ELLCONV(A,B,THETA,PHI,'ab2xy') converts from %   semi-major / semi-minor axis form to X/Y component form. %  %   More specifically, a time-varying ellipse in a complex-valued time %   series is represented in one of three ways%   %         'xy'  Z = X COS(PHIX) + i Y * COS(PHIY)  %         'pn'  Z = P EXP(i PHIP) + i N EXP(-i PHIN) %         'ab'  Z = EXP(i THETA) [A COS(PHI) + i B * SIN(PHI)] %%   where 'i' is the square root of negative one.  For additional details,%   see Lilly and Gascard (2006).  %   _____________________________________________________________________%%   Ridge conversion%%   ELLCONV also converts a wavelet transform along a wavelet ridge into%   one of the ellipse forms listed above, as follows:   %%   [...]=ELLCONV(WX,WY,'xy2--') converts to the ellipse form specifed by %    '--' from WX and WY, the complex-valued analytic wavelet transforms %   of the X and Y time series components along a ridge.%  %   [...]=ELLCONV(WP,WN,'pn2--') similarly converts from WP and WN, the %   complex-valued analytic and anti-anaytic wavelets, respectively, of a %   complex-valued time series Z=X+iY.  %%   The wavelet transform pairs WX and WY, and WP and WN, are defined in%   Lilly and Gascard (2006) and are implemented by TRANSCONV.%   _____________________________________________________________________%%   See also ECCONV, ELLDIFF.%  %   Usage:  [a,b,theta,phi]=ellconv(U,V,phiu,phiv,'xy2ab');%           [P,N,phip,phin]=ellconv(U,V,phiu,phiv,'xy2pn'); %           [a,b,theta,phi]=ellconv(wu,wv,'xy2ab');%%   'ellconv --t' runs some tests%   __________________________________________________________________%   This is part of JLAB --- type 'help jlab' for more information%   (C) 2005--2006 J.M. Lilly --- type 'help jlab_license' for details    %         'kl'  [KAP,LAM,THETA,PHI]  Kappa / Lambda form%         'kl'  Z = KAP EXP(i THETA)     if strcmp(i1,'--t')  ellconv_test;returnendif nargin==5   str=i5;elseif nargin==3   str=i3;endellconv_checkstr(str);if nargin==5   eval(['[o1,o2,o3,o4]=ellconv_' str '(i1,i2,i3,i4);'])elseif nargin==3   eval(['[o1,o2,o3,o4]=ellconv_' str '(i1,i2);'])endfunction[]=ellconv_checkstr(str)ii=findstr(str,'2');if isempty(ii)  error(['STR should be of the form ''xy2pn''.'])endstr1=str(1:ii-1);str2=str(ii+1:end);if ~aresame(str1,'ab') && ~aresame(str1,'pn') && ~aresame(str1,'xy')  error(['Input parameter name ' str1 ' is not supported.'])endif ~aresame(str2,'ab') && ~aresame(str2,'pn') && ~aresame(str2,'xy')  error(['Output parameter name ' str2 ' is not supported.'])end function[X,Y,phix,phiy]=ellconv_xy2xy(X,Y,phix,phiy)%Sort out input arguments  if nargin==2     wu=X;   wv=Y;   X=abs(wu);   Y=abs(wv);   phix=angle(wu);   phiy=angle(wv);endphix=unwrangle(phix);phiy=unwrangle(phiy);function[P,N,phip,phin]=ellconv_pn2pn(P,N,phip,phin)  %Sort out input argumentsif nargin==2     wp=P;   wn=N;   P=abs(wp)*sqrt(2);   N=abs(wn)*sqrt(2);   phip=angle(wp);   phin=-angle(wn);endphip=unwrangle(phip);phin=unwrangle(phin);function[A,B,theta,phi]=ellconv_ab2ab(A,B,theta,phi)phi=unwrangle(phi);theta=unwrangle(theta);%Do nothing  function[A,B,theta,phi]=ellconv_xy2ab(X,Y,phix,phiy)if nargin==2  [X,Y,phix,phiy]=ellconv_xy2xy(X,Y);end[P,N,phip,phin]=ellconv_xy2pn(X,Y,phix,phiy);[A,B,theta,phi]=ellconv_pn2ab(P,N,phip,phin);function[X,Y,phix,phiy]=ellconv_ab2xy(A,B,theta,phi)[P,N,phip,phin]=ellconv_ab2pn(A,B,theta,phi);[X,Y,phix,phiy]=ellconv_pn2xy(P,N,phip,phin);function[P,N,phip,phin]=ellconv_ab2pn(A,B,theta,phi)  P=frac(A+B,2);N=frac(A-B,2);phip=unwrangle(phi+theta);phin=unwrangle(phi-theta);function[A,B,theta,phi]=ellconv_pn2ab(P,N,phip,phin)  if nargin==2  [P,N,phip,phin]=ellconv_pn2pn(P,N);endA=P+N;B=P-N;theta=unwrangle(phip/2-phin/2);phi=  unwrangle(phip/2+phin/2);function[X,Y,phix,phiy]=ellconv_pn2xy(P,N,phip,phin)if nargin==2  [P,N,phip,phin]=ellconv_pn2pn(P,N);endtheta=unwrangle(phip/2-phin/2);phi=  unwrangle(phip/2+phin/2);X=sqrt(squared(P)+squared(N)+2.*P.*N.*cos(2*theta));Y=sqrt(squared(P)+squared(N)-2.*P.*N.*cos(2*theta));phixprime=imlog(P.*rot(theta)+N.*rot(-theta));phiyprime=imlog(P.*rot(theta)-N.*rot(-theta))-pi/2;phix=unwrangle(phi+phixprime);phiy=unwrangle(phi+phiyprime);function[P,N,phip,phin]=ellconv_xy2pn(X,Y,phix,phiy)if nargin==2  [X,Y,phix,phiy]=ellconv_xy2xy(X,Y);endphia=(phix+phiy+pi/2)/2;phid=(phix-phiy-pi/2)/2;P=frac(1,2)*sqrt(squared(X)+squared(Y)+2.*X.*Y.*cos(2*phid));N=frac(1,2)*sqrt(squared(X)+squared(Y)-2.*X.*Y.*cos(2*phid));phip=unwrangle(phia+imlog(X.*rot(phid)+Y.*rot(-phid)));phin=unwrangle(phia+imlog(X.*rot(phid)-Y.*rot(-phid)));function[]=ellconv_testL=1000;Xo=abs(randn(L,1));Yo=abs(randn(L,1));phixo=2*pi*rand(L,1);phiyo=2*pi*rand(L,1);tol=1e-6;%XY to AB and back[a,b,theta,phi]=ellconv(Xo,Yo,phixo,phiyo,'xy2ab');[X,Y,phix,phiy]=ellconv(a,b,theta,phi,'ab2xy');mat=[X Y angle(rot(phix)) angle(rot(phiy))];mato=[Xo Yo angle(rot(phixo)) angle(rot(phiyo))];bool(1)=aresame(mat,mato,tol);%XY to PN and back[P,N,phip,phin]=ellconv(Xo,Yo,phixo,phiyo,'xy2pn');[X,Y,phix,phiy]=ellconv(P,N,phip,phin,'pn2xy');mat=[X Y angle(rot(phix)) angle(rot(phiy))];mato=[Xo Yo angle(rot(phixo)) angle(rot(phiyo))];bool(2)=aresame(mat,mato,tol);                 %XY to AB to PN and back[a,b,theta,phi]=ellconv(Xo,Yo,phixo,phiyo,'xy2ab');[P,N,phip,phin]=ellconv(a,b,theta,phi,'ab2pn');[a,b,theta,phi]=ellconv(P,N,phip,phin,'pn2ab');[X,Y,phix,phiy]=ellconv(a,b,theta,phi,'ab2xy');mato=[Xo Yo angle(rot(phixo)) angle(rot(phiyo))];bool(3)=aresame(mat,mato,tol);%XY to PN to AB and back[P,N,phip,phin]=ellconv(Xo,Yo,phixo,phiyo,'xy2pn');[a,b,theta,phi]=ellconv(P,N,phip,phin,'pn2ab');[P,N,phip,phin]=ellconv(a,b,theta,phi,'ab2pn');[X,Y,phix,phiy]=ellconv(P,N,phip,phin,'pn2xy');mato=[Xo Yo angle(rot(phixo)) angle(rot(phiyo))];bool(4)=aresame(mat,mato,tol);reporttest('ELLCONV four-parameter conversions', all(bool)); clear bool%XY to XY [X,Y,phix,phiy]=ellconv(Xo.*rot(phixo),Yo.*rot(phiyo),'xy2xy');mato=[Xo Yo angle(rot(phixo)) angle(rot(phiyo))];bool(1)=aresame(mat,mato,tol);%XY to PN to PN and back [X,Y,phix,phiy]=ellconv(Xo.*rot(phixo),Yo.*rot(phiyo),'xy2xy');[P,N,phip,phin]=ellconv(Xo,Yo,phixo,phiyo,'xy2pn');[P,N,phip,phin]=ellconv(P.*rot(phip)/sqrt(2),N.*rot(phin)/sqrt(2),'pn2pn');[X,Y,phix,phiy]=ellconv(P,N,phip,phin,'pn2xy');mato=[Xo Yo angle(rot(phixo)) angle(rot(phiyo))];bool(2)=aresame(mat,mato,tol);reporttest('ELLCONV two-parameter conversions', all(bool));