0
$\begingroup$

I am looking over Blackman-Tukey method for Autopower.

It uses the DFT after a window is applied on a autocorrelation.

$$ Power Spectrum = \frac1{2\pi} \ \sum_{(k=-(N-1))}^{N-1} w[k] R[k] e^{-i\omega k} $$

  1. Where does the $ 1/(2\pi) $ come from?
    It's not from the Autocorrelation $R[k]$ or window $w[k]$ and it's not an inverse transform, it's forwards, so why require it here?
  2. $e^{-i\omega k}$ why is this not the DFT exponential term of $e^{-i 2\pi kn/N }$
    We seem to be missing an $n$?
  3. Is this autocorrelation method 'biased' by division of $N$ (length of raw original signal) after the autocorrelation calculation?
    I assume it must be to give a Power spectrum because if there was no division of N it would be giving Energy spectrum.
$\endgroup$
1
$\begingroup$

The equation looks reasonable to me: The scaling by $\frac{1}{2\pi}$ is to have the result in units of normalized frequency $f$ instead normalized angular frequency $\omega$. What may be confusing is using the index $k$ from the autocorreation instead of $n$ since it would be the time domain variable for $R[k]$ and $w[k]$, while $k$ is often associated with the frequency index in the DFT. The OP is seeking the frequency variable to be in the exponent, in this case suggesting $n$, but instead the unit used is actually $\omega$ representing the fundamental normalized angular frequency (which for a DFT over $N$ samples would be given as the expected $2\pi n/N$ as in the OP's 2nd question).

In this case, there are $2N-1$ total samples consistent with a linear correlation over $N$ samples apparently given as:

$$R[k] = \sum_{n=0}^{N-1}x[n]x^*[n-k], \space\space\space k \in [-N-1, \ldots N-1]\tag{1}\label{1}$$

The Discrete Fourier Transform of the windowed autocorrelation function could then be (using a frequency index $\ell$ and assuming there are also $2N-1$ frequency samples):

$$X(\ell) = \sum_{k=N-1}^{N-1}w[k]R[k]e^{-j2\pi k \ell /(2N-1)}, \space\space\space \ell \in [-N-1, \ldots N-1]\tag{2}\label{2}$$

This will grow by a factor of $2N-1$, which we could normalize if desired by dividing by $2N-1$. As in any DFT, this is just a scaling and its used would be dictated by defining the units of the output.

Finally the normalized radian frequency (the fundamental frequency when $\ell = 1$) in this case is $\omega = 2\pi \ell/(2N-1)$, if substitute this in $\ref{2}$ we get a function of $\omega$:

$$X(\omega) = \sum_{k=N-1}^{N-1}w[k]R[k]e^{-j\omega k}\tag{3}\label{3}$$

Or if instead if we want a function of $f$, given $\omega = 2\pi f$:

$$X(f) = \frac{1}{2\pi}\sum_{k=N-1}^{N-1}w[k]R[k]e^{-j\omega k}\tag{4}\label{4}$$

I was not able to locate an actual copy of Blackman and Tukey's 1958 paper (which predated the FFT) but did find this additional variation from "Noninvasive Instrumentation in Medical Diagnosis" by Robert B. Northrop. In this variant, I would think that either the exponent should be divided by $(N+1)$ or the summation and number of output frequencies should be done over $N$ samples:

Blackman Tukey Autopower

The OP asked in comments if this is best practice for estimating power spectral density. The power spectral density is given as the Fourier Transform of the auto-correlation function, and doing a linear correlation (which can be done by zero-padding the time domain sequence and doing a circular correlation) and windowing the resulting time domain sequence prior to computing a DFT are indeed best practice for estimating the Fourier Transform. The linear correlation avoids time domain aliasing, and the windowing reduces spectral leakage error. So this is standard practice.

One precaution not yet mentioned with windowing and spectral estimation: the equivalent noise bandwidth of a rectangular window is 1 DFT bin, which means when the spectrum in generally equally distributed the total power is equivalent to being distributed neatly in each DFT bin without double counting - and we can then use Parseval's Theorem to measure the total power (or the power across a band of frequencies) by summing the power in each bin. Windowing increases the equivalent noise bandwidth of each DFT bin (as given by the kernel, or discrete-time Fourier Transform, of the window), so that if we added up the power in each bin double counting would result and the total power would be over-estimated. The window itself also removes power from the signal since we are decreasing the signal over much of the time interval, resulting in a net loss (the loss from the window can be predicted based on the sum-square of the window samples, and the net loss is this loss plus the gain from the increase in kernel bandwidth). A single tone that is not spread out over multiple bins would only be effected by the power loss but not the gain from double counting, so the SNR is effected. This is called the processing gain of the window as detailed by fred harris in his classic paper "On The Use of Windowing in the Fourier Transform" and further detailed in this other posting:

How to calculate resolution of DFT with Hamming/Hann window?

$\endgroup$
5
  • 1
    $\begingroup$ That's very well explained. Thank you very much. A book I am reading recommended this is a 'safe' and best approach instead of performing a standard FFT. But you seem an expert and perhaps I am now thinking this method is not so common? $\endgroup$ – Natalie Johnson May 12 '20 at 11:43
  • $\begingroup$ The only difference that makes it "safe" and "best" is windowing the time data before taking the FFT. This is good and best practice. However what it not yet explained is that you then need to factor in the windowing loss and the effect of resolution bandwidth from the window. I will add details on this at the bottom of my answer as it is critical for doing actual spectral power measurements. $\endgroup$ – Dan Boschen May 12 '20 at 11:48
  • $\begingroup$ On second thought. The 1/2pi division on an fourier transform is from the fact we are assuming periodicity by letting the period go to infinity and it has to go on either the forwards or reverse. The DFT is divided by 1/N because of tthis. This is a transform, why is it not divided by 1/(2N-1) since now the length 2N-1 after autocorrelation? $\endgroup$ – Natalie Johnson May 12 '20 at 11:49
  • $\begingroup$ Nevermind, you address this in your answer but saying we can normalise by this if required. $\endgroup$ – Natalie Johnson May 12 '20 at 12:03
  • $\begingroup$ Right. I typically do scale my transforms by $M$ where $M$ is the total number of samples, just to put the results back to the same general range as my input: $Ae^{j\omega t}$ would map to a bin with magnitude $A$ and I personally like that. It is all just math and like I said as long as you are clear with what units the output represents, you can scale to your preference.) $\endgroup$ – Dan Boschen May 12 '20 at 12:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.