# The autocorrelation vector of frequency data

I am trying to find the correlation of audio frequency data resulting from an FFT with itself. It would be interesting to hear what you all think on the following question :

Should one use the complex data from each bin to calculate autocorrelation through a sliding dot product, or should one use the real magnitude of each bin to do so.

My guess is that using the magnitude will only provide one with an array of coefficients corresponding to the similarity of magnitude levels in each bin of the FFT compared to its magnitude in the 0th bin.

What kind of autocorrelation would result using the complex data is more vague to me.

Thanks!

• I think I might not understand but my guess is: you want the dual to Parseval's theorem; i.e. convoluting a frequency function with itself corresponds to integration of pointwise(inner) multiplication of the original signal. You could try page 16 of princeton.edu/~cuff/ele301/files/lecture8_2.pdf and go backwards to see what he is saying. More generally Formula (21) in Laplace xforms in authors.library.caltech.edu/43489/1/Volume%201.pdf i.e. pdf page 131 Jul 19 at 20:15
• why do you want to autocorrelate your frequency-domain data in the first place? Jul 23 at 3:39

You would use the complex conjugate when correlating any complex data, such as an FFT result. Thus if a signal has a magnitude and phase, when correlated we multiply by the same magnitude and opposite phase-- this cancels out the phase in each product contribution the maximum to the overall correlation result (as a sum of products, or equivalently in the continuous time domain and integration of products).

Consider a simple example of correlating a 3 sample complex sequence:

$$x[n] = [e^{j\phi _1}, e^{j\phi _2}, e^{j\phi _3}]$$

Correlation is a sum of complex-conjugate products, specifically for an auto-correlation (not the auto-correlation function where we would repeat this at each time offset, but just the correlation at time index = 0):

$$\rho_{xx} = \sum_{n=0}^2 x[n]x^*[n]$$

We see here how due to the complex conjugate product, each product will equal 1 in this case and we would get the maximum correlation of 3:

$$\rho_{xx} = e^{j\phi _1}e^{-j\phi _1}+e^{j\phi _2}e^{-j\phi _2}+e^{j\phi _3}e^{-j\phi _3} = 3$$

If we did not correctly use a complex conjugate product, then for any phase values other than $$0$$ and $$\pi$$ (which means the sample is real and thus is already its own complex conjugate), the summation would be less than the maximum 3 that we achieved.

UPDATE Regarding Correlating with Just Magnitudes

Vatatia asked in the comments more specifically about using the product of the magnitude of complex waveforms (such as the FFT result) for the correlation rather than the complex conjugate products. This can be done and will provide a correlation result, but will have significantly more noise than operating on the complex values directly. This is consistent with the difference between "pre-detection" and "post-detection averaging, a great example is a spectrum analyzer for those that are familiar with that instrument: adjusting "Resolution Bandwidth" is adjusting the bandwidth of a filter prior to detection (taking the magnitude) and results in a decrease in the total power of the measured noise floor. In contrast, adjusting "Video Bandwidth" is adjusting the bandwidth of a filter post detection (filtering the output of a power detector) and only serves to smooth the measured noise but does not decrease the measured noise floor. The noise on the complex samples for many processes tends to be additive white Gaussian noise, which when zero mean will average toward zero, while this same noise once detected (magnitude of the complex signal) will be Rayleigh distributed which is NOT zero mean and will therefore accumulate in the average. Correlation for discrete time waveforms is a process of multiplying and adding (accumulating). As we add zero mean noise, sometimes it will add to the result and sometimes it will take away given it can be positive and negative; but if we take the magnitude first, it will always be positive and the result for the noise components will always go up. I can demonstrate this with a simple example in MATLAB of a sinusoidal tone in the presence of additive white Gaussian noise, and the resulting autocorrelation in frequency with and without taking the magnitude in the correlation.

wvfrm = randn(1, 2**16)+sin(2*pi*.01 .* [1:2**16])


