1
$\begingroup$

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).

$\endgroup$
2
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for the fast response. I am not so worried about the instability at 0 or pi, because I'm using Clay Turner's Magic Circle that has different poles. The data I'm working with would be variable in length because there might be ten-thousands to millions of samples, but in any case only needing several frequencies does really favor approaches like this. I had already seen the FFTW page on pruning; While I only need several bins, they will not be the first several (they're not going to be at all close unless doing a very small DFT). I don't understand what is meant by bin centering, above. $\endgroup$ – LGS May 22 at 22:19
  • $\begingroup$ Also I'm marking this as answered because the main question was can I compute the spectral density in this way. I'm just trying to clarify (in the summary of my question provided by the responder) that I'm interested in other Goertzel like algorithms, not specifically Goertzel's algorithm. $\endgroup$ – LGS May 22 at 22:21
  • $\begingroup$ Sounds good. To answer your question on bin centering, in the DFT all bins are centered on $F_s/N$. And each bin acts as a bandpass filter for that frequency thus would be maximized if the actual frequency was aligned with the center of the band---so this applies when you are looking for a very specific frequency, with the approach above you can center any bin exactly on the frequency of interest. $\endgroup$ – Dan Boschen May 22 at 23:38
  • $\begingroup$ @LGS I also want to make sure it was clear to you that computing the DFT directly would be an alternative algorithm to using the Goertzel $\endgroup$ – Dan Boschen May 23 at 1:12
  • $\begingroup$ So literally how would that change the computation? Would I just use the freqency in radians per sample in the computation of k above? $\endgroup$ – LGS May 24 at 1:51

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.