StollLab · stestoll · Mar 2, 2021 · Mar 5, 2021 · Mar 5, 2021 · Mar 6, 2021
diff --git a/docsrc/fftso3.html b/docsrc/fftso3.html
@@ -0,0 +1,127 @@
+<!DOCTYPE html>
+<html>
+<head>
+   <meta charset="utf-8">
+   <link rel="icon" href="img/eslogo196.png">
+   <link rel="stylesheet" type="text/css" href="style.css">
+   <link rel="stylesheet" href="highlight/matlab.css">
+   <script src="highlight/highlight.pack.js"></script>
+   <script>hljs.initHighlightingOnLoad();</script>
+   <title>fftso3</title>
+</head>
+
+<body>
+
+<header>
+<ul>
+<li><img src="img/eslogo42.png">
+<li class="header-title">EasySpin
+<li><a href="index.html">Documentation</a>
+<li><a href="references.html">Publications</a>
+<li><a href="http://easyspin.org" target="_blank">Website</a>
+<li><a href="http://easyspin.org/academy" target="_blank">Academy</a>
+<li><a href="http://easyspin.org/forum" target="_blank">Forum</a>
+</ul>
+</header>
+
+<section>
+
+<div class="functitle">fftso3</div>
+
+<p>
+Fast Fourier transform over the rotation group SO(3).
+</p>
+
+<!-- ============================================================= -->
+<div class="subtitle">Syntax</div>
+
+<pre class="matlab">
+c = fftso3(fcn,Lmax)
+[LMK,vals] = fftso3(fcn,Lmax)
+[c,LMK,vals] = fftso3(fcn,Lmax)
+</pre>
+
+<!-- ============================================================= -->
+<div class="subtitle">Description</div>
+
+<p>
+<code>fftso3</code> takes the function specified by the function handle <code>fcn</code> and calculates its Fourier coefficients over the space of rotations.
+</p>
+
+<p>
+<code>fcn</code> is a function handle for a function <code>fcn(alpha,beta,gamma)</code> defined over the three <a href="eulerangles.html">Euler angles</a> <code>alpha</code>, <code>beta</code>, and <code>gamma</code>, in radians. <code>fcn</code> should be able to handle array inputs for <code>alpha</code> and <code>gamma</code>.
+</p>
+
+<p>
+The Fourier expansion is in terms of Wigner functions
+</p>
+
+<div class="eqn">
+    <img src="eqn/fftso31.png" alt="[eqn]"><!--MATH
+    \begin{equation*}
+      f(\alpha,\beta,\gamma) =
+      \sum_{L=0}^\infty
+      \sum_{M=-L}^{L}
+      \sum_{K=-L}^{L}
+      f^L_{MK}
+      D^L_{MK}(\alpha,\beta,\gamma)
+    \end{equation*}
+    -->
+</div>
+
+<p>
+<code>Lmax</code> is the maxmium value of L.
+</p>
+
+<p>
+<code>c</code> is a cell array with <code>Lmax</code> elements. <code>c{L+1}</code> contains the (2L+1)x(2L+1) matrix of coefficients for a given L, with both M and K in ascending order. In other words, <img src="eqn/fftso32.png" alt="[eqn]"><!--MATH$f^L_{MK}$--> is given by <code>c{L+1}(L+1+M,L+1+K)</code>.
+</p>
+
+<p>
+<code>LMK</code> is a list of (L,M,K) triplets of the coefficients with non-zero values, with one triplet per row. <code>vals</code> is the corresponding list of values.
+</p>
+
+
+<!-- ============================================================= -->
+<div class="subtitle">Examples</div>
+
+<p>
+Here is a simple example that takes a Wigner function and returns the Fourier transform. Since the function is a single Wigner function, only a single coefficient is non-zero.
+</p>
+
+<pre class="matlab">
+f = @(a,b,c)10*wignerd([1 1 -1],a,b,c);
+c = fftso3(f,2);
+c{2}
+</pre>
+<pre class="mloutput">
+ans =
+    0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
+    0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
+   10.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
+</pre>
+
+
+<!-- ============================================================= -->
+<div class="subtitle">Algorithm</div>
+
+<p>
+See Kostelec & Rockmore, FFTs on the Rotation Group, J. Fourier Anal. Appl., 2008, 14, 145-179, <a href="https://doi.org/10.1007/s00041-008-9013-5">https://doi.org/10.1007/s00041-008-9013-5</a>.
+</p>
+
+<!-- ============================================================= -->
+<div class="subtitle">See also</div>
+
+<p>
+<a class="esf" href="erot.html">erot</a>,
+<a class="esf" href="eulang.html">eulang</a>,
+<a class="esf" href="wignerd.html">wignerd</a>
+</p>
+
+<hr>
+</section>
+
+<footer></footer>
+
+</body>
+</html>
diff --git a/docsrc/funcsalphabet.html b/docsrc/funcsalphabet.html
@@ -81,6 +81,7 @@
 <tr><td><a href="fdaxis.html">fdaxis</a></td><td>Frequency domain axis for FFT</td></tr>
 <tr><td><a href="fieldmod.html">fieldmod</a></td><td>Field modulation of EPR absorption spectra</td></tr>
 <tr><td><a href="gammae.html">gammae</a></td><td>Electron gyromagnetic ratio</td></tr>