Below is a plot of the test waveform; a real sine wave of amplitude 1 buried in real Gaussian noise with a variance of 1: Below if a plot of the magnitude of the FFT of this waveform; which is a complex result: The xcorr command will do a linear cross correlation (or auto-correlation when we only enter one waveform, or the same waveform twice), and as explained in the documentation for the function, it will take the complex conjugate product when processing. Below shows a plot of the result for doing the autocorrelation of the FFT result:

[autocorr, lag] = xcorr(fft(wvfrm))


Plot of magnitude of autocorrelation: Now observe what happens when we do the same operation, except with the magnitude of the fft instead of the complex result:

[autocorr2, lag] = xcorr(abs(fft(wvfrm))) What we see here is the significant accumulation of the noise floor when we used the magnitude of the waveform to compute the correlation rather than the complex result. (The floor is a triangular shape since this command used performs a linear correlation; if a circular correlation is done as described further below, the noise floor would be flat at the higher level in vicinity of the main correlation peak). Bottom line, using a magnitude to compute the correlation will lead to significant noise enhancement, even though the result will correctly show the same correlation peaks.

The Fourier Transform equation itself is a good example of the correlation of complex signals; for each frequency considered (represented as a single complex tone given as $$e^{j\omega t}$$) with a frequency $$\omega$$ radians/sec, we multiply our complex time domain signal ($$x(t)$$) by the complex conjugate for that frequency tone ($$e^{-j\omega t}$$), and integrate the product, resulting in the correlation at that specific frequency $$\omega$$:

$$X(\omega) = \int_{-\infty}^{\infty}x(t)e^{-j\omega t}$$

We also see the complex conjugate put to practice when we use the FFT to perform fast correlation; circular cross-correlation in the time domain is equivalent to the inverse FFT of the complex-conjugate product of the FFT's in the frequency domain:

$$\rho_{xy} = \text{IFFT}\{\text{FFT}\{y\}\text{FFT}^*\{x\}\}$$

Where $$\text{FFT}^*\{x\}$$ represents the complex conjugate of the FFT of x[n].

Similarly if we wanted to correlate FFT's directly (as a correlation function vs frequency offset index $$m$$) we could do this as:

$$\rho_{XY}[m] =\sum_{n=0}^{N-1}X^*[k]Y[k+m] \tag{1} \label{1}$$

Where:

$$X[k], Y[k]$$: are the FFT for $$x[n], y[n]$$

$$k$$ is the frequency index from $$0$$ to $$N-1$$

$$n$$ is the time index from $$0$$ to $$N-1$$

for $$N$$ samples in time and frequency.

$$\rho_{XY}[m]$$ is the correlation of the FFT results versus frequency offset index $$m$$.

Changing which variable gets indexed and which gets complex conjugated would change the interpretation of the shift in frequency for each correlation result. Importantly and not accurately represented in the simplified formulation above, $$Y(k+m)$$ is a circular shift of the sequence $$Y(k)$$.

Note that given the time frequency duality properties of the DFT (and Fourier Transform in general) we can instead compute the circular correlation in frequency in the time domain using the same formulation as \ref{1} but changing time for freq:

$$\rho_{XY} = \text{FFT}\{(y x^*\} \tag{2} \label{2}$$

To be clear on the use of the above, if $$x$$ and $$y$$ were similar such that they would be correlated if they were perfectly aligned in time and frequency, then in the case that there was no time misalignment but a frequency offset existed (such as carrier offset error between transmitter and receiver in a wireless comm link), the above would determine that frequency offset by providing a function of correlation versus frequency. If there was both a frequency and time offset (as typically the case), then the correlation would need to be done in both dimensions, as I have demonstrated with FFT's here.

As a simple demonstration of the above, consider this example with a 11-chip barker sequence in the time domain as BPSK symbols (+1 or -1):

$$x[n] = [1,-1, 1, 1, -1, 1, 1, 1, -1, -1, -1]$$

This sequence has a desirable time domain auto-correlation which is plotted below. When the sequence is aligned with itself in time, it's autocorrelation (as a sum of products) is maximum and equal to the length of the sequence (11). For all other values the magnitude will be 1 (since the sequence length is odd the sum of products can't be zero for any given shift). Given this property, it is useful for time domain resolution tasks such as synchronization and ranging. If there was a sufficient frequency offset in one of the sequences, then the correlation would be significantly degraded until the offset was corrected. For example, given a frequency index $$k$$ from 0 to 11 in the FFT (for the 11 sample sequence), I introduce a frequency offset of k_{offset} = 3.3 by a time domain product as follows:

