Fourier Transform: Inverse FFT of Positive Definite Power Spectrum Not Giving Positive Definite Function

Hello there,

I am attempting to write a Fourier transform "round trip" in 2D to obtain a real, positive definite covariance function. This is the following workflow:

1. Calculate covariance map. This map will be symmetric but it will NOT have all positive values.
2. Take the FFT of the covariance. This will give a power spectrum that is REAL, but NOT positive definite because the covariance is real but not positive definite.
3. Do some smoothing on the power spectrum so that all values are either positive or zero. This power spectrum is now positive definite and symmetric.
4. Take the IFFT of the smoothed, positive definite, real power spectrum. This should give a back-transformed, positive definite, real covariance map.

The purpose of this work is to apply Bochner's theorem, which states that a function positive in the frequency domain should be positive definite in the spatial domain. However, after making all power spectrum values either positive or zero, the back-transformed Covariance map is NOT positive definite, as it still has some negative values. Why is this the case? By ensuring that all power spectrum values are real and non-negative, shouldn't the back-transformed Covariance also be real and non-negative?

Here is my round-trip code. h1 and h2 are the x and y lags, Cov is the input covariance, and k is the degree of nxn smoothing on the spectrum (need not focus on s and countpairs, as they are for sparse data).

function [Cov,CovIFFT,fomega,fomega_s,fomega_ss] = FFT2D(h1,h2,Cov,k,s,countpairs) %k is degree of smoothing, h1 is x lags, h2 is y lags (enter -h:h)

if s==1 %if sparse data, smooth covariances first
    [X_smoothing, Y_smoothing] = meshgrid(h1,h2);
    Cov = SmoothingCovNaN(Cov,countpairs,X_smoothing,Y_smoothing,0.1);
end

Cov_shift=ifftshift(Cov); %shift coordinates so FFT is real
fomega_unshifted=fft2(Cov_shift); %take FFT
fomega_unshifted_real=real(fomega_unshifted); %remove negligible imaginary component 
fomega=fftshift(fomega_unshifted_real); %shift FFT to center; this is the data we will do post-processing on


%SCHEME 1: WINDOW SMOOTHING, SET NEGATIVE VALUES TO 0
% if k==0 %dont perform smoothing if k==0
%     fomega_s=fomega;
% else
%     fomega_s=SmoothingMovingAverage(fomega,k); %kxk window smoothing
% end

%SCHEME 2: SHIFT ALL VALUES BY MIN SO SMALLEST VALUE IS 0
[fomega_s]=NoNegShift(fomega);

[fomega_ss] = SmoothingSum1DS(fomega_s,h1,h2); %exponential standardize to 1 smoothing

% fomega_s=fomega; %REMOVE -- JUST FOR TESTING PURPOSES
% fomega_ss=fomega_s; %REMOVE -- JUST FOR TESTING PURPOSES


fomega_ss_unshifted=ifftshift(fomega_ss); %go back to unshifted image after post-processing to ensure IFFT is real

CovIFFT_unshifted=ifft2(fomega_ss_unshifted); %back transform unshifted image to get real covariance result

CovIFFT=real(fftshift(CovIFFT_unshifted)); %now that unshifted image has undergone IFFT, shift it back. This is the final output! (okay to take real because imaginary component is small)
end

Here are some example result figures:

1. Covariance map

2. Density spectrum (unsmoothed, not positive definite)

3. Smoothed density spectrum (positive definite, done by shifting all values by minimum value)

4. Smoothed density spectrum (positive definite, unit sum)

5. Back-transformed covariance map (SHOULD be positive definite because smoothed density spectrum is positive definite, but is not!)

Any help would be dearly appreciated. Thank you!

Corey

[Image Analyst]

What's your code? You forgot to include it. Are you smoothing the power spectrum, which is the magnitude squared of the Fourier Transform? Or just the real part of the FT? Or both the real and imaginary part of the FT?

To do a round trip intact you need to not touch the FT while it's in the frequency domain. If you inverse FT just the real part, you will get a Hermitian back in the spatial (or time) domain. But if you inverse FT the complex FT you MAY get close to the original back but it depends on how much you altered the signal while in the frequency domain.

[Corey Hoydic]

Hello there,

I have added the round-trip code to my original post. I can provide individual function codes if you would like, but the descriptions should suffice.

Since the input covariance function is real, the FFT of the covariance only has a real component. Thus, I am smoothing the real-valued spectrum values.

