How to choose a phase for the deconvolution of an autocorrelation?

Question

Say I have a function, $C=C\left(x\right)$, whose fourier transform is denoted by $c=c\left(k\right)$, i.e. $C\left(x\right)=\sum_{k=-\infty}^{\infty}c\left(k\right)\chi\left(x\right)$, where $\chi\left(x\right)$ is some complex valued basis function.

I know that the autocorrelation is given by $M\left(\Delta x\right)=\sum_{k=-\infty}^{\infty}m\left(k\right)\chi\left(\Delta x\right)$, where $\Delta x = x_2-x_1$ and by the convolution theorem I know that $m\left(k\right)=\left(c\left(k\right)\right)^*c\left(k\right)$. This means that knowing the coefficients of the autocorrelation immediately gives me the squared magnitude of the original coefficients. This implies that any two functions $C'\left(x\right)=C\left(x-x'\right)$, that differ only by a translation, will possess the same autocorrelation, $M\left(\Delta x\right)$, and so if I am trying to obtain the $c\left(k\right)$ from the $m\left(k\right)$, then I can choose an arbitrary phase for the $c\left(k\right)$. However, I don't think I can choose the phase of $c\left(k\right)$ independently for each $k$, it seems like there needs to be some way to consistently choose the phase for all $k$ (which corresponds to specifying a translation for $C\left(x\right)$).

My question is: How do I choose a consistent phase for each of the coefficients for the original function, i.e. for all $k$?

Answer 1

For a real-valued signal $C$, the autocorrelation function $M$ is a real-valued even function and the power spectral density $m$ (Fourier transform of $M$) is a real-valued nonnegative even function. Now, $m(k) = c(k)c^*(k)$ where $c$ is the Fourier transform of $C$, as you correctly assert, but given only $m$ and no other information about $c$ (or $C$), it is not possible to determine $c$ at all. Note that $m(k)$ also equals $[c(k)e^{j\theta_k}][c(k)e^{j\theta_k}]^* = c(k)c^*(k) e^{j\theta_k}e^{-j\theta_k}$ for arbitrary choice of $\theta_k$, and the inverse Fourier transform of $c(k)e^{j\theta_k}$ is also a real-valued signal if we are careful to preserve conjugate symmetry so that we choose $\theta_{-k}$ as $-\theta_k$. In other words, there are infinitely many quite different signals that share the same autocorrelation function, and it is not the case that only $C$ and its time-delayed versions have autocorrelation function $M$. As an example, all PN signals, that is, pulse trains generated from a maximal-length linear feedback shift register (LFSR), have an inverted thumbtack autocorrelation, and the distinct PN signals (from different LFSRs) are even two-level signals, a property not shared by all of the infinitely many real signals that have inverted thumbtack autocorrelation.

Answer 2

This is impossible in general as stated in the accepted answer.

The conditions under which it becomes possible to reconstruct the original signal given its autocorrelation need not be very restrictive however. If the signal that is to be reconstructed is sparse in some basis so that the constraint-ratio is larger than one certain phase retrieval algorithms can be employed to reconstruct the signal except for certain symmetries (shift, inversion).

In the simplest case assume the signal has extent a which is less than 1/2 the extent of the autocorrellation function in each direction. The simplest reconstruction algorithm is an iterative projection between the Fourier amplitude constraint and the known extent. This approach is known as error reduction. Starting from random phases and the amplitudes given by the Fourier transform of the autocorrelation signal, this signal is inversely Fourier-transformed. Then everything outside the known extent of the signal, the support, is imposed, everything outside is set to zero. This approximation is then Fourier-transformed and the amplitudes of the signal are set to the known amplitudes of the autocorrelation signal. This algorithm will eventually converge and produce an estimate of the signal with some probability greater than 0 , so usually it is run several times to generate a set of possible solutions from which the best one is chosen.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

n    = 256
mask = np.arange(n)<16
orig = np.random.rand()*mask
data = abs(np.fft.fft(orig))
recon= np.random.rand(n)
for i in range(4096):
  recon =np.fft.fft(b)
  recon =np.where(abs(recon)<1e-60,1e-60,recon)
  recon*=data/abs(recon)
  recon =mask*np.fft.ifft(recon)
plt.plot(orig)
plt.plot(recon)
plt.show()