2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ This answer could be helpful. I think it's about the same problem, but it could be that Gardner goes a bit further than that. $\endgroup$ – Matt L. Jan 1 '17 at 16:55
  • $\begingroup$ Gardner appears to imply that while doing "efficient convolution" that this is one of the optimizations to employ - seeing as we already have spectral information for both halves of the previous block of signal. I'm making do without it for now but would be interested in trying to get it implemented eventually. $\endgroup$ – Mark Jan 2 '17 at 14:07
  • $\begingroup$ The trick is neat, but there are some arguments against using it. If the sizes of the partitions in an optimal partitioning are not successive powers of two, then the trick becomes less effective. Often FFT libs are available that are highly optimized for the platform. Those cannot be used with the trick. If you do the platform optimization yourself, the code becomes less portable. $\endgroup$ – Olli Niemitalo Jan 4 '17 at 17:44
2
$\begingroup$

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

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. I'm still mulling this.

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ For a newbie to DSP, Bills paper is remarkably readable - probably just a little bit over my head maths wise though - I need to improve my skills in that area. The two areas that are causing me a challenge at the moment are scaling of the DFT outputs and the DFT-ODD implementation. I am working on the scaling issue at the moment before I try and get my head around the DFT-ODD optimisation. The output of my code sounds quite smooth until I loop around one filter block when it starts to stutter a bit. Still not sure why yet. Might post a soundcloud for comment.... $\endgroup$ – Mark Jan 4 '17 at 1:10
  • $\begingroup$ my answer will spell it out when i get done with it. there are some subtle but necessary issues (regarding real time constraints) that i am not sure @MattL. is dealing with. and there are optimization issues about it that Bill hasn't dealt with, which have to do with the segment FIR length, the buffer length, the FFT length, and the segment delay. (dividing the FFT of these segments in half and progressively doubling the segment size does not work.) $\endgroup$ – robert bristow-johnson Jan 4 '17 at 1:40
  • $\begingroup$ BTW, this must be said: I have never seen a 4 month latency between "submission" (as a convention paper) to publication in the AES Journal. Bill's paper went from first public exposure to publication faster than any other paper in the modern history of the Journal. I think that is an indication of the creativity and novelty of thought regarding the paper. $\endgroup$ – robert bristow-johnson Jan 4 '17 at 1:45
2
$\begingroup$

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
| improve this answer | |
$\endgroup$

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.