Sunday, 5 June 2016

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


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?




No comments:

Post a Comment

readings - Appending 内 to a company name is read ない or うち?

For example, if I say マイクロソフト内のパートナーシップは強いです, is the 内 here read as うち or ない? Answer 「内」 in the form: 「Proper Noun + 内」 is always read 「ない...