+<tr><td><a href="fftso3.html">fftso3</a></td><td>Fast Fourier transform on the rotation group SO(3)</td></tr>
 <tr><td><a href="garlic.html">garlic</a></td><td>Liquid-state cw EPR spectra</td></tr>
 <tr><td><a href="gaussian.html">gaussian</a></td><td>Gaussian line shape</td></tr>
 <tr><td><a href="gfree.html">gfree</a></td><td>g value of the free electron</td></tr>

diff --git a/docsrc/funcscategory.html b/docsrc/funcscategory.html
@@ -199,6 +199,7 @@
 
 <tr><td colspan="2" class="grouptitle">Angular momentum</td></tr>
 <tr><td><a href="clebschgordan.html">clebschgordan</a></td><td>Clebsch-Gordan coefficients</td></tr>
+<tr><td><a href="fftso3.html">fftso3</a></td><td>Fast Fourier transform on the rotation group SO(3)</td></tr>
 <tr><td><a href="plegendre.html">plegendre</a></td><td>Legendre polynomials and Associated Legendre polynomials</td></tr>
 <tr><td><a href="spherharm.html">spherharm</a></td><td>Spherical harmonics</td></tr>
 <tr><td><a href="wigner3j.html">wigner3j</a></td><td>Wigner 3-j symbols</td></tr>

diff --git a/docsrc/releases.html b/docsrc/releases.html
@@ -61,6 +61,7 @@ <h1>Changes from release to release</h1>
 <li>The syntax for obtaining subspectra from simulation function has changed, the field name in the <code>Opt</code> structure has changed from <code>Opt.Output</code> to <code>Opt.separate</code> and has more options: <code>'components'</code>, <code>'transitions'</code>, <code>'sites'</code> and <code>'orientations'</code>. There is extended support for returning subspectra across all EasySpin simulation functions.</li>
 <li>For baseline correction with <a class="esf" href="basecorr.html">basecorr</a>, a fitting region can now be provided.</li>
 <li><a class="esf" href="eulang.html">eulang</a> now calculates Euler angles for a near-orthogonal matrix by first orthogonalizing it using singular-value decomposition.</li>
+<li><a class="esf" href="fftso3.html">fftso3</a> performs a Fast Fourier Transform over the space of rotations, SO(3).</li>
 </ul>
 
 <p>Major bug fixes</p>

diff --git a/easyspin/chili.m b/easyspin/chili.m
@@ -1007,7 +1007,11 @@
     % Calculate sqrt(Peq) vector
     Opt_.useSelectionRules = Opt.useStartvecSelectionRules;
     Opt_.PeqTolerances = Opt.PeqTol;
-    [sqrtPeq,nInt] = chili_eqpopvec(Basis,Potential,Opt_);
+    if Opt.useLMKbasis
+      [sqrtPeq,nInt] = chili_eqpopvec_fast(Basis,Potential,Opt_);
+    else
+      [sqrtPeq,nInt] = chili_eqpopvec(Basis,Potential,Opt_);
+    end
     % Set up in full product basis, then prune
     StartVector = kron(sqrtPeq,SdetOp(:)/norm(SdetOp(:)));
     StartVector = StartVector(keep);

