Inferring parent distribution using phase of 2D discrete fourier transform of random image

Question

Suppose I have an NxM matrix with real entries $S_{n,m}$ which are iid with zero mean. Say $n \in [0, 1, ..., N-1], m \in [0, 1, ..., M-1]$.

From $S_{n,m}$ I can compute the 2D discrete fourier transform: $$\eta_{j,l} = \sum_{n=0}^{N-1} \sum_{m=0}^{M-1} S_{n,m}~e^{-2\pi i (nj/N + ml/M)} $$ which has amplitude matrix $F_{j,l}$ and phase matrix $\Phi_{j,l}$, defined from: $$ \eta_{j,l} = F_{j,l} e^{i\Phi_{j,l}}$$

Overall Question: Does the phase matrix $\Phi_{j,l}$ alone contain significant information about the distribution from which the $S_{n,m}$ were drawn? I have a situation where I can only know $\Phi_{j,l}$, and want to infer something about the distribution of the $S_{n,m}$.

Concrete Example (this is the kind of thing I'd like to do):

Suppose I have two phase matrices $\Phi_{j,l}^{A}$ and $\Phi_{j,l}^{B}$. I know that one of them was generated from an $S_{n,m}$ matrix drawn from the t-distribution with high degrees-of-freedom = 1000 (so it's close to the normal distribution), while the other was drawn from a t-distribution with degrees-of-freedom = 2 (very heavy tails). I don't have the associated amplitude matrices.

Is there any hope of deriving a procedure to determine which distribution was used to produce which phase matrix, with a better than 50% success rate? If so, what kind of methods can be used?

EDIT (following some further investigation):

Experimentally the answer seems to be 'yes, at least in some cases' -- but I would appreciate an explanation of why.

Firstly, since the phase is unaffected by positive rescaling of $S_{n,m}$, we could never hope to detect differences in the scale of two distributions used to generate them. But it seems we can detect other differences.

Experimentally, say $Y_{n,m}$ is the inverse DFT of $e^{i\Phi_{j,l}}$ (like looking at the original image with all amplitudes set to 1). The code below investigates the distribution of this in the example above. It appears the Kurtosis of $Y_{n,m}$ is generally much higher for the low degrees-of-freedom t-distribution -- so for this particular problem, we would have a good chance of determining the correct parent distribution from the phase alone.

Can anyone explain why the disribution of $Y_{n,m}$ retains some information on the kurtosis of the distribution? This is not obvious to me.

kurtosis<-function(x){
    n = length(x)
    n * sum((x - mean(x))^4)/(sum((x - mean(x))^2)^2)
}


N = 20
M = 20

nreps = 20

kurtosis1000_store = rep(NA, nreps)
kurtosis2_store = rep(NA, nreps)

pdf('Inverse_dft_phase_plots.pdf', width=8, height=8)
for(i in 1:nreps){
    # This S uses the t-distribution with a large df = 1000
    S_t1000 = matrix( rt(400, df=1000), ncol = M, nrow = N)

    # This S uses the t-distribution with a small df = 2 (finite mean but not
    # variance)
    S_t2 = matrix( rt(400, df=2), ncol = M, nrow = N)

    # Phase
    Phase_t1000 = Arg(fft(S_t1000))
    Phase_t2 = Arg(fft(S_t2))

    # If we only have these phase matrices, could we make a statistical test to
    # determine which came from the low-df / high-df distributions?
    par(mfrow=c(2,2))

    # Try to look at the inverse DFT exp(i*Phase)
    # There seem to be obvious differences (the one associated with df=2 has
    # heavier outliers, greater kurtosis, etc)
    t1000 = Re(fft(exp(1i*Phase_t1000), inverse=T)/(N*M))
    t2 = Re(fft(exp(1i*Phase_t2), inverse=T)/(N*M)) 
    image(t1000, main='df = 1000')
    image(t2, main='df = 2')
    hist(t1000, main='df = 1000')
    hist(t2, main='df = 2')

    kurtosis1000_store[i] = kurtosis(t1000)
    kurtosis2_store[i] = kurtosis(t2)
}
dev.off()

This leads to

> summary(kurtosis1000_store)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.680   2.817   2.921   2.946   3.029   3.471 

> summary(kurtosis2_store)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  5.548  19.880  56.040  86.500 120.500 310.400 

Answer 1

I thought about this some. I don't think I've figured everything out, but I may as well get my thoughts on paper.

Let me start with the case of a continuous signal, and then reason by analogy. So let's say that for a function $f: R \rightarrow C$, I write $\hat{f}$ for its fourier transform, and $\check{f}$ for its inverse transform.

Now, I'm going to model your situation by taking two real valued functions $f_1, f_2 : R \rightarrow R$, and analyzing what the relationship must be if the phase of their fourier transforms is always equal. That is, if:

$$ Arg(\hat{f_1}(\xi)) = Arg(\hat{f_2}(\xi)) $$

for all frequencies $\xi$. The relationship above holds if and only if:

$$ \hat{f_1}(\xi) = r(\xi)\hat{f_2}(\xi) $$

for $r$ a real valued function. Taking the inverse transform of this relationship, and using the convolution identity:

$$ f_1(t) = (\check{r} * f_2)(t) $$

Now, $h = \check{r}$ is the inverse transform of a real valued function, and hence must satisfy $h(-t) = \overline{h(t)}$. It's clear that all such functions, when decomposed into their real and imaginary parts:

$$ h(t) = e(t) + o(t)i $$

must have $e$ and even function and $o$ and odd function. But then:

$$ f_1(t) = (e * f_2)(t) + i(o * f_2)(t) $$

Since $f_1$ is real valued, the imaginary part of this expression must vanish, and hence we are led to the relationship:

$$ f_1(t) = (e * f_2)(t) $$

So $f_1$ and $f_2$ are related by convolution with an even function. It's easy to reverse these steps, and convince yourself that this is actually an equivalence:

$$ Arg(\hat{f_1}(\xi)) = Arg(\hat{f_2}(\xi)) \iff f_1(t) = (e * f_2)(t) $$

The upside is that you can tell some things about a function $f$ from only its phase, without knowing its amplitudes. The things you can discern are exactly those properties that are preserved when convolving with an even function. A couple that seem geometrically evident are:

• If $f$ has thick tails, so will its convolution with an even function.
• If $f$ is skew, so will be its convolution with an even function.

Ok, so to the discrete case. The two essential properties I used in the argument were the convolution property, and the fact that the transform of a real valued function is hermitian (satisfies $h(-t) = \overline{h(t)}$). Both of these have discrete analogs. The dft convolution property is well known, and the analog of the hermitian property is

$$ \hat{x}_{N-k} = \overline{\hat{x}_{k}} $$

for real sequences $x_k$. Using these properties, the argument should go through, and we can conclude that sequences whose discrete fourier transforms agree up to a phase must differ by convolution with an even sequence.

Now I'll take a stab at explaining your kurtosis observation, but I'm not so sure it's correct. When you randomly fill in elements of the array by selecting from the distribution with long tails, it is probable that there is an asymmetry in the array caused by an extreme observation on one side or another of the center. Convolving with an "even" array preserves this asymmetry. This is why high kurtosis measures are preserved.