ParaHermIDFT

Rt = ParaHermIDFT(Rf,Nfft);

  Rt = ParaHermIDFT(Rf) returns the IDFT of a parahermitian matrix stored in
  Rf, with indices running from DC to 2pi. The matrix Rf is of dimension MxMxL,
  where M is the spatial dimension, and L the number of sample points on the
  unit circle.

  The function returns an odd-length time-domain version, whereby the length
  is appropriate trimmed in case the time domain support is shorter than the 
  number of sample points.

  Input parameters:
    Rf         cross-spectral density matrix
 
  Output parameter:
    Rt         space-time covariance-type matrix

0001 function Rt = ParaHermDFT(Rf);
0002 %Rt = ParaHermIDFT(Rf,Nfft);
0003 %
0004 %  Rt = ParaHermIDFT(Rf) returns the IDFT of a parahermitian matrix stored in
0005 %  Rf, with indices running from DC to 2pi. The matrix Rf is of dimension MxMxL,
0006 %  where M is the spatial dimension, and L the number of sample points on the
0007 %  unit circle.
0008 %
0009 %  The function returns an odd-length time-domain version, whereby the length
0010 %  is appropriate trimmed in case the time domain support is shorter than the
0011 %  number of sample points.
0012 %
0013 %  Input parameters:
0014 %    Rf         cross-spectral density matrix
0015 %
0016 %  Output parameter:
0017 %    Rt         space-time covariance-type matrix
0018 
0019 % S. Weiss, 18/4/2021
0020 
0021 [M,~,L] = size(Rf);
0022 Threshold = 10^(-14);
0023 
0024 % apply inverse DFT
0025 R_td = ifft(Rf,L,3);
0026 
0027 % rearrange time domain data
0028 if mod(L,2)==0,             % even length
0029   dummy = zeros(M,M,L+1);
0030   dummy(:,:,L/2+1:L+1) = R_td(:,:,1:L/2+1);
0031   dummy(:,:,1:L/2) = R_td(:,:,L/2+1:L);
0032 else                        % odd length
0033   dummy = zeros(M,M,L);
0034   dummy(:,:,(L+1)/2:L) = R_td(:,:,1:(L+1)/2);
0035   dummy(:,:,1:(L-1)/2) = R_td(:,:,(L+1)/2+1:L);
0036 end;
0037   
0038 % trim leading and trailing zeros
0039 L = size(dummy,3);
0040 Power = zeros((L+1)/2,1);
0041 for l = 1:(L+1)/2,
0042    Power(l) = norm(dummy(:,:,l),'fro');
0043 end;
0044 Index = max(find(cumsum(Power)<Threshold));
0045 Rt = dummy(:,:,Index+1:L-Index);