$$y[n] = x[n] e^{j 2 \pi k_{offset} n/N}$$

A plot of the Discrete-Time Fourier Transform (DTFT) magnitude confirms the frequency translation that was introduced (note that I scaled the continuous frequency axis $$\omega$$ that extends from $$0$$ to $$2\pi$$ to show index $$k$$: Below is an updated time domain cross-correlation, where we see that the frequency offset has destroyed the nice time domain correlation properties we desired: Given in this set-up example we only have a frequency offset and not a time offset, we can use \ref{2} directly to determine the frequency (carrier) offset. Further, we can zero pad the FFT to resolve the fractional index frequency offset:

Below I show the result in Python using an FFT of the 11 sample product padded out to 256 samples:

rho_freq = fft.fft(y * np.conj(x), 256) (Note that we could arbitrarily conjugate y instead of x, as in rho_freq = fft.fft(x * np.conj(y), 256). This would simply change which waveform is the reference and we would get a similar result but with a circular negative frequency shift instead of the positive one indicated in the plot above (the peak would be at k=11-3.3, instead of +3.3). This is simply a sign change and we just need to be aware of which convention we are using when accessing the result to determine the frequency error and applying the correction.)

From the above we can sufficiently resolve the frequency error, correct it by multiplying our time domain signal with the complex conjugate, and recover the nice time domain correlation that we desire:

$$x\hat[n] = y[n] e^{-j 2 \pi k_{offset} n/N}$$

A utility of this is correcting for carrier offsets in a communication link where transmitter and receiver are not yet synchronized: by using such codes in an acquisition header, we can get our initial estimates for carrier offset (as I showed here) as well as timing offset, which can then be maintained with carrier and timing loops that track these offsets as they change by monitoring the rest of the message using sample by sample tracking techniques.

• A great explanation! I do want to ask you three things: 1. What resources are out there for learning complex correlation? 2. You mention in the three-element complex array example that the autocorrelation at time lag zero will be the dot product of the array -> so the autocorrelation at every lag value can be found using a sliding dot product -> my question is: will the resultant array always only represent the correlation of element 0 with the rest of the array? 3. I am guessing the FT equation is an example of correlation rather than suggesting correlating FFT data at each f using x(t)? Jul 21 at 19:07
• @Vatatia I teach three online courses that are all focused on providing the intuitive view of signal processing with complex signals. I may be biased but I think they are the best resources out there! You can find more info on my profile and watch these links for course dates: dsprelated.com/courses and ieeeboston.org/2022-courses The course that introduces correlation in most detail is "DSP for Wireless Communications" which will be offered the October through the IEEE link (not posted yet), or in January as posted on dsprelated.com. Jul 21 at 19:21
• Re your question- yes a dot product is a correlation. It will be the correlation of the first array with the second array (not just the first element). The slide represents the time shift for that correlation. The FT equation is an optimized correlation using a Fourier Transform property. The result is the same as doing a circular rotation and taking the dot product for each shift. See this post for further details: dsp.stackexchange.com/questions/22877/… Jul 21 at 19:25
• Interesting! I will surely keep your course in mind seeing as I have not come across many resources going into complex correlation in DSP. As a reply: I understand that the slide represents the time shift for that correlation, but would you argue that a localized autocorrelation could be more effective in the case of audio frequencies due to the trends that exist in audio data? I ask this because I am using normalized autocorrelation across all available bins, however, my resulting correlation vector seems to oscillate quite a bit. Jul 22 at 1:23
• Hmmm I am not really an expert in audio (my background in wireless Comm) but there are plenty of audio experts here so consider posting that as a question. Did I answer your original question as posted to your satisfaction ? Jul 22 at 1:30