Evaluating discrete spectral density at only a few frequencies

Question

I'm trying to obtain the spectral density at three particular frequencies for a computational chemistry problem that I'm working on (if you are curious, it has to do with the estimation of Nuclear Overhauser Effects from molecular simulations). Three is a small enough number, especially given the length of the signals being sampled, that it seemed helpful to use a pruned DFT to do this (a Goertzel-like algorithm). Now I'm wondering if I can get away with that, or am somehow thinking about the problem incorrectly. I'm not a signal processing expert by any stretch of the imagination, so I am worried about missing something obvious.

In the chemistry literature, these effects are always functions of the spectral density of the autocorrelation function of some order parameter (in this case, if it matters, the dipole interaction tensor between an internuclear vector and the laboratory magnetic field, which I'm just going to be treating as 'the signal' below). From wikipedia's article on autocorrelation, for a signal $X(t)$: $$ \begin{align} F_R(f) &= \mathrm{DFT}[X(t)] \\ S(f) &= F_R(f)F_R^*(f) \\ R(\tau) &= \mathrm{IDFT}[S(f)] \end{align} $$

Since the spectral densities are given as the $\mathrm{DFT}[R(\tau)]$ in the textbooks/papers I've found discussing this physical effect, and since I only need the spectral density at three frequencies (0, 600, and 1200 MHz, if that matters) I thought the most straightforward solution was to use a Goertzel-like algorithm (this especially because much computation and disk IO not related to the DFT needs to be done to obtain each sample, so the single pass character of Goertzel-like algorithms is nice for this application).

$$ \begin{align} S(f) &= \mathrm{DFT}[R(\tau)] \\ F_R(f_0) &= \mathrm{Goertzel}_{f_0}[X(t)] \\ S(f_0) &= F_R(f_0)F_R^*(f_0) \end{align} $$

To reiterate the question, is the above reasoning valid? I've gotten funny results that have been hard to troubleshoot so I'm going back to questioning my underlying assumptions. It's also quite possible that all this works, but I'm applying it incorrectly, so if there are restrictions on which frequencies can be monitored in this way (or other things that would be obvious to someone in the field) that'd be helpful to hear.

The algorithm I'm using is sourced from Clay Turner's 'Oscillator Theory' approach to single frequency DFTs: here's C-like sample code he provided for it on a comp.dsp post from a while back:

// The input data is in x[], the data has N samples. And the bin number is w.

y1=0;
y2=0;
k=2*sin(pi*w/N);             // not 2 pi !!

for (j=0;j<N;j++) {
      y2=y2-k*y1+x[j];
      y1=y1+k*y2;
    }

// And the energy is simply

E = y1*y1 + y2*y2 - k*y1*y2;

Note that I chose this after doing some looking when I read Gentleman's 1969 paper suggesting that Goertzel's algorithm performs poorly for low frequencies (I need zero, so that seemed bad). I am working in C++, but because my 'signal' is matrix-valued my code does not look exactly like this. (Yes I can post it, but it might be TMI so I won't unless someone feels that might help).

Answer 1

I believe the OP's question is simplified to the following (confirming I didn't actually miss the salient question):

Given we can compute a power spectral density from the DFT using the conjugate product as:

$$S(k) = |X(k)|^2 = X(k)X^*(k)$$

Can we use a more efficient algorithm, such as the Goertzel, to compute a subset of $S(f)$ when only a few points are needed. The OP also wanted to understand if there was a limitation with using the Goertzel specifically, especially when only the lowest frequencies are of interest.

The answer is yes we can use a more efficient algorithm and specifically we can approximate DFT bins directly using the Goertzel. A windowing function is also recommended in either case (full DFT) or with the Goertzel to minimize aliasing issues (spectral leakage between bins).

The Goertzel is more efficient than the FFT when $M<log_2(N)$ where N is the total number of samples and $M$ is the total number of bins to compute. Other advantages are we can center each bin of interest at exact frequency values.

The Goertzel algorithm does have an error bound proportional to $N^2$ that is most pronounced near $\omega = 0$ and $\omega= 2\pi$. An alternate approach that would be identical to FFT results is to compute the DFT bins directly using the DFT formula, with the windowing suggested above also shown:

With windowing applicable to any approach: $x_{win}(n) = w(n)x(n)$, where $w(n)$ is a window, such as Kaiser, Hamming, Blackman, etc.

$$X(k) = \sum_{n=0}^{N-1}x_{win}(n)e^{-j2\pi nk/N} \tag{1}\label{1}$$

Which assuming a complex $x_w(n)$ requires a total of $4N$ real multiplies and $4N$ real adds for each bin computed. In comparison, the FFT would compute all $N$ bins with $2Nlog_2(N)$ real multiplies and $2Nlog_2(N)$ additions.

Therefore the total number of bins where it there is savings in computing directly using the DFT is given by $M$ for:

$$4NM < 2Nlog_2(N)$$

Resulting in:

$$M < \frac{log2(N)}{2}$$

So for example, if we had a 1024 point DFT, if less than $log_2(1024)/2 = 5$ bins were needed, it would be more efficient to compute them with the DFT equation directly rather than use the FFT. The Goertzel extends this threshold to 10 bins.

Further, as we can in the Goertzel, we can modify the DFT equation to center the bin on any frequency $\omega_o$ for $\omega_o \in [0, 2\pi)$ as follows:

$$X(\omega_o) = \sum_{n=0}^{N-1}x_w(n)e^{-j\omega_o n} \tag{2} \label{2}$$

In the DFT equation, $\omega$ is limited to discrete frequencies given by $k\omega_o$ where $\omega_o = 2\pi n/N$, thus it should be clear how $\ref{2}$ and $\ref{1}$ are related and how the form in $\ref{2}$ allows us to place the bin on any frequency along the continuous $\omega$ axis rather than the discrete axis for the DFT.

If windowing, please also refer to the last paragraph in this post to properly compensate noise and signal terms after using a window:

Blackman-Tukey Autopower equation

The FFTW page also provides a useful reference on Pruned FFTs and specifically an approach when the interest is in only the first K outputs which may be applicable to the OP's challenge:

http://www.fftw.org/pruned.html