I have also edited my original post to include some samplefigures for (1) covariance map, (2) density spectrum, (3) positive definite density spectrum, (4) positive definite unit sum density spectrum, and (5) back-transformed covariance map, which should have all positive values but does not. Please do let me know if you need any other information to help. I appreciate it.

[David Goodmanson]

Hi Corey,

A = cov(rand(20,6))
Ashift = ifftshift(A);
B = fft2(Ashift)

A =
    0.0729   -0.0110   -0.0062   -0.0128   -0.0136   -0.0162     % real, symmetric
   -0.0110    0.0784   -0.0212    0.0250    0.0177   -0.0076
   -0.0062   -0.0212    0.0683    0.0040    0.0038    0.0127
   -0.0128    0.0250    0.0040    0.0850    0.0247    0.0142
   -0.0136    0.0177    0.0038    0.0247    0.0439    0.0101
   -0.0162   -0.0076    0.0127    0.0142    0.0101    0.0960

B =

  Columns 1 through 5

   0.4915 + 0.0000i   0.1058 - 0.0460i  -0.0161 + 0.0024i   0.1698 + 0.0000i  -0.0161 - 0.0024i
   0.1058 - 0.0460i   0.0675 + 0.0060i   0.0267 + 0.0272i   0.0722 - 0.0644i  -0.1007 + 0.0257i
  -0.0161 + 0.0024i   0.0267 + 0.0272i  -0.0918 - 0.0785i  -0.0619 - 0.0891i   0.4740 + 0.0000i
   0.1698 + 0.0000i   0.0722 - 0.0644i  -0.0619 - 0.0891i   0.4595 + 0.0000i  -0.0619 + 0.0891i
  -0.0161 - 0.0024i  -0.1007 + 0.0257i   0.4740 - 0.0000i  -0.0619 + 0.0891i  -0.0918 + 0.0785i
   0.1058 + 0.0460i   0.3840 + 0.0000i  -0.1007 - 0.0257i   0.0722 + 0.0644i   0.0267 - 0.0272i

  Column 6

   0.1058 + 0.0460i
   0.3840 + 0.0000i
  -0.1007 - 0.0257i
   0.0722 + 0.0644i
   0.0267 - 0.0272i
   0.0675 - 0.0060i

B is definitely complex. What am I missing?

[Corey Hoydic]

Hello David,

Please take a quick look at the following:

https://blogs.mathworks.com/steve/2010/07/16/complex-surprises-from-fft/

Since I am dealing with spatial data ranging from 2h lags -h to +h, I need to first use ifftshift() on the covariance function so that the resultant fft(Cov) is real-valued. It should hold that a real-valued covariance function should have a real-valued Fourier transform.

Corey

[David Goodmanson]

Hi Corey, A real symmetric function, as in Steve's example. has a real symmetric fourier transform. I assume your cov matrix has a similar form to A above, which is symmetric in the matrix sense but has no addtional symmetry, is that the case? If it is, fft2 operates on all the rows, then all the columns (or vice versa). Neither the individual rows or individual columns have the necessary form to produce a real fft.

[Corey Hoydic]

Hello David,

Please see an example matrix for the symmetries. You can see that the result is indeed real

A =

    0.5997    0.6909    0.7844    0.8178    0.7765
    0.6854    0.7865    0.8904    0.8872    0.8089
    0.7736    0.8854    1.0000    0.8854    0.7736
    0.8089    0.8872    0.8904    0.7865    0.6854
    0.7765    0.8178    0.7844    0.6909    0.5997
    
fft2(ifftshift(A))
ans =

   19.7733    0.9676    0.0200    0.0200    0.9676
    0.8892    0.7805   -0.1590    0.2828   -0.3881
    0.0192   -0.1706    0.1510   -0.0676    0.2883
    0.0192    0.2883   -0.0676    0.1510   -0.1706
    0.8892   -0.3881    0.2828   -0.1590    0.7805f

fftshift(ans)

ans =

    0.1510   -0.1706    0.0192    0.2883   -0.0676
   -0.1590    0.7805    0.8892   -0.3881    0.2828
    0.0200    0.9676   19.7733    0.9676    0.0200
    0.2828   -0.3881    0.8892    0.7805   -0.1590
   -0.0676    0.2883    0.0192   -0.1706    0.1510

[David Goodmanson]

Hi Corey, I agree that this kind of symmetry, A = A + fliplr(flipud(A)), does give a real result.