Understanding DFT-ODD operation in Gardner Efficient Convolution paper

Question

As a newbie to DSP I am trying to understand and implement an efficient convolution engine as per JAES V43.3 1995/03 - Gardner - Efficient Convolution without I/O Delay.

I'm pretty much there with a working proof of concept, but I am struggling to get my head around appendix two - which talks about DFT-ODD operations and computing spectra using the spectra of two existing blocks of results.

Can someone point me to some additional (readable for newbie) documentations that perhaps covers this topic. This is the last piece in my jigsaw puzzle.

Answer 1

I'm actually in a discussion with Bill Gardner about this at this very moment (or month). I gave him a pdf document that outlined my thoughts. I wish I could effortlessly translate this to $\LaTeX$.

Well nothing is effortless. Let's remember that Bill's algorithm is about doing low-latency convolution reverb, using fast-convolution (like overlap-save or overlap-add) in real time. This means that there are definite relationships and limits between buffer size and delay for the various FIR segments of the convolver.

Convolution output is

$$\begin{align} y[n] &= h[n] \circledast x[n] \\ \\ &= \sum\limits_{i=-\infty}^{+\infty} h[i] \, x[n-i] \\ \end{align}$$

Where the finite impulse response is of length $L$ and is delayed by $D_0$ samples

$$ h[n] = 0 \quad \quad \text{for } n < D_0 \text{ or } n \ge D_0 + L $$

and is partitioned into $M$ adjacent and non-overlapping segments, each segment having length $L_m$.

$$ h[n] = \sum\limits_{m=0}^{M-1} h_m[n−D_m] $$

where

$$ h_m[n − D_m] \triangleq \begin{cases} h[n] \quad & D_m \le n < D_{m+1} \\ 0 & \text{otherwise} \end{cases} $$

and

$$ D_{m+1} = D_m + L_m \qquad \text{for } 0 \le m < M \qquad \qquad \qquad \qquad \textsf{(Eq 1)}$$

$D_0$ is the pre-delay of the entire FIR, the total length of the FIR (including pre-delay) is $D_M = D_0 + L$ and the nonzero length of the entire FIR is the sum of the lengths, $L_m$, of all $M$ segments:

$$ L = \sum\limits_{m=0}^{M-1} L_m $$

Each segment FIR, $h_m[n]$, is undelayed, zero-padded (after sample $L_m$), and transformed by the DFT to be the transfer function or frequency response, $H_m[k]$, for that $m$-th segment.

$$ h_m[n] \triangleq \begin{cases} h[n+D_m] \quad & 0 \le n < L_m \\ 0 & \text{otherwise} \end{cases} $$

$$\begin{align} H_m[k] &= \sum\limits_{n=0}^{N_m-1} h_m[n] \ e^{-j2\pi nk/N_m} \\ &= \sum\limits_{n=0}^{L_m-1} h_m[n] \ e^{-j2\pi nk/N_m} \\ \end{align}$$

$L_m$ is the length of the $m$-th segment. $D_m$ is the left boundary of the $m$-th segment and is the delay of that segment. $D_{m+1}-1$ is the right boundary of the $m$-th segment. The time-advanced impulse response of the $m$-th segment is

$$ h_m[n] \triangleq \begin{cases} h[n+D_m] \quad & 0 \le n < L_m \\ 0 & \text{otherwise} \end{cases} $$

Then the final output of the entire composite filter, $y[n]$, is the sum of the $M$ FIR filters:

$$\begin{align} y[n] &= h[n] \circledast x[n] \\ \\ &= \sum\limits_{i=-\infty}^{+\infty} h[i] \, x[n-i] \\ \\ &= \sum\limits_{i=D_0}^{D_M-1} h[i] \, x[n-i] \\ \\ &= \sum\limits_{m=0}^{M-1} \quad \sum\limits_{i=D_m}^{D_{m+1}-1} h[i] \, x[n-i] \\ \\ &= \sum\limits_{m=0}^{M-1} \quad \sum\limits_{i=D_m}^{D_m+L_m-1} h_m[i-D_m] \, x[n-i] \\ \\ &= \sum\limits_{m=0}^{M-1} \quad \sum\limits_{i=0}^{L_m-1} h_m[i] \, x[n-D_m-i] \\ \\ &= \sum\limits_{m=0}^{M-1} \quad h_m[n] \circledast x[n-D_m] \\ \\ &= \sum\limits_{m=0}^{M-1} \quad y_m[n-D_m] \\ \end{align}$$