diff --git a/easyspin/fftso3.m b/easyspin/fftso3.m
@@ -0,0 +1,185 @@
+% fftso3   Fast Fourier transform over the rotation group SO(3)
+%
+%   c = fftso3(fcn,Lmax)
+%   [LMK,vals] = fftso3(fcn,Lmax)
+%   [c,LMK,vals] = fftso3(fcn,Lmax)
+%   ___ = fftso3(fcn,Lmax,mode)
+%
+% Takes the function fcn(alpha,beta,gamma) defined over the Euler angles
+% 0<=alpha<2*pi, 0<=beta<=pi, 0<=gamma<2*pi and calculates the coefficients
+% c^L_MK of its tuncated expansion in terms of Wigner D-functions, D^L_MK:
+%
+%    fcn(alpha,beta,gamma) = sum_(L,M,K) c^L_MK * D^L_MK(alpha,beta,gamma)
+%
+% where L = 0,..,Lmax, M = -L,..,L, K = -L,..,L.
+%
+% Input:
+%   fcn   function handle for function fcn(alpha,beta,gamma), with angles
+%         in units of radians; it should be vectorized and accept array
+%         inputs for alpha and gamma
+%   Lmax  maximum L for expansion
+%   mode  beta grid mode, either 'GL' for Gauss-Legendre (Lmax+1 points)(default)
+%         or 'DH' for Driscoll-Healy (2*Lmax+2 points)
+%
+% Output:
+%   c     (Lmax+1)-element cell array of Fourier coefficients, L = 0 .. Lmax
+%         c{L+1} contains the (2L+1)x(2L+1)-dimensional matrix of expansion
+%         coefficients c^L_MK, with M and K ordered in ascending order: M,K =
+%         -L,-L+1,..,L.
+%   LMK   L, M, and K values of the non-zero coefficients
+%   vals  values of the non-zero coefficients
+
+% References:
+% 1) D. K. Maslen, D. N. Rockmore
+%    Generalized FFT's - A survey of some recent results
+%    Groups and Computation II
+%    DIMACS Series in Discrete Mathematics and Theoretical Computer
+%    Science, vol.28
+%    L. Finkelstein, W. M. Kantor (editors)
+%    American Mathematical Society, 1997
+%    https://doi.org/10.1007/s00041-008-9013-5
+% 2) P. J. Kostelec, D. N. Rockmore
+%    FFTs on the Rotation Group
+%    J.Fourier.Anal.Appl. 2008, 14, 145/179
+%    https://doi.org/10.1007/s00041-008-9013-5
+% 3) Z. Khalid, S. Durrani, R. A. Kennedy, Y. Wiaux, J. D. McEwen
+%    Gauss-Legendre Sampling on the Rotation Group
+%    IEEE Signal Process. Lett. 2016, 23, 207-211
+%    https://doi.org/10.1109/LSP.2015.2503295
+
+function varargout = fftso3(fcn,Lmax,betamode)
+
+if nargin==0
+  help(mfilename);
+end
+
+if nargin~=2 && nargin~=3
+  error('Two inputs are required (fcn, Lmax).');
+end
+
+if nargin<3
+  betamode = 'GL';
+end
+
+if numel(Lmax)~=1 || mod(Lmax,1) || Lmax<0
+  error('Second input (Lmax) must be a non-negative integer.');
+end
+
+B = Lmax + 1; % bandwidth (0 <=L < B)
+
+% Set up sampling grid over (alpha,gamma)
+alpha = 2*pi*(0:2*B-2)/(2*B-1); % radians
+gamma = 2*pi*(0:2*B-2)/(2*B-1); % radians
+
+% Get knots and weights for integration over beta
+switch betamode
+  case 'DH' % Driscoll-Healy
+    beta = pi*((0:2*B-1)+1/2)/(2*B); % radians
+    nbeta = numel(beta);
+    % Calculate beta weights
+    q = 0:B-1; % summation index
+    w = zeros(1,nbeta);
+    for b = 1:nbeta
+      w(b) = 2/B * sin(beta(b)) * sum(sin((2*q+1)*beta(b))./(2*q+1));
+    end
+  case 'GL' % Gauss-Legendre
+    [z,w] = gauleg(-1,1,B);
+    beta = acos(z);
+end
+nbeta = numel(beta);
+
+% Calculate 2D inverse FFT along alpha and gamma, one for each beta
+[alpha,gamma] = ndgrid(alpha,gamma); % 2D grid for function evaluations
+S = cell(1,nbeta);
+for b = 1:nbeta
+  f = fcn(alpha,beta(b),gamma); % evaluate function over (alpha,gamma) grid
+  S_ = (2*pi)^2*ifft2(f);
+  S{b} = fftshift(S_); % reorder to M,K = -(B-1), ..., 0, ..., B-1
+end
+
+% Calculate Wigner d-function values for all L and beta, with matrix exponential
+d = cell(Lmax+1,nbeta);
+for L = 0:Lmax
+  v = sqrt((1:2*L).*(2*L:-1:1))/2i; % off-diagonals of Jy matrix
+  Jy = diag(v,-1) - diag(v,1); % Jy matrix (M,K = -Lmax,-Lmax+1,...,+Lmax)
+  [U,mj] = eig(Jy);
+  mj = diag(mj).'; % equal to -L:L
+  for b = 1:nbeta
+    e = exp(-1i*beta(b)*mj);
+    d{L+1,b} = (U .* e) * U'; % equivalent to U*diag(e)*U', but slightly faster
+  end
+end
+
+% Calculate Fourier coefficients
+c = cell(Lmax+1,1);
+for L = 0:Lmax
+  c_ = zeros(2*L+1);
+  idx = (Lmax+1)+(-L:L); % to access the range M,K = -L:L in S
+  for b = 1:nbeta % run over all beta values
+    c_ = c_ + w(b)*d{L+1,b}.*S{b}(idx,idx);
+  end
+  c{L+1} = (2*L+1)/(8*pi^2)*c_;
+end
+
+% Remove numerical noise
+maxcoeff = cellfun(@(x)max(abs(x),[],'all'),c);
+threshold = 1e-13*max(maxcoeff);
+for L = 0:Lmax
+  idx = abs(c{L+1})<threshold;
+  c{L+1}(idx) = 0;
+end
+
+% Return full matrices or list of non-zero coefficients
+returnList = nargout>1;
+if returnList
+  % Convert to [L M K value] list
+  LMK = [];
+  vals = [];
+  for L = 0:Lmax
+    [idxK,idxM,vals_] = find(c{L+1}.'); % transpose to assure lexicographic order of L,M,K
+    nNonzero = numel(vals_);
+    if nNonzero>0
+      LMK = [LMK; repmat(L,nNonzero,1) idxM-L-1 idxK-L-1];
+      vals = [vals; vals_];
+    end
+  end
+end
+
+% Organize output
+switch nargout
+  case 1, varargout = {c};
+  case 2, varargout = {LMK,vals};
+  case 3, varargout = {c,LMK,vals};
+end
+
+end
+
+% Calculate Gauss-Legendre grid points and weights over inteval [x1,x2]
+% Press et al, Numerical Recipes in C, 2nd ed., p.152
+function [x,w] = gauleg(x1,x2,n)
+
+m = (n+1)/2;
+xm = (x2+x1)/2;
+xl = (x2-x1)/2;
+for i = 1:m
+  z = cos(pi*(i-0.25)/(n+0.5));
+  z1 = inf;
+  while abs(z-z1)>eps
+    p1 = 1;
+    p2 = 0;
+    for j = 1:n
+      p3 = p2;
+      p2 = p1;
+      p1 = ((2*j-1)*z*p2 - (j-1)*p3)/j;
+    end
+    pp = n*(z*p1-p2)/(z^2-1);
+    z1 = z;
+    z = z1 - p1/pp;
+  end
+  x(i) = xm - xl*z;
+  x(n+1-i)= xm + xl*z;
+  w(i) = 2*xl/((1-z^2)*pp^2);
+  w(n+1-i) = w(i);
+end
+
+end
diff --git a/easyspin/private/chili_eqpopvec.m b/easyspin/private/chili_eqpopvec.m
@@ -38,7 +38,7 @@
   Mp = Potential.M;
   Kp = Potential.K;
   lambda = Potential.lambda;
-else
+elseif isnumeric(Potential)
   if size(Potential,2)~=4
     error('Potential must have four columns (L, M, K, lambda)');
   end
@@ -48,6 +48,8 @@
   lambda = Potential(:,4);
 end
 
+% Process potential
+%--------------------------------------------------------------------------
 % Assure that potential contains only terms with nonnegative K, and nonnegative M
 % for K=0. The other terms needed to render the potential real-valued are
 % implicitly supplemented in the subfunction U(a,b,c).
@@ -92,8 +94,7 @@
   if numel(idx0)~=1
     error('Exactly one orientational basis function with L=M=K=0 is allowed.');
   end
-  sqrtPeq = zeros(nOriBasis,1);
-  sqrtPeq(idx0) = 1;
+  sqrtPeq = sparse(idx0,1,1,nOriBasis,1);
   return
 end
 
@@ -187,10 +188,11 @@
     sqrtPeq(b) = Int;
   end
 
-  sqrtPeq = sparse(sqrtPeq);
 
 end
 
+sqrtPeq = sparse(sqrtPeq);
+
   % General orientational potential function (real-valued)
   % (assumes nonnegative K, and nonnegative M for K=0, and real lambda for
   % M=K=0; other terms are implicitly supplemented)