0

I am calculating an autocorrelation function of a signal using Numpy's FFT by first calculating the power spectrum. I believe that autocorrelation functions are positive definite. However, when I test if my simulated autocorrelation function is positive definite, it often fails. Here's some example code:

import numpy as np

def is_pos_def(x):
    return np.all(np.linalg.eigvals(x) > 0)

n = 5
signal = np.random.random(size=(n,n)) #create a real valued signal
F = np.fft.fftn(signal)               #FFT of real signal is conjugate symmetric
A = np.abs(F)                         #amplitudes/magnitudes
ps = A**2                             #power spectrum is square of magnitudes
acf = np.fft.ifftn(ps)                #autocorrelation function is ifft of the power spectrum
acf = acf.real                        #.real to clean up since the acf is real anyways by definition, and indeed the imaginary terms are zero.

print(is_pos_def(acf))

This outputs False. Am I missing something?

tomerg
  • 353
  • 2
  • 12
  • I have no idea what the answer to your question is, but I do note that several of the eigenvalues you get is `is_pos_def` are complex. Is that expected behavior? – Frank Yellin Jan 20 '22 at 21:21
  • hmm, I did not get complex values when I ran it. However, running it a few Moree times I did occasionally get complex values, which I guess is just dependent on the random seed. – tomerg Jan 20 '22 at 22:01
  • “Autocorrelation function is positive definite”, but you link to a question about the autocorrelation matrix, which is not the same thing. “Positive definite” is a property of a matrix, not of a function. Note that if the input were a 1D function or a 3D function, then you wouldn’t be able to do your computations. – Cris Luengo Jan 20 '22 at 23:53
  • @CrisLuengo Oh I didn't realize they were different. I guess I mean matrix then, rather than function, since I'm working with numpy arrays (does it matter if it's a numpy array rather than a numpy matrix?). I will update the question to say matrix. – tomerg Jan 21 '22 at 00:00
  • No, your array is a discrete function, and it’s not an autocorrelation matrix. The autocorrelation matrix follows from a random vector, not from a signal or function. https://en.wikipedia.org/wiki/Autocorrelation#Auto-correlation_of_random_vectors — We’re talking different concepts here. There is no reason to expect the eigenvalues of your image matrix to be positive, or even to exist… – Cris Luengo Jan 21 '22 at 00:05
  • I'm confused. In the code above I am starting with a random array and calculating the autocorrelation as the fft of the power spectrum. You're saying that is not the same as an autocorrelation matrix? – tomerg Jan 21 '22 at 00:10
  • @CrisLuengo so I'm going to have to try and wrap my head around what you said. But to start, is this what you're getting at, in particular the section near the bottom of the page: https://www.gaussianwaves.com/2015/05/auto-correlation-matrix-in-matlab/ ? – tomerg Jan 21 '22 at 00:37

1 Answers1

0

Regarding positive definite functions

After some rounds in the comment I suspect that you want to check that the autocorrelation is a positive definite function

import scipy.linalg;
n = 5; # number of independent signals
m = 10; # number of samples each vector
signal = np.random.random(size=(m)) #create a real valued signal
# auto-correlation function
acf = np.fft.ifft(np.abs(np.fft.fft(signal))**2)
# check function spectrum of the matrix a[i,j] = f(x[i]-x[j])
spectrum = np.linalg.eigvals(scipy.linalg.toeplitz(acf))

The auto-correlation function positive definiteness is confirmed by

assert np.allclose(spectrum, spectrum.real) # real
assert np.all(spectrum.real > 0) # positive

The original answer

First, the auto correlation of a signal is ifft(abs(fft(x)**2), not ifft(fft(x)). From the convolution theorem we know that given A=fft(A), B=fft(b), ifft(A*B) = conv(a,b).

This autocorrelation will compute the correlation of a signal with different delayed versions of it.

The autocorrelation matrix on the other hand, computes the correlation between each pair of columns in a matrix.

n = 5; # number of independent signals
m = 10; # number of samples each vector
signal = np.random.random(size=(m,n)) #create a real valued signal
acf = signal.T.conj() @ signal;

Then acf is guaranteed to be a symmetric semidefinite positive matrix, i.e. all the eigenvectors are real and non-negative (>=0)

Bob
  • 13,867
  • 1
  • 5
  • 27
  • I do take ifft(abs(fft(x))**2) in my code, i.e. the power spectrum. However, the result of my code gives an acf whose eigenvalues are not positive, hence my question. Also, I'm confused about your last line. The acf is not just the conjugate transpose of the signal (I'm not sure what the @ signal is, is that just a typo meant to be # signal for a comment?). – tomerg Jan 26 '22 at 13:20
  • `@` is the [matrix multiplication operator](https://www.python.org/dev/peps/pep-0465/) – Bob Jan 26 '22 at 14:20
  • Notice that for complex values `A**2` is not the same as `abs(A)**2` – Bob Jan 26 '22 at 14:22
  • Yes, but I do A = np.abs(F) in the line above that – tomerg Jan 26 '22 at 15:07
  • Sorry, I missed that detail on my first read and I never read again hehe – Bob Jan 26 '22 at 15:08
  • no worries. Also I didn't know about the @ operator (other than as a decorator), I usually work with ndarrays and not matrices. Now I know! – tomerg Jan 26 '22 at 18:23
  • Did you get the point point is that the [auto correlation matrix](https://en.wikipedia.org/wiki/Autocorrelation#Matrix) is different from what you are calculating. – Bob Jan 26 '22 at 18:41
  • Yes, I now understand that the autocorrelation matrix is not the same as the autocorrelation function I am calculating, and indeed this page was helpful too: https://www.gaussianwaves.com/2015/05/auto-correlation-matrix-in-matlab/. However, should I still expect my autocorrelation function to be positive definite in my case? I ask because the fft(acf) result would obviously in a positive power spectrum (since I started with a positive PS), which requires the acf to be positive definite from what I've read on Bochner's theorem, and other posts on SO, so I'm confused why my acf is not pos def. – tomerg Jan 27 '22 at 20:50
  • Thank you for the update, which is quite helpful. Note that this does not appear to work when the signal is (n,n) shape as in the original question, rather than (n,) shape. Also, I'm not sure why you are constructing a toeplitz matrix first (which I had to look up)? I need the acf to be positive definite, not the toeplitz(acf), i.e. I need the fft(acf) to be positive and real. – tomerg Feb 02 '22 at 23:50
  • You wanted to test positive-definiteness of a function, that is done by testing the positive-definiteness of a matrix constructed with `a[i,j] = f(x[i] - x[j])`, if `x` is from a equally spaced grid, that matrix is a toeplitz matrix. – Bob Feb 03 '22 at 11:38