function[m,n,k,l]=morsemom(p,ga,be)
%MORSEMOM  Frequency-domain moments of generalized Morse wavelets.
%
%   [MP,NP]=MORSEMOM(P,GAMMA,BETA) computes the Pth order frequency-
%   domain moment M and energy moment N of the lower-order generalized 
%   Morse wavelet specified by parameters GAMMA and BETA.
%
%   The Pth moment and energy moment are defined as
%
%           mp = 1/(2 pi) int omega^p  psi(omega)     d omega 
%           np = 1/(2 pi) int omega^p |psi(omega)|.^2 d omega 
%
%   respectively, where omega is the radian frequency.  These are 
%   evaluated using the normalization that max(abs(psi(omega)))=2.
%
%   The input parameters must either be matrices of the same size, or
%   some may be matrices and the others scalars.   
%
%   [MP,NP,KP,LP]=MORSEMOM(...) also returns the Pth order cumulant KP and
%   the Pth order energy cumulant LP.
%
%   For details see
%
%       Lilly and Olhede (2009).  Higher-order properties of analytic 
%           wavelets.  IEEE Trans. Sig. Proc., 57 (1), 146--160.
%
%   See also MORSEWAVE, MORSEDERIV, MOM2CUM.
%
%   'morsemom --t' runs some tests.
%
%   Usage:  mp=morsemom(p,ga,be);
%           [mp,np]=morsemom(p,ga,be);
%           [mp,np,kp,lp]=morsemom(p,ga,be);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2007--2009 J.M. Lilly --- type 'help jlab_license' for details

if strcmp(p, '--t')
    morsemom_test,return
end


m=morsemom1(p,ga,be);

if nargout>1
    n=frac(morsea(ga,be).^2,morsea(ga,2*be)).*frac(1,2.^(frac(2*be+1+p,ga))).*morsemom1(p,ga,2*be);
    n=frac(2,2.^(frac(1+p,ga))).*morsemom1(p,ga,2*be);
end

if nargout>2
     for i=0:maxmax(p)
         mcell{i+1}=morsemom1(i,ga,be);
     end
     kcell=mom2cum(mcell);
     if length(p)==1 
        k=kcell{p+1};
     else
        k=zeros(size(m));
        for i=1:length(m)
            k(i)=kcell{i};
        end
     end
end


if nargout>3
     for i=0:maxmax(p)
          ncell{i+1}=frac(2,2.^frac(1+i,ga)).*morsemom1(i,ga,2*be);
     end
     lcell=mom2cum(ncell);
     if length(p)==1 
        l=lcell{p+1};
     else
        l=zeros(size(m));
        for i=1:length(m)
            l(i)=lcell{i};
        end
     end
end



function[m]=morsemom1(p,ga,be)
m=morsea(ga,be).*morsef(ga,be+p);

function[f]=morsef(ga,be)
%MORSEF  Returns the generalized Morse wavelet first moment "f".
%
%   F=MORSEF(GAMMA,BETA) returns the normalized first frequency-
%   domain moment "F_{BETA,GAMMA}" of the lower-order generalized 
%   Morse wavelet specified by parameters GAMMA and BETA.

f=frac(1,2*pi*ga).*gamma(frac(be+1,ga));

function[]=morsemom_test
morsemom_test_expression;
morsemom_test_numerical;
morsemom_test_cumulant;
morsemom_test_energycumulant;

function[]=morsemom_test_energycumulant 
ga=[3 0.5 9];be=[2.0000 15.1200 0.751];

[a,sigt,sigo]=morsebox(ga,be);
om=morsefreq(ga,be);
[m2,n2,k2,l2]=morsemom(2,ga,be);
reporttest('MORSEMOM second energy cumulant matches MORSEBOX',aresame(sqrt(l2./om./om),sigo,1e-6))

function[]=morsemom_test_expression
ga=(2:.1:10);
be=(1:.1:10);
[ga,be]=meshgrid(ga,be);

clear bool
for p=1:10
    [mp,np]=morsemom(p,ga,be);
    mp2=2*frac(exp(1)*ga,be).^(be./ga).*frac(1,2*pi*ga).*gamma(frac(be+1+p,ga));
    bool(p)=aresame(mp,mp2,1e-10);
end

reporttest('MORSEMOM analytic expression',all(bool))

function[]=morsemom_test_numerical

ga1=(2:2:12);
be1=(1:2:10);

%n=0;
omi=(0:.01:20)';
dom=omi(2)-omi(1);
                
[m1c,n1c,m2c,n2c]=vzeros(length(ga1),length(be1));
clear m1c n1c m2c n2c
for i=1:length(ga1)
    for j=1:length(be1)
        ompsi=morsefreq(ga1(i),be1(j));
        
        psi=morsea(ga1(i),be1(j)).*(omi.*ompsi).^be1(j).*exp(-(omi.*ompsi).^ga1(i));
        %n=n+1;
        %psiall(:,n)=psi;
        
        m1c(j,i)=frac(dom,2*pi).*ompsi.^2.*vsum(omi.*psi,1);
        n1c(j,i)=frac(dom,2*pi).*ompsi.^2.*vsum(omi.*psi.^2,1);
        m2c(j,i)=frac(dom,2*pi).*ompsi.^3.*vsum(omi.^2.*psi,1);
        n2c(j,i)=frac(dom,2*pi).*ompsi.^3.*vsum(omi.^2.*psi.^2,1);
    end
end

[ga,be]=meshgrid(ga1,be1);
[m1,n1]=morsemom(1,ga,be);
[m2,n2]=morsemom(2,ga,be);

tol=2*10.^(-2);
bool=aresame(m1,m1c,tol)&&aresame(n1,n1c,tol)&&aresame(m2,m2c,tol)&&aresame(n2,n2c,tol);

reporttest('MORSEMOM numerical calculation',all(bool))

function[]=morsemom_test_cumulant

ga1=(2:2:12);
be1=(1:2:10);
[ga,be]=meshgrid(ga1,be1);

[mo,no,ko]=morsemom(0,ga,be);
[m1,n1,k1]=morsemom(1,ga,be);
[m2,n2,k2]=morsemom(2,ga,be);
[m3,n3,k3]=morsemom(3,ga,be);
    
kob=log(mo);
k1b=frac(m1,mo);
k2b=frac(m2,mo)-frac(m1,mo).^2;
k3b=2*frac(m1.^3,mo.^3)-3.*frac(m1,mo).*frac(m2,mo)+frac(m3,mo);


reporttest('MORSEMOM first three cumulants',aresame(ko,kob,1e-10)&&aresame(k1,k1b,1e-10)&&aresame(k2,k2b,1e-10)&&aresame(k3,k3b,1e-10))
