Forums

combining small CONSECUTIVE FFTs to make a larger FFT

Started by Intuitive June 16, 2010
On Thu, 17 Jun 2010 10:38:17 -0500, "Intuitive"
<andrea3891@n_o_s_p_a_m.gmail.com> wrote:

>I don't understand the approach proposed by Greg too. Can you give a better >explanation? (Possible step by step).
Use a fixed-width font. The DFT value in frequency bin "k" of the data set starting at sample "n0=0" is: 1023 X(k,n0) = SUM [x(n)*exp(-j2PIkn/1024)] n=0 The DFT in frequency bin "k" of the data set starting at sample "n1=128" is: 1151 X(k,n1) = SUM [x(n)*exp(-j2PIkn/1024)] n=128 1023 = SUM [x(n)*exp(-j2PIkn/1024)] n=0 1151 + SUM [x(n)*exp(-j2PIkn/1024)] n=1024 127 - SUM [x(n)*exp(-j2PIkn/1024)] n=0 1023 = SUM [x(n)*exp(-j2PIkn/1024)] n=0 127 + SUM [x(n+1024)*exp(-j2PIk(n+1024)/1024)] n=0 127 - SUM [x(n)*exp(-j2PIkn/1024)] n=0 1023 = SUM [x(n)*exp(-j2PIkn/1024)] n=0 127 + SUM [x(n+1024)*exp(-j2PIk(n+1024)/1024) - x(n)*exp(-j2PIkn/1024)] n=0 1023 = SUM [x(n)*exp(-j2PIkn/1024)] n=0 127 + SUM {[x(n+1024)-x(n)]*exp(-j2PIkn/1024)} n=0 127 = X(k,n0) + SUM {[x(n+1024)-x(n)]*exp(-j2PIkn/1024)} n=0 So the DFT at of the data set starting at sample 128 is equal to the DFT of the data set starting at sample 0 plus something resembling a 128-point DFT of the difference between the "newest" 128 samples and the "oldest" 128 samples. (Note the "1024" in the denominator, however.) Since you are building the DFT gradually, you can define the "oldest" 128 samples to be zeroes: 127 X(k,n1) = X(k,n0) + SUM {x(n+1024)*exp(-j2PIkn/1024)} n=0 Similarly, the DFT in frequency bin "k" of the data set starting at sample "n2=256" is: 1279 X(k,n2) = SUM [x(n)*exp(-j2PIkn/1024)] n=256 1151 = SUM [x(n)*exp(-j2PIkn/1024)] n=128 1279 + SUM [x(n)*exp(-j2PIkn/1024)] n=1152 255 - SUM [x(n)*exp(-j2PIkn/1024)] n=128 1151 = SUM [x(n)*exp(-j2PIkn/1024)] n=128 255 + SUM [x(n+1024)*exp(-j2PIk(n+1024)/1024)] n=128 255 - SUM [x(n)*exp(-j2PIkn/1024)] n=128 1151 = SUM [x(n)*exp(-j2PIkn/1024)] n=128 255 + SUM [x(n+1024)*exp(-j2PIk(n+1024)/1024) - x(n)*exp(-j2PIkn/1024)] n=128 1151 = SUM [x(n)*exp(-j2PIkn/1024)] n=128 255 + SUM {[x(n+1024)-x(n)]*exp(-j2PIkn/1024)} n=128 255 = X(k,n1) + SUM {[x(n+1024)-x(n)]*exp(-j2PIkn/1024)} n=128 So the DFT at of the data set starting at sample 256 is equal to the DFT of the data set starting at sample 128 plus something resembling a 128-point DFT of the difference between the "newest" 128 samples and the "oldest" 128 samples. (Note the "1024" in the denominator, however.) Since you are building the DFT gradually, you can define the "oldest" 128 samples to be zeroes: 255 X(k,n2) = X(k,n1) + SUM {x(n+1024)*exp(-j2PIkn/1024)} n=128 And so on.
Intuitive <andrea3891@n_o_s_p_a_m.gmail.com> wrote:
(snip)
 
