What is the adjoint of the Goertzel algorithm?

Question

The adjoint of the DFT is the IDFT. But suppose we don't need to calculate a full DFT (just a subset of frequencies), and so are using the Goertzel algorithm instead. What is the adjoint of the Goertzel algorithm?

One way to think about this might be to model the Goertzel operation as a full FFT followed by truncation to the frequency indices of interest. Then the adjoint of that operation would be:

1. Zero pad back up to the full frequency vector (fill in the missing frequency bins with zero values)

2. Take the inverse FFT

This would give the adjoint operation, but calculating a full FFT and throwing away many of the values is not efficient. Can anyone provide insight on a more efficient adjoint Goertzel operator?

Edit

The definition of "adjoint" I'm using here is that of the Hermitian adjoint from linear operator theory. It is a generalization of the complex conjugate transpose of a matrix to linear operators.

More information on adjoints can be found at this answer.

Answer 1

My assumption is the OP wants to reconstruct the time domain samples corresponding to $M$ arbitrarily chosen DFT bins of an $N$ sample DFT.

As given by the inverse DFT formula directly, each complex coefficient $c_k$ in the frequency domain DFT result represents the magnitude and starting phase for a complex spinning phasor in the time domain given as:

$$x_k = c_k e^{j k \omega_o n}$$

Where $\omega_o$ is the fundamental frequency given as $f_s/N$ with $f_s$ as the sampling rate, $k$ is the frequency index (harmonic number) and $n$ is the time index extending from $0$ to $N-1$ for an $N$ sample DFT.

Thus if there are only a small number $M$ of frequency coefficients (sparse DFT) with a DFT of length $N$, the time domain can be recomputed directly with $MN$ complex multiplications. A full Inverse DFT in comparison requires $(N/2)log_2(N)$ complex multiplications when using the inverse FFT algorithm.

If we equate the two:

$$MN = (N/2)log2(N)$$ $$M = \frac{log2(N)}{2}$$

So for example, with a 1024 length DFT, $M=log2(1024)/2 = 5$, so if we have fewer than 5 frequency bins in this case, it would be more efficient to recompute the time domain waveform directly as:

$$x[n] = \sum_k c_k e^{jk\omega_o n}$$

Otherwise for greater than 5 bins, it would be more efficient to use the Inverse FFT algorithm. In comparison, for a 1024 length waveform, the Goertzel algorithm is more efficient than using the FFT directly when less than 10 frequency bins are needed (For more details of that case, see this post).

A concept I haven't tried as I just thought of it now promises to be as efficient as the Goertzel at the cost of quarter cycle memory storage and phase accumulators for each bin. This approach is similar to an NCO, except we maintain a separate phase accumulator for each frequency bin at a starting phase as given by that frequency bin. Each accumulator output (representing phase) is converted to complex magnitude with a quarter-cycle storage look-up table, each complex magnitude is scaled corresponding to the magnitude for that bin (requiring two real multipliers) and the scaled magnitudes are summed resulting in each time domain sample. Thus there are 2NM real multipliers, or half that given with the direct inverse DFT approach described above. Any of the techniques used for memory reduction in NCO's (Hutchinson's Algorithm, Sunderland's Algorithm, sample interpolation, CORDIC rotator etc would also apply here if further memory storage efficiency is needed).

This algorithm is diagrammed below. Real data paths are represented with thinner traces, and complex data paths (as $I+jQ$) with wider traces. The diagram shows the inverse DFT construction with three arbitrary frequency bins given as real indexes $k_1$, $k_2$ and $k_3$. Each bin has a complex weight given as $A_ke^{j\theta_k}$ with $A_k$ and $\theta_k$ as real numbers. The phase $\theta_k$ is the starting index within the fixed precision accumulator size corresponding to a range of $[0, 2\pi)$ with the accumulator set to wrap on overflow. For each output sample, the commutators go to each accumulator and use a single look up table to compute $e^{j\theta} = \cos(\theta)+j\sin(\theta)$. Note that this can be done with one real look-up table with a quarter cycle storage for the $\sin(\theta)$ result. Each complex result is then weighted by real $A_k$ and the results for all are summed for each time domain output sample $x[n]$. The accumulator then advances to the next value (by adding $k_n$ to the accumulator total) and the process is repeated.

I have also considered the idea of using guaranteed stable quadrature oscillators such as these described by Rick Lyons, using a separate oscillator for each bin, but they all require multipliers and the complex outputs would need to be scaled for each bin prior to combining the time domain result. The scaling alone is 2 real multipliers per time sample $n$, and thus would not be more efficient than the approach described above. (Note the structure of the quadrature oscillators is very similar to that of the transfer function of the Goertzel, essentially a complex pole on the unit circle with features to ensure rounding does not deviate the pole from stable oscillation).

Answer 2

Goertzel's Algorithm has more than one definition running around. I'm going to assume here that you're using the definition where -- excepting numerical errors -- it behaves exactly like one bin of the DFT.

1. Zero pad back up to the full frequency vector (fill in the missing frequency bins with zero values)

2. Take the inverse FFT

That will give you the least square error reconstruction of the original signal. So from your definition above, it will give you an adjoint operator on the Goertzel, on multiple parallel Goertzels of different frequencies, or on any other sparse vector of DFT results.

Note that depending on the signal content and your goals, this could be a pretty bad reconstruction. If what you need isn't well encompassed by a least square error reconstruction, there may be other things that you'd need to do.