Home > manopt > manifolds > complexcircle > realphasefactory.m

realphasefactory

PURPOSE ^

Returns a manifold struct to optimize over phases of fft's of real signals

SYNOPSIS ^

function M = realphasefactory(n, z0, zmax)

DESCRIPTION ^

 Returns a manifold struct to optimize over phases of fft's of real signals

 function M = realphasefactory(n)
 function M = realphasefactory(n, z0)
 function M = realphasefactory(n, z0, zmax)

 If x is a real vector of length n, then y = fft(x) is a complex vector
 which obeys certain symmetries. Specifically, for any integer k,

  y(1+mod(k, n)) = conj(y(1+mod(n-k, n)))

 The same holds for the phases of the Fourier transform z = sign(y).

 This factory returns a Manopt manifold structure which represents the set
 of complex vectors z of length n which could be the phases of the Fourier
 transform of a real signal of length n:

   abs(z) = 1   and   z(1+mod(k, n)) = conj(z(1+mod(n-k, n))) for each k.

 For k = 1, this readily implies that z(1) is +1 or -1, so that the set of
 possible z's is disconnected. To choose which connected component to work
 with, set the second input z0 to +1 or -1 (this is the sign of the mean
 of x). By default, z0 = 1.

 Furthermore, if n is even, then k = n/2 implies z(1+n/2) is +1 or -1 as
 well, thus further disconnecting the set of acceptable z's. To choose
 which component to work with, set the third input zmax to +1 or -1. By
 default, it is +1.

 The Riemannian manifold structure is the Riemannian submanifold
 structure from the embedding space R^2 x ... x R^2, i.e., the complex
 circles are identified with the unit circle in the real plane.
 Concretely, this means the inner product is <u, v>_z = real(u'*v).
 Tangent vectors at z are complex vectors of length n which notably
 satisfy z(1+0) = 0 and, if n is even, z(1+n/2) = 0.

 n must be integer and n >= 3 (for n = 1:2 the manifold has dimension 0).

 Extra functions available in M include M.up, M.down and M.downup. They
 allow to capture the symmetries concisely, as:

    M.up(z) == conj(M.down(z)).

 See in code for more details.

 See also complexcirclefactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = realphasefactory(n, z0, zmax)
0002 % Returns a manifold struct to optimize over phases of fft's of real signals
0003 %
0004 % function M = realphasefactory(n)
0005 % function M = realphasefactory(n, z0)
0006 % function M = realphasefactory(n, z0, zmax)
0007 %
0008 % If x is a real vector of length n, then y = fft(x) is a complex vector
0009 % which obeys certain symmetries. Specifically, for any integer k,
0010 %
0011 %  y(1+mod(k, n)) = conj(y(1+mod(n-k, n)))
0012 %
0013 % The same holds for the phases of the Fourier transform z = sign(y).
0014 %
0015 % This factory returns a Manopt manifold structure which represents the set
0016 % of complex vectors z of length n which could be the phases of the Fourier
0017 % transform of a real signal of length n:
0018 %
0019 %   abs(z) = 1   and   z(1+mod(k, n)) = conj(z(1+mod(n-k, n))) for each k.
0020 %
0021 % For k = 1, this readily implies that z(1) is +1 or -1, so that the set of
0022 % possible z's is disconnected. To choose which connected component to work
0023 % with, set the second input z0 to +1 or -1 (this is the sign of the mean
0024 % of x). By default, z0 = 1.
0025 %
0026 % Furthermore, if n is even, then k = n/2 implies z(1+n/2) is +1 or -1 as
0027 % well, thus further disconnecting the set of acceptable z's. To choose
0028 % which component to work with, set the third input zmax to +1 or -1. By
0029 % default, it is +1.
0030 %
0031 % The Riemannian manifold structure is the Riemannian submanifold
0032 % structure from the embedding space R^2 x ... x R^2, i.e., the complex
0033 % circles are identified with the unit circle in the real plane.
0034 % Concretely, this means the inner product is <u, v>_z = real(u'*v).
0035 % Tangent vectors at z are complex vectors of length n which notably
0036 % satisfy z(1+0) = 0 and, if n is even, z(1+n/2) = 0.
0037 %
0038 % n must be integer and n >= 3 (for n = 1:2 the manifold has dimension 0).
0039 %
0040 % Extra functions available in M include M.up, M.down and M.downup. They
0041 % allow to capture the symmetries concisely, as:
0042 %
0043 %    M.up(z) == conj(M.down(z)).
0044 %
0045 % See in code for more details.
0046 %
0047 % See also complexcirclefactory
0048 
0049 % This file is part of Manopt: www.manopt.org.
0050 % Original author: Nicolas Boumal, Feb. 2, 2017.
0051 % Contributors: joint work with Tamir Bendory, Zhizhen Zhao and Amit Singer
0052 % Change log:
0053 %
0054 %   July 20, 2017 (NB)
0055 %       The distance function is now more accurate. Improved logarithm
0056 %       accordingly.
0057 
0058     assert(n == round(n) && n >= 3, 'n must be an integer >= 3.');
0059     
0060     even_n = (round(n/2) == n/2);
0061     
0062     if ~exist('z0', 'var') || isempty(z0)
0063         z0 = 1;
0064     end
0065     if ~exist('zmax', 'var') || isempty(zmax)
0066         zmax = 1;
0067     end
0068     
0069     assert(z0 == 1 || z0 == -1, 'z0 must be +1 or -1.');
0070     assert(zmax == 1 || zmax == -1, 'zmax must be +1 or -1.');
0071 
0072     if even_n
0073         M.name = @() sprintf('Phases of fft''s of real signals of length %d (z0 = %d, zmax = %d)', n, z0, zmax);
0074     else
0075         M.name = @() sprintf('Phases of fft''s of real signals of length %d (z0 = %d)', n, z0);
0076     end
0077     
0078     M.dim = @() floor((n-1)/2);
0079     
0080     M.inner = @(z, v, w) real(v'*w);
0081     
0082     M.norm = @(z, u) norm(u);
0083     
0084     M.dist = @(z1, z2) norm(real(2*asin(.5*abs(z1 - z2))));
0085     
0086     M.typicaldist = @() pi*sqrt(n/2);
0087     
0088     % Special functions to ease working with the symmetries.
0089     down = @(u) u;
0090     up = @(u) u([1 ; (n:-1:2)']);
0091     downup = @(u) (down(u) + conj(up(u)))/2;
0092     M.down = down;
0093     M.up = up;
0094     M.downup = downup;
0095     
0096     M.proj = @proj;
0097     function pu = proj(z, u)
0098         duu = downup(u);
0099         pu = duu - real(duu .* conj(z)).*z;
0100         % Note that there is no need to enforce pu(1) = 0 or (if n is even)
0101         % pu(1+n/2) = 0 manually, since the IEEE standard ensures that the
0102         % above operation will be exact for those entries provided z(1)
0103         % (and possibly z(1+n/2) is +1 or -1, as should be the case.
0104     end
0105     
0106     M.tangent = M.proj;
0107     
0108     % For Riemannian submanifolds, converting a Euclidean gradient into a
0109     % Riemannian gradient amounts to an orthogonal projection.
0110     M.egrad2rgrad = M.proj;
0111     
0112     M.ehess2rhess = @ehess2rhess;
0113     function rhess = ehess2rhess(z, egrad, ehess, zdot)
0114         rhess = M.proj(z, ehess - real(downup(egrad) .* conj(z)).*zdot);
0115     end
0116     
0117     M.exp = @exponential;
0118     function y = exponential(z, v, t)
0119         
0120         if nargin == 2
0121             % t = 1;
0122             tv = v;
0123         else
0124             tv = t*v;
0125         end
0126 
0127         y = zeros(n, 1);
0128 
0129         nrm_tv = abs(tv);
0130         
0131         % We need to be careful for zero steps.
0132         mask = (nrm_tv > 0);
0133         y(mask) = z(mask).*cos(nrm_tv(mask)) + ...
0134                   tv(mask).*(sin(nrm_tv(mask))./nrm_tv(mask));
0135         y(~mask) = z(~mask);
0136         
0137     end
0138     
0139     M.retr = @retraction;
0140     function y = retraction(z, v, t)
0141         if nargin == 2
0142             % t = 1;
0143             tv = v;
0144         else
0145             tv = t*v;
0146         end
0147         y = sign(z+tv);
0148     end
0149 
0150     M.log = @logarithm;
0151     function v = logarithm(x1, x2)
0152         v = M.proj(x1, x2 - x1);
0153         di = real(2*asin(.5*abs(x1 - x2)));
0154         nv = abs(v);
0155         factors = di ./ nv;
0156         factors(di <= 1e-6) = 1;
0157         v = v .* factors;
0158     end
0159     
0160     M.hash = @(z) ['z' hashmd5( [real(z(:)) ; imag(z(:))] ) ];
0161     
0162     M.rand = @random;
0163     function z = random()
0164         z = sign(downup(randn(n, 1) + 1i*randn(n, 1)));
0165         z(1) = z0;
0166         if even_n
0167             z(1 + n/2) = zmax;
0168         end
0169     end
0170     
0171     M.randvec = @randomvec;
0172     function v = randomvec(z)
0173         v = M.proj(z, randn(n, 1) + 1i*randn(n, 1));
0174         v = v / norm(v);
0175     end
0176     
0177     M.lincomb = @matrixlincomb;
0178     
0179     M.zerovec = @(z) zeros(n, 1);
0180     
0181     M.transp = @(z1, z2, u) M.proj(z2, u);
0182     
0183     M.pairmean = @pairmean;
0184     function z = pairmean(z1, z2)
0185         z = sign(z1+z2);
0186     end
0187 
0188     % This vec/mat pair is an isometry which allows to switch between the
0189     % classical representation of tangent vectors---as complex vectors of
0190     % length n---to real vectors of length M.dim() whose entries are the
0191     % coordinates of the tangent vector in the basis 1i*z, for the first
0192     % half. A scaling of sqrt(2) is applied to ensure isometry, since
0193     % tangent vectors are represented with only half of their entries.
0194     I = 2 : floor((n+1)/2);
0195     if even_n
0196         middle = 0;
0197     else
0198         middle = [];
0199     end
0200     M.vec = @(z, u_mat) sqrt(2)*real(u_mat(I) .* conj(1i*z(I)));
0201     M.mat = @(z, u_vec) [0 ; u_vec.*(1i*z(I)) ; middle ; ...
0202                              flipud(conj(u_vec.*(1i*z(I))))]/sqrt(2);
0203     M.vecmatareisometries = @() true;
0204     
0205 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005