where each segment FIR output is

$$\begin{align} y_m[n] &\triangleq h_m[n] \circledast x[n] \\ \\ &= \sum\limits_{i=0}^{L_m-1} h_m[i] \, x[n-i] \\ \end{align}$$

and is shown not yet delayed by the delay, $D_m$, for that segment.

The circular FFT convolution will convolve $h_m[n]$ with the undelayed $x[n]$ and the compensating delay of $D_m$ will occur as a consequence of processing in a double-buffered I/O context. Double buffering is necessary so that the thread for the $m$-th segment FIR can be guaranteed to have all $B_m$ samples of input $x[n]$ valid at the beginning of the thread and need not guarantee correctness of the $B_m$ output samples of $y_m[n−D_m]$ until the end of the thread (before looping and repeating).

The consequential minimum delay of double buffering the $m$-th segment FIR is at least twice the buffer or block size

$$ D_m \ \ge \ 2 \, B_m \qquad \qquad \qquad \qquad \qquad \quad \textsf{(Eq 2)}$$

Using the FFT to perform circular convolution, the length of the result of linear convolution of the FIR segment of length $L_m$ with sample block of length $B_m$ is the sum of the lengths less one. This must be no longer than the FFT size, $N_m$, which is also the period of the circular convolution result.

$$ N_m \ \ge \ B_m + L_m − 1 \qquad \qquad \qquad \qquad \textsf{(Eq 3)}$$

where

$$\begin{align} H_m[k] &= \sum\limits_{n=0}^{N_m-1} h_m[n] \ e^{-j2\pi nk/N_m} \\ \\ X_m[k] &= \sum\limits_{n=-L_m+1}^{N_m-L_m} x[n] \ e^{-j2\pi nk/N_m} \qquad \qquad \qquad \textsf{(overlap-save)} \\ \\ Y_m[k] &= H_m[k] \cdot X_m[k] \\ \\ y_m[n] &= \frac{1}{N_m} \sum\limits_{k=0}^{N_m-1} Y_m[k] \ e^{+j2\pi nk/N_m} \\ &= \frac{1}{N_m} \sum\limits_{k=0}^{N_m-1} H_m[k] X_m[k] \ e^{+j2\pi nk/N_m} \qquad 0 \le n < B_m \qquad \textsf{(overlap-save)} \\ \end{align}$$

and finally, we add up the outputs of all of the $M$ FIR segments

$$ y[n] = \sum\limits_{m=0}^{M-1} y_m[n-D_m] $$

Now here we must note that the delay $D_m$ of the $m$-th segment is what happens naturally or as a direct consequence of the double-buffering mentioned above. If equality in Eq 2 is taken (no extra delay added) the length of this delay is precisely twice the length of the buffer or block of samples, $B_m$. I think it's best to use the equality in Eq 2 and the inequality ($\ge$) in Eq 3 as illustrated in the example below.

From Eqs. 1-3 (taking equality),

$$\begin{align} D_{m+1} &= D_m + L_m \\ &= D_m + (N_m − B_m + 1) \\ &= D_m + (N_m − \tfrac{1}{2}D_m + 1) \\ &= \tfrac{1}{2}D_m + N_m + 1 \\ \end{align}$$

Suppose we were to double the buffer size (which also doubles the delay) for each succeeding $m$-th segment, then

$$ D_{m+1} = \tfrac{1}{2}D_m + N_m + 1 = 2 \, D_m $$

which leads to

$$\begin{align} D_m &= \tfrac{2}{3}(N_m + 1) \\ B_m &= \tfrac{1}{3}(N_m + 1) \\ L_m &= N_m - B_m + 1 = \tfrac{2}{3}(N_m + 1) \\ \end{align} $$

This results in a rather small buffer size, about a third of the FFT size. Normally fast convolution is more efficient with a buffer size significantly larger than the FIR length.

Suppose, instead, we double the buffer size and FFT size in the segment after the immediately following adjacent segments:

$$\begin{align} D_{m+1} &= \tfrac{1}{2}D_m + N_m + 1 \\ D_{m+2} &= \tfrac{1}{2}D_{m+1} + N_{m+1} + 1 \\ &= 2 \, D_m \\ \end{align} $$

However, for compact FFT code, we will only use an FFT of a size that is a power of 2

$$ N_m = 2^{p_m} \qquad \qquad \text{for } p_m \in \mathbb{Z} $$

