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
Reply by kevinjmcee●March 11, 20142014-03-11
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
Reply by Eric Jacobsen●March 11, 20142014-03-11
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
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
Reply by DonGateley●March 10, 20142014-03-10
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
Reply by Rick Lyons●June 20, 20102010-06-20
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.
>
On Jun 16, 4: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 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
Reply by Greg Heath●June 18, 20102010-06-18
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
Reply by glen herrmannsfeldt●June 17, 20102010-06-17
> 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
Reply by Greg Berchin●June 17, 20102010-06-17
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.