Get a N-FFT with two N/2-FFT already computed

Question

After somme researchs on the web, I don't find the answer of my problem (or I don't understand it) and I hope this post will succeed. I'm working on a real-time FFT convolution algorithm (C++) which split the impulse response in several increasing size blocks (64/128/256/...). To realise the convolution, I have to split the incoming signal in blocks of the same sizes (64/128/...) to compute its FFts (when there is enough input signal to do this).

My idea is to use the FFts already computed to save some CPU time : with 2 FFts of 64 samples, I want to simulate the result of a FFT of 128 samples, ...

During the 1st step, we usually compute the FFT (128 samples (64 input samples and 64 zeropadding samples)

During the 2nd step, we usually compute the FFT (128 samples), and we join the 2 128-FFTs to create the 256-FFT

During the 3rd step, we usually compute the FFT (128 samples)

During the 4th step, we usually compute the FFT (128 samples), and we join the 2 128-FFTs to create the 256-FFT, and we join the 2 256-FFTs to create the 512-FFT

... and so on.

Little example:

1) x1 = {a, b, c, d, 0, 0, 0, 0}

X1 = fft(x1)

2) x2 = {e, f, g, h, 0, 0, 0, 0}

X2 = fft(x2)

   X  = X1+X2

X might be equal to : fft({a, b, c, d, e, f, g, h, 0, 0, 0, 0, 0, 0, 0, 0}

Is it even possible???

In my algorithm I use the AudioFFT library (https://github.com/HiFi-LoFi/AudioFFT) (Accelerate for Mac, FFTW for Windows and OouraFFT for Linux). Its class return an array for the real part and an array for the imaginary part (the size of each array is N/2+1, where N is the FFT lenght).

For now, I succeeded to join the result of the two smalls FFT (FFT1, FFT2) to obtain the even bins of the large FFT (FFT0) :

                    FFT0[0] = FFT1[0] + FFT2[0]

                    FFT0[2] = FFT1[1] - FFT2[1]

                    FFT0[4] = FFT1[2] + FFT2[2]

                    FFT0[6] = FFT1[3] - FFT2[3]

                    ...

But, how can I get the odd bins??????

Thank you in advance and I hope I was clear enough (despite my bad English).

Please, ask me for more details if necessary.

Thanks.

---

Thanks for your explanations. I thought I have the answer with that but finally no :'( The problem is the original signal is zeropaded. So, in your explanation user14819, the equation will look like : fft[a,b,c,d,e,f,g,h,0,0,0,0,0,0,0,0] = fft[a,b,c,d,0,0,0,0] + fft[0,0,0,0,e,f,g,h]

-> x1 (a,b,c,d,0,0,0,0) is actually the concatenation of the first and the third quarter of x(a,b,c,d,e,f,g,h,0,0,0,0,0,0,0,0) (or the zeropaded first half of the original signal)

->x2 (0,0,0,0,e,f,g,h) is the concatenation of the fourth and the second quarter of x

To show that, I used your Matlab code, Matt L., and I change it a little (with difficulties I admit):

N = 16;
N2 = N/2;
N4 = N/4;
n = (0:N2-1)';
g = e.^(-j*2*pi/N*n);
G = fft(g);
x = rand(N2,1); %. . . . . . . . . . . . original signal

x0 = x; %. . . . . . . . . . . . . . . . . . .zeropaded signal
x0(end+(N2)) = 0;

%your algorithm
x1 = x0(1:N2);%. . . . . . . . . . . . . first part of the zeropaded signal
x2 = x0(N2+1:N); %. . . . . . . . . . second part of the zeropaded signal

X1 = fft(x1);
X2 = fft(x2);

F1 = fft(ifft(X1).*g); %. . . . . . . . . or ifft(fft(X1).*fft(G))/N2
F2 = fft(ifft(X2).*g); %. . . . . . . . . or ifft(fft(X2).*fft(G))/N2

Xe = X1 + X2;
Xo = F1 - F2; %. . . . . . . . . . . . . .or F1 (because x2 is a zero array)
X = [Xe.'; Xo.'];
X = X(:);

max(abs(fft(x0)-X)) % near of epsilon

%my problem is x1 and x2 are the two zeropaded halfs of x (and not the two halfs of x0)
x1 = x(1:N4); %. . . . . . . .. . . . . . zeropaded first part of the original signal
x1(end+N4) = 0;

x2 = x(N4+1:N2); %. . . . . . . . . . zeropaded second part of the original signal
x2(end+N4) = 0;

X1 = fft(x1);
X2 = fft(x2);

F1 = fft(ifft(X1).*g);
F2 = fft(ifft(X2).*g);

Xe = X1 + X2;
Xo = F1 - F2;
X = [Xe.'; Xo.'];
X = X(:);

max(abs(fft(x0)-X)) %. . . . . . . . random because of the signal but too high

Do I make myself clear?

Answer 1

You can compute a length $N$ FFT from two length $N/2$ FFTs of consecutive parts of the data, but it can't be done efficiently, i.e. it is more efficient to go back to the time domain, combine the data, and compute the total FFT.

But anyway, the way it can be done is as follows. For convenience, I repeat the two formulas below from this answer:

$$X[2k]=\sum_{n=0}^{N/2-1}\left(x[n]+x[n+N/2]\right)W_{N/2}^{nk}\\ X[2k+1]=\sum_{n=0}^{N/2-1}\left(x[n]-x[n+N/2]\right)W_N^nW_{N/2}^{nk}\tag{1}$$

where $X[k]$ is the DFT of the complete data block, and where I assume that the length of the data block $N$ is even. If we denote the two halves of $x[n]$ by $x_1[n]$ and by $x_2[n]$ i.e.

$$\begin{align}x_1[n]&=x[n],&n=0,1,\ldots,N/2-1\\ x_2[n]&=x[n+N/2],&n=0,1,\ldots,N/2-1\end{align}$$

and if $X_1[k]$ and $X_2[k]$ are the length $N/2$ DFTs of $x_1[n]$ and $x_2[n]$, respectively, then we can rewrite (1) as

$$\begin{align}X[2k]&=X_1[k]+X_2[k]\\ X[2k+1]&=\text{DFT}_{N/2}\{W_N^nx_1[n]\}[k]-\text{DFT}_{N/2}\{W_N^nx_2[n]\}[k]\end{align}\tag{2}$$

If we let

$$G[k]=\text{DFT}_{N/2}\{W_N^n\}[k]\tag{3}$$

and remembering the DFT property

$$\text{DFT}_{N/2}\{u[n]v[n]\}=\frac{2}{N}(U\star V)[k]$$

where $\star$ denotes circular convolution, we can write the equation for the odd indices in (2) as

$$X[2k+1]=\frac{2}{N}\left((G\star X_1)[k]-(G\star X_2)[k]\right)\tag{4}$$

In sum, according to Eq. (2), the even indices of the DFT of the complete data block are simply given by the sum of the two shorter DFTs. The odd indices can be computed from (4), where circular convolution with the function $G[k]$ is necessary. This circular convolution performs interpolation in the frequency domain. From its definition (3), the function $G[k]$ is given by

$$G[k]=\frac{e^{j\pi (2k+1)/N}}{j\sin(\pi (2k+1)/N)},\quad k=0,1,\ldots,N/2-1$$

This short Matlab/Octave script shows how (and that) it works:

N = 128;            % (total) FFT length (even)
N2 = N/2;           % length of shorter FFTs
x = rand(N,1);      % some signal ...
x1 = x(1:N2);       % ... and its two parts
x2 = x(N2+1:N);
X1 = fft(x1);
X2 = fft(x2);
n = (0:N2-1)';
g = e.^(-j*2*pi/N*n);
G = fft(g);
F1 = ifft(fft(X1).*fft(G))/N2;   % circular convolution of X1 and G
F2 = ifft(fft(X2).*fft(G))/N2;   % circular convolution of X2 and G

% final result
Xe = X1 + X2;       % even indices of X
Xo = F1 - F2;       % odd indices of X
X = [Xe.'; Xo.'];
X = X(:);

% compare with direct computation of FFT
max(abs(fft(x)-X))  % 4.6204e-15

Answer 2

The FFT and DFT are linear operations ( fft[x+y] = fft[x] + fft[y] ).

Using an 8 point example:

fft[a,b,c,d,e,f,g,h] = fft[a,b,c,d,0,0,0,0] + fft[0,0,0,0,e,f,g,h]

In the above, you'll be replacing one 8 point fft with two zero padded eight point ffts, but since half of the 8 point inputs to the ffts are zero, you can prune the graphs to reduce the computations required. And if the input is real, you can save more. But you will still do slightly more computation than a single fft.

There are 16 basic types of fft: 8 DIT and 8 DIF, characterized by input order (bit reversed or sequential), output order, and geometry (in-place, isogeometric, constant input, constant output). There are a great many more non-basic FFT algorithms (mixed radix, sliding fft, special algorithms for hardware processors, etc., etc.). ALL of them require inputs to be staggered (ie: in all sixteen 8 point examples, input point 0 is ALWAYS added/subtracted to input 4, as are inputs 1 and 5, etc., no matter what the decimation, geometry, or input/output order).

You might indeed save some time if you can process the first half your data in less time than it takes for the second half to come in, but in some very high speed systems, that is simply not possible. And you will still be doing more operations than a single fft (and in inefficient algorithms/designs, much more).

Answer 3

The FFT butterflys depend on the accumulated state derived from all of the inputs for the last calculation. Thus all of the N must be FFT'd to get the last iteration. Two N/2 FFTs lose all of that correlated information.