> At this moment I'm following the advice of Vladimir. I tested > the second approach in matlab and It works. However I don't > understand the his first proposed apporach.
> What is the meaning of: > "The best and the simplest is get the data back by inverse FFTs, > then compite the large FFT." > Can you give me a more simple explanation? (Possible step by step).
That sounds obvious enough for me, but you don't want to do it. It comes from a theorist who doesn't care how long it takes. Given eight FFTs of 128 samples each, do an inverse FFT on each, concatenate the result, and do a 1024 point FFT. Slow, but it does work.
> I don't understand the approach proposed by Greg too. > Can you give a better explanation? (Possible step by step).
-- glen
On Jun 16, 7:44 am, "Intuitive" <andrea3891@n_o_s_p_a_m.gmail.com>
wrote:
> Hello everybody. > > I have got a DSP question for all of you. > I'd like to know if it is possible to obtain a large FFT (for example 1024 > samples) combining small SEQUENTIAL FFTs(128 samples). > > It is really important for me that the input samples are consecutive (so it > is not possible to apply the usual trick of even and odd samples). > I'm working on a real time system with a RT constraint of 128 samples. > Actually I'm accumulating samples for 8 iterations (128*8=1024) until to > obtain the 1024 samples to perform the FFT of the input signal. However I'd > like to decompose the FFT in small sequential FFTs in order to distribute > the compute for each frame. > > If you have any question don't hesitate to contact me. > > I'm looking forward to hearing from you, Andrea.
It is readily accomplished by accounting for time delays via a linear phase shift and zeropadding transforms to length N=8*N0 for N0 = 128.. The following example has N0 = 4 and N = 4*N0 = 16. close all, clear all, clc N0 = 4 N = 4*N0 Fs = 16 dt = 1/Fs t = dt*(0:N-1); t0 = N0*dt ind1 = (1:N0); ind2 = (N0+1:2*N0); ind3 = (2*N0+1:3*N0); ind4 = (3*N0+1:N); rand('state',0); f0 = Fs/8 x = cos(2*pi*(f0*t+rand)); df = Fs/N f = df*(0:N-1); X = fft(x); X1 = fft(x(ind1),N); X2 = exp(-i*2*pi*t0*f).*fft([x(ind2)],N); X3 = exp(-i*4*pi*t0*f).*fft([x(ind3)],N); X4 = exp(-i*6*pi*t0*f).*fft([x(ind4)],N); Xc = X1 + X2 + X3 + X4; check = maxabs(X-Xc) % = 1.0782e-014 Hope this helps. Greg
On Jun 16, 4:44&#2013266080;am, "Intuitive" <andrea3891@n_o_s_p_a_m.gmail.com>
wrote:
> Hello everybody. > > I have got a DSP question for all of you. > I'd like to know if it is possible to obtain a large FFT (for example 1024 > samples) combining small SEQUENTIAL FFTs(128 samples). > > It is really important for me that the input samples are consecutive (so it > is not possible to apply the usual trick of even and odd samples). > I'm working on a real time system with a RT constraint of 128 samples. > Actually I'm accumulating samples for 8 iterations (128*8=1024) until to > obtain the 1024 samples to perform the FFT of the input signal. However I'd > like to decompose the FFT in small sequential FFTs in order to distribute > the compute for each frame.
If you are willing to give up latency and memory usage for a more even distribution of FFT computation over time, you could try a pipelined approach. After collecting 8 groups of 128 samples in 8 time slots, you could start doing a 1024 pt fft in 8 pipe stages, 1 pipe stage per time slot, each pipe stage consisting of 1.25 layers of 1024 FFT butterflys (e.g. 1280 butterflys per time slot). After 8 more pipe stages you will end up with all 10 butterfly layers of your 1st 1024 pt FFT, just as you've finished collecting the data to start your next 1024 point FFT. The latency from 1st data group to FFT result will be 16 very computationally balanced pipe stages, instead of 8 pipe stages with a huge summation lump the 8th stage. The memory usage will also be approximately double since you will be gathering data for one FFT while still computing the previous one. This is a typical technique in hardware design (ASIC or FPGA). IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M http://www.nicholson.com/rhn/dsp.html
On Wed, 16 Jun 2010 06:44:18 -0500, "Intuitive"
<andrea3891@n_o_s_p_a_m.gmail.com> wrote:

>Hello everybody. > >I have got a DSP question for all of you. >I'd like to know if it is possible to obtain a large FFT (for example 1024 >samples) combining small SEQUENTIAL FFTs(128 samples). > >It is really important for me that the input samples are consecutive (so it >is not possible to apply the usual trick of even and odd samples). >I'm working on a real time system with a RT constraint of 128 samples. >Actually I'm accumulating samples for 8 iterations (128*8=1024) until to >obtain the 1024 samples to perform the FFT of the input signal. However I'd >like to decompose the FFT in small sequential FFTs in order to distribute >the compute for each frame. > >If you have any question don't hesitate to contact me. > >I'm looking forward to hearing from you, Andrea. >
Hi, You might want to have a look at: http://www.dsprelated.com/showarticle/63.php Good Luck, [-Rick-]
For some reason, no matter how often or where the question is asked, no one
ever seems to post the full solution to this problem, so here it is:

/////////////////////////////////////////////////////////////////////////////////////////////////////
function Y = fftx(X)
% Y = fftx(X)
%Extended fft by combining smaller fft's recursively.
%URL:
http://www.mathworks.com/company/newsletters/articles/faster-finite-fourier-transforms-matlab.html
%Modified to end recursion at max fft size by Don Gateley

%Maximum allowed fft size (use whatever is appropriate for you.)
	max = 512;

	X = X(:);
	n = length(X);

	if(n <= max)
	%Down to chunks and end of recursion. Use native fft() for chunks or
less.
		Y = fft(X);
	else
	% Recursive divide and conquer.
		w = exp(2 * 1i * pi * (0: n / 2 - 1) / n)';
		u = fftx(X(1: 2: n-1));
		v = fftx(X(2: 2: n)) .* w;
		Y = [u + v; u - v];
	end