And, say for even $m$, the FFT size for segment $m$ and segment $m+1$ are the same

$$ N_{m+1} = N_m = N_0 \cdot 2^{m/2} $$

For any integer $m$,

$$ N_m = N_0 \cdot 2^{\lfloor m/2 \rfloor} $$

where $\lfloor \cdot \rfloor$ is the floor(⋅) operator. Then, for even $m$,

$$\begin{align} D_{m+2} &= \tfrac{1}{2}D_{m+1} + N_{m+1} + 1 \\ &= \tfrac{1}{2}(\tfrac{1}{2}D_m + N_m + 1) + N_{m+1} + 1 \\ &= \ (\tfrac{1}{4}D_m + \tfrac{1}{2}N_m + \tfrac{1}{2}) + N_m + 1 \\ &= 2 \, D_m \\ \end{align} $$

This leads to

$$\begin{align} D_m &= \tfrac{6}{7} (N_m + 1) \\ B_m &= \tfrac{3}{7} (N_m + 1) \\ L_m &= \tfrac{4}{7} (N_m + 1) \\ \\ D_{m+1} &= \tfrac{10}{7} (N_m + 1) \\ B_{m+1} &= \tfrac{5}{7} (N_m + 1) \\ L_{m+1} &= \tfrac{2}{7} (N_m + 1) \\ \end{align} $$

Instead of increasing the delay by a factor of $\sqrt{2}$ each segment, the even-numbered segment increases the delay by a factor of $\tfrac{5}{3} > \sqrt{2}$ and the odd-numbered segments increase the delay by a factor of $\tfrac{6}{5} < \sqrt{2}$. Each combined adjacent segment pair doubles the delay, the buffer size, the FFT size, and FIR segment lengths. An interesting note is that the odd-numbered segment FIR length, $L_{m+1}$, is shorter than the FIR length, $L_m$ , of the preceding segment.

Of course these FIR lengths and buffer lengths and delays are integer values, so the floor() function must be applied to something. It's a little tricky, since Eq 1 is an equality and we must use the inequalities of Eq 2 and Eq 3 to absorb the rounding.

Here's an example: Suppose your first two FIR segments use an FFT of length 256. Then $N_0 = N_1 = 256$.

$$\begin{align} B_0 &= \big\lfloor \tfrac{3}{7} (N_0 + 1) \big\rfloor = 110 \\ D_0 &= 2 \cdot B_0 = 220 \\ L_0 &= 146 \le N_0 - B_0 + 1 \\ \\ D_1 &= D_0 + L_0 = 366 \\ B_1 &= \tfrac{1}{2}D_1 = 183 \\ L_1 &= 74 \le N_1 - B_1 + 1 \\ \end{align} $$

After that, for each adjacent pair of FIR segments, double the numbers for the pair just preceding.

$$\begin{align} N_{m+2} &= 2 \cdot N_m \\ D_{m+2} &= 2 \cdot D_m \\ B_{m+2} &= 2 \cdot B_m \\ L_{m+2} &= 2 \cdot L_m \\ \end{align} $$

The reason I chose $L_0=146$ instead of $L_0=147$ (which would be the longest $L_0$ allowed), is that I needed $L_0$ to be even, so when it's added to $D_0$, which is even (being twice $B_0$), that $D_1=D_0+L_0$ is even so that $B_1=\tfrac{1}{2}D_1$ is an integer. It all works out nice. $D_2=D_1+L_1$ is exactly twice $D_0$ as it should be.

It might be even more efficient to double the FFT length after three adjacent FIR segments. I.e. $ N_{m+3} = 2 \, N_m $. But this is getting complicated enough and my salient point is that the three numbered equations (Eqs 1, 2, 3) are fundamental and must be satisfied. When $B_m = L_m = \tfrac{1}{2} N_m$, then I am not confident that they are satisfied, which is why I am talking with Bill about this.

In another post, I show my reasons for saying that the sample buffer or block length $B_m$ should be the greater portion of the FFT buffer length $N_m$ for purposes of efficiency in fast convolution. If you're gonna really implement this in a real-time context, you have to worry about the details presented in this. That's another reason I have been discussing this with Bill.

Answer 2

The task is to compute the DFT of a sequence from the DFTs of its two halves:

$$X[k]=\text{DFT}_N\{x[n]\}=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$ $$\begin{align}x_1[n]&=x[n],\quad n=0,\ldots,N/2-1\\x_2[n]&=x[n],\quad n=N/2,\ldots,N-1\end{align}$$ $$X_1[k]=\text{DFT}_{N/2}\{x_1[n]\},\quad X_2[k]=\text{DFT}_{N/2}\{x_2[n]\}$$

