For some reason, no matter how often or where the question is asked, no one
ever posts 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 ●March 10, 20142014-03-10
For some reason, no matter how often or where the question is asked, no one ever posts 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
/////////////////////////////////////////////////////////////////////////////////////////////////////
Reply by Afinko●January 14, 20102010-01-14
>On Jun 1, 8:06 pm, "digintu" <digi...@gmail.com> wrote:
>you're missing a step or two. if you're gonna stuff two halves of
>your length 128 data into two length 64 FFTs and post-process the
>results, you have to put the even-indexed samples in one and the odd-
>index samples in another. do two 64-pt FFTs, and post process the
>results with 64 butterflies that have coefficients of
>
> W(k) = exp(j*2*pi*k/N) for 0 <= k < N/2
>
>where N is 128. this is called the Decimation-In-Time algorithm.
>
>if you want to do this in normal order and toss the "first half" in
>one FFT and the "last half" in the other FFT, that's the Decimation-In-
>Frequency algorithm and instead of post-processing the results, you
>*pre*-process what goes in with similar looking butterflies. either
>the DIT or DIF can be done in-place where the input is in normal or
>bit-reversed order and the output would then naturally end up bit-
>reversed or normal, respectively. in fact, the DIT with normal input,
>bir-reversed output, can be compared directly to the DIF with bit-
>reversed input and normal order output. a butterfly in the final
>stage (or "pass") of the DIT simply has its function inverted with the
>corresponding butterfly in the first pass of the DIF. (assume the
>first DIT alg is a "forward DFT" and the second DIF alg is set for
>"inverse DFT". "j" in the first becomes "-j" in the second.)
>
>r b-j
Hi Robert bristow-johnson,
Can you pleas give me more information about preprocessing by DIT or DIF,
I want to implement it in DSP. The code for MATLAB could be great. I need
to know what are the steps for preprocessing, and then after FFT how to
combine them.
Thanks in advance.
Afi
Reply by Rick Lyons●August 16, 20092009-08-16
On Sun, 01 Jun 2008 19:06:00 -0500, "digintu" <digintu@gmail.com>
wrote:
>Hello,
>
>When combining two smaller FFTs to make a larger one, what are th
>techniques. Example..I want to combine two 64 pt FFTs into a 128 pt FFT.
>have fed the first 64 data samples into the first FFT and the second 6
>into the 2nd FFT. The results are output from the two FFTs into a 12
>array. The problem is that the 128 pt FFT result has the same peaks as th
>64 pt FFT. What is the problem here? How do I correct this in order to ge
>the same results as the standard 128 pt FFT?
>
>Thank you
>DJ
>>Then, X[k] = X_e[k] + exp(-j*2*pi*k/128) * X_o[k], k = 0,...,127,
>where
>>X_e and X_o are the 64-point DFTs of x_e and x_o, respectively. You
>should
>>take the mod(64) values of k on the right side, i.e. rem(k,64) in
Matlab
>>language, since X_e and X_o are periodic with 64.
>
>Just for the heck of it, here is the verifying Matlab code:
>
>>> x=randn(128,1);
>>> xe = x(1:2:end);
>>> xo = x(2:2:end);
>>> X = fft(x);
>>> Xe = fft(xe);
>>> Xo = fft(xo);
>>> Xel = [Xe;Xe] + exp(-j*2*pi*[0:127].'/128) .* [Xo;Xo];
>>> max(abs(Xel-X))
>
>ans =
>
> 8.8818e-15
>
>
Reply by emre●June 17, 20082008-06-17
digintu asks:
I'm trying to implement this in an FPGA so it takes a bit of time. Another
question...If I had a 256 pt FFT, would it be easier to get the results for
a 128 pt result out of it...ie have two 128 pt results.
any ideas?
>I initially had them sequential x[n] n = 1...64 for the first and n =
>65..128 for the second.
>
>I can change them to suit any method.
Bob, if you take the even and odd samples as your 64-point sequences, the
result is elegantly simple.
Say x_e[n] = x[2n], and x_o = x[2n+1], n = 0,...,63.
Then, X[k] = X_e[k] + exp(-j*2*pi*k/128) * X_o[k], k = 0,...,127, where
X_e and X_o are the 64-point DFTs of x_e and x_o, respectively. You should
take the mod(64) values of k on the right side, i.e. rem(k,64) in Matlab
language, since X_e and X_o are periodic with 64.
You can verify this by partitioning the sum for X[k], the 128-point DFT,
into even and odd samples.
Hope this helps.
Emre
Reply by digintu●June 12, 20082008-06-12
>It is not clear from your posting how the two 64 pt sequences are
related
>in time.
>
>Are they sequential, i.e., {x[n]} n=1,...,64 for the first, and
>n=65,...,128 for the second?
>
>Or are they even and odd samples of a 128 point sequence?
>
>emre
>
I initially had them sequential x[n] n = 1...64 for the first and n =
65..128 for the second.
I can change them to suit any method. If you look at my post from earlier
today in response to Robert, I currently have the input indexes bit
reversed and passed through 64 radix2 butterflys. I'm now looking for the
correct way to present that data to the radix4 ffts.
Bob
Reply by emre●June 12, 20082008-06-12
>When combining two smaller FFTs to make a larger one, what are the
>techniques. Example..I want to combine two 64 pt FFTs into a 128 pt FFT.
I
>have fed the first 64 data samples into the first FFT and the second 64
>into the 2nd FFT. The results are output from the two FFTs into a 128
>array. The problem is that the 128 pt FFT result has the same peaks as
the
>64 pt FFT. What is the problem here? How do I correct this in order to
get
>the same results as the standard 128 pt FFT?
It is not clear from your posting how the two 64 pt sequences are related
in time.
Are they sequential, i.e., {x[n]} n=1,...,64 for the first, and
n=65,...,128 for the second?
Or are they even and odd samples of a 128 point sequence?
emre