end
/////////////////////////////////////////////////////////////////////////////////////////////////////

And the inverse:

/////////////////////////////////////////////////////////////////////////////////////////////////////
function T = ifftx(X)
%ifftx() using fftx()
		n = size(X, 1);
		T = fftx(imag(X) + 1i * real(X));
		T= imag(T) + 1i * real(T);
		T = T / n;
end
/////////////////////////////////////////////////////////////////////////////////////////////////////

	 

_____________________________		
Posted through www.DSPRelated.com
For some reason, no matter how often or where the question is asked, no one
ever seems to post the full solution to this problem, so here it is:

/////////////////////////////////////////////////////////////////////////////////////////////////////
function Y = fftx(X)
% Y = fftx(X)
%Extended fft by combining smaller fft's recursively.
%URL:
http://www.mathworks.com/company/newsletters/articles/faster-finite-fourier-transforms-matlab.html
%Modified to end recursion at max fft size by Don Gateley

%Maximum allowed fft size (use whatever is appropriate for you.)
	max = 512;

	X = X(:);
	n = length(X);

	if(n <= max)
	%Down to chunks and end of recursion. Use native fft() for chunks or
less.
		Y = fft(X);
	else
	% Recursive divide and conquer.
		w = exp(2 * 1i * pi * (0: n / 2 - 1) / n)';
		u = fftx(X(1: 2: n-1));
		v = fftx(X(2: 2: n)) .* w;
		Y = [u + v; u - v];
	end
end
/////////////////////////////////////////////////////////////////////////////////////////////////////

And the inverse:

/////////////////////////////////////////////////////////////////////////////////////////////////////
function T = ifftx(X)
%ifftx() using fftx()
		n = size(X, 1);
		T = fftx(imag(X) + 1i * real(X));
		T= imag(T) + 1i * real(T);
		T = T / n;
end
/////////////////////////////////////////////////////////////////////////////////////////////////////	


_____________________________		
Posted through www.DSPRelated.com
On Mon, 10 Mar 2014 21:20:48 -0500, "DonGateley" <99427@dsprelated>
wrote:

> >For some reason, no matter how often or where the question is asked, no one >ever seems to post the full solution to this problem, so here it is:
Well, you've posted the same thing four times now with three different subject lines, so I'm guessing it's covered.
>///////////////////////////////////////////////////////////////////////////////////////////////////// >function Y = fftx(X) >% Y = fftx(X) >%Extended fft by combining smaller fft's recursively. >%URL: >http://www.mathworks.com/company/newsletters/articles/faster-finite-fourier-transforms-matlab.html >%Modified to end recursion at max fft size by Don Gateley > >%Maximum allowed fft size (use whatever is appropriate for you.) > max = 512; > > X = X(:); > n = length(X); > > if(n <= max) > %Down to chunks and end of recursion. Use native fft() for chunks or >less. > Y = fft(X); > else > % Recursive divide and conquer. > w = exp(2 * 1i * pi * (0: n / 2 - 1) / n)'; > u = fftx(X(1: 2: n-1)); > v = fftx(X(2: 2: n)) .* w; > Y = [u + v; u - v]; > end >end >///////////////////////////////////////////////////////////////////////////////////////////////////// > >And the inverse: > >///////////////////////////////////////////////////////////////////////////////////////////////////// >function T = ifftx(X) >%ifftx() using fftx() > n = size(X, 1); > T = fftx(imag(X) + 1i * real(X)); > T= imag(T) + 1i * real(T); > T = T / n; >end >///////////////////////////////////////////////////////////////////////////////////////////////////// > > >_____________________________ >Posted through www.DSPRelated.com
Eric Jacobsen Anchor Hill Communications http://www.anchorhill.com
On Mon, 10 Mar 2014 21:20:48 -0500, "DonGateley" <99427@dsprelated> 
wrote: 

>For some reason, no matter how often or where the question is asked, no one >ever seems to post the full solution to this problem, so here it is:
.... u = fftx(X(1: 2: n-1)); v = fftx(X(2: 2: n)) .* w; Y = [u + v; u - v]; ... I don't use MATLAB, so I'm a little confused by the above three statements. Now I may be reading it wrong, but don't these statements perform FFTs on non-consecutive elements of the vector X and then combine them into one FFT called Y? And, as the OP stated:
>It is really important for me that the input samples are consecutive (so it >is not possible to apply the usual trick of even and odd samples).
Kevin McGee
On Mon, 10 Mar 2014 12:53:20 -0500, "DonGateley" <99427@dsprelated>
wrote:

>For some reason, no matter how often or where the question is asked, no one >ever seems to post the full solution to this problem, so here it is: > >///////////////////////////////////////////////////////////////////////////////////////////////////// >function Y = fftx(X) >% Y = fftx(X) >%Extended fft by combining smaller fft's recursively.
<snip> Where do these posts come from? Is this forum trolling on dsprelated? Mark