i.e., compute $X[k]$ from $X_1[k]$ and $X_2[k]$. This exact question is treated in this answer, where it is shown (cf. Eq.$(2)$) that the even indices of $X[k]$ are easily computed as

$$X[2k]=X_1[k]+X_2[k]\tag{1}$$

but where it is also shown that there is no efficient way of computing the odd indices $X[2k+1]$ from $X_1[k]$ and $X_2[k]$, even though it is of course possible, albeit inefficiently.

For this reason, in his paper, Gardner suggests another method for computing $X[2k+1]$, called "DFT-ODD", which evaluates the DFT at odd indices. First, we need two equivalent representations of $X[2k+1]$:

$$\begin{align}1.\;X[2k+1]&=\text{DFT-ODD}_N\{x[n]\}=\text{DFT}_{N/2}\{(x_1[n]-x_2[n])W_N^n\}\\ 2.\;X[2k+1]&=\text{DFT-ODD}_{N/2}\{x[2n]\}+W_N^{2k+1}\text{DFT-ODD}_{N/2}\{x[2n+1]\}\end{align}\tag{2}$$

where $\text{DFT-ODD}_N\{x[n]\}$ denotes the DFT of the length $N$ signal $x[n]$ evaluated at odd indices (i.e., its length is $N/2$), and $W_N=e^{-j2\pi/N}$. The first representation is identical to the one given in the above mentioned answer, and the second one follows immediately from the definition of the DFT, after splitting the signal in its even and odd samples.

We need one more expression relating the DFT-ODD of a complex-valued sequence of length $N/2$ to the DFT-ODDs of its real and imaginary parts, respectively:

$$\begin{align}g[n]&=g_R[n]+jg_I[n]\\ \text{DFT-ODD}_{N/2}\{g[n]\}&=G[2k+1]=G'[k]\\ \text{DFT-ODD}_{N/2}\{g_R[n]\}&=\frac12\left(G'[k]+G'^*[N/4-k-1]\right)\\ \text{DFT-ODD}_{N/2}\{g_I[n]\}&=\frac{1}{2j}\left(G'[k]-G'^*[N/4-k-1]\right)\end{align}\tag{3}$$

where $*$ denotes complex conjugation.

The DFT-ODD of the original sequence $x[n]$ is now computed as follows:

1. generate a complex sequence from the even and odd indices of $x[n]$: $g[n]=x[2n]+jx[2n+1]$
2. compute DFT-ODD of $g[n]$ using representation $1$ of Eq. $(2)$: $G'[k]=\text{DFT}_{N/4}\{(g[n]-g[n+N/4])W_{N/2}^n\}$
3. compute the DFT-ODDs of $x[2n]$ and $x[2n+1]$ from $G'[k]$ according to Eq. $(3)$: $\text{DFT_ODD}_{N/2}\{x[2n]\}=\frac12\left(G'[k]+G'^*[N/4-k-1]\right)\\\text{DFT_ODD}_{N/2}\{x[2n+1]\}=\frac{1}{2j}\left(G'[k]-G'^*[N/4-k-1]\right)$
4. compute $X[2k+1]$ from these two DFT-ODDs according to the second representation in Eq. $(2)$.

This little Matlab/Octave script shows how it works:

N = 256;    % some multiple of 4
x = rand(N,1);
N2 = N/2; N4 = N/4;
x1 = x(1:N2); x2 = x(N2+1:N);
X1 = fft(x1); X2 = fft(x2);
Xeven = X1 + X2;

xe = x(1:2:N); xo = x(2:2:N); g = xe + 1i*xo;
G = fft( ( g(1:N4) - g(N4+1:N2) ) .* exp( -1i*2*pi/N2*(0:N4-1)' ) );  % DFT_odd of g
Xe = .5 * ( G + conj( G(N4:-1:1) ) );                          % DFT_odd of x[2n]
Xo = -1i*.5 * ( G - conj( G(N4:-1:1) ) );                      % DFT_odd of x[2n+1]

Xodd = [Xe;Xe] + exp(-1i*2*pi/N*(1:2:N-1)') .* [Xo;Xo];

X = [Xeven.'; Xodd.']; X = X(:);

% check result
X2 = fft(x);
max(abs(X2-X))    % 2.8422e-14