# Calculating an N point FFT by calculating 4X N/4 point FFTs in parallel

I have an FPGA based application where I need to perform 4096 point FFTs in real time on a 1GS/s data stream. Data comes to the FFT from an A/D converter as 4 samples in parallel at 250Mhz. My data consists entirely of real values. I would like the FFT to process 4 real samples per clock. Rather than starting from scratch, I would like to use four 1024 point FFT cores in parallel, and then write some VHDL code to combine the results from the four FFTs into a single 4096 point FFT.

I found this post which has an excellent example: Perform non-power-of-two FFT using ARM CMSIS library

I was able to easily modify that example code to work with 4X 1024 point FFTs rather than 5X 256 point FFTs. At a high level, I understand how this works.


Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = n;             % Length of signal
t = (0:L-1)*T;        % Time vector

x = 0.7*sin(2*pi*50*t);

figure;
plot(x);

fx = fft(x);

% Break down into five signals of 1024 points each, interleaved
p = x(1:4:end);
q = x(2:4:end);
r = x(3:4:end);
s = x(4:4:end);

% FFT each of those. This is a 1024 point power-of-two standard FFT
fp = fft(p);
fq = fft(q);
fr = fft(r);
fs = fft(s);

fp4 = [fp fp fp fp];
fq4 = [fq fq fq fq];
fr4 = [fr fr fr fr];
fs4 = [fs fs fs fs];

fp4 = reshape(fp4,n,1);
fq4 = reshape(fq4,n,1);
fr4 = reshape(fr4,n,1);
fs4 = reshape(fs4,n,1);

% calculate the 4096 twiddle factors
k4 = (0:n-1)';
W4 = exp(-1i*2*pi*k4/n);

% assemble the result

fy4 = fp4 + W4.*fq4 + W4.^2.*fr4 + W4.^3.*fs4;

figure;
plot(abs(fx(1:n/2)));

figure;
plot(abs(fy4(1:n/2)));



I am having trouble understanding, and coming up with a hardware efficient implementation for the complex arithmetic step "fy4 = fp4 + W4.*fq4 + W4.^2.*fr4 + W4.^3.*fs4;" from that example.

This statement does not translate directly to hardware very easily, and I suspect that there are some optimizations that could be done to reduce the computational complexity. I would greatly appreciate it if someone could help me understand how to re-write that step of the algorithm into a form that would translate more easily into hardware.

I am looking for an explanation that is similar to how the radix-2 butterfly is described below, but for the butterfly that I need to implement to combine four N/4 point FFTs into a single N point FFT.

Thank you!

Update:

Below is a version of Hilmar's code that generates two samples per loop. I also separated out the real and imaginary components since the hardware implementation can only handle real arithmetic.

I plan to calculate power spectra from the FFT results, so I only need to keep the first N/2 points from the FFT. Therefore I only need to calculate two output points for every four input points.

This works, and it is in a state where I can translate it to VHDL. It uses 10 lookup tables (5 sine, 5 cosine), and 24 multiplies per loop. Because the lookup tables will be implemented in FPGA block RAM, I cannot really take advantage of the circular addressing trick. I need all of the twiddle factors to be available on every clock cycle.

I still have a suspicion that there is a more efficient way to do this. Are there simplifications that would reduce the number of operations, and reduce the number of twiddle factor lookup tables that I need?

I would also like to understand if this operations is the same as a Radix-4 butterfly. The references that I have seen on the radix-4 butterfly indicate that it uses fewer lookup tables and fewer multiplications than this solution, but I do not understand how to get from one to the other.

n = 4096;

Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = n;                % Length of signal
t = (0:L-1)*T;        % Time vector

x = 0.7*sin(2*pi*50*t)*(2^16);

figure;
plot(x);

% calculate FFT using MATLAB native fft() function.
% We'll use this as a reference to prove it works
fx = fft(x);

% Break down into four signal of 1024 points each, interleaved
p = x(1:4:end);
q = x(2:4:end);
r = x(3:4:end);
s = x(4:4:end);

% FFT each of those. This is a 1024 power-of-two standard FFT
fp = fft(p);
fq = fft(q);
fr = fft(r);
fs = fft(s);

fp4 = [fp fp fp fp];
fq4 = [fq fq fq fq];
fr4 = [fr fr fr fr];
fs4 = [fs fs fs fs];

fp4 = reshape(fp4,n,1);
fq4 = reshape(fq4,n,1);
fr4 = reshape(fr4,n,1);
fs4 = reshape(fs4,n,1);

% calculate the 4096 twiddle factors
k4 = (0:n-1)';
W4 = exp(-1i*2*pi*k4/n);

% assemble the result

fy4 = fp4 + W4.*fq4 + W4.^2.*fr4 + W4.^3.*fs4;

figure;
plot(abs(fy4(1:n/2)));

%use sines and cosine instead of exp
C  = cos(2*pi*k4/n);
C2 = cos(2*pi*k4*2/n);
C3 = cos(2*pi*k4*3/n);

S  = -sin(2*pi*k4/n);
S2 = -sin(2*pi*k4*2/n);
S3 = -sin(2*pi*k4*3/n);

fy4a = 0*fy4;
fy4b = 0*fy4;

s = 2^20;   %Scaling factor for integer lookup tables

for i = 1:n/4

fy4a(i) = fp4(i) + W4(i)*fq4(i) + W4(i)^2*fr4(i) + W4(i)^3*fs4(i);

xa = real(fp4(i));
ya = imag(fp4(i));
xb = real(fq4(i));
yb = imag(fq4(i));
xc = real(fr4(i));
yc = imag(fr4(i));
xd = real(fs4(i));
yd = imag(fs4(i));

War = round(C(i)*s);
Wai = round(S(i)*s);
Wbr = round(C2(i)*s);
Wbi = round(S2(i)*s);

Wcr = round(C3(i)*s);
Wci = round(S3(i)*s);

War2 = round(C(i+n/4)*s);
Wai2 = round(S(i+n/4)*s);

%Can resuse the C2 value from the first calculation
%Saves two lookup tables.
%Wbr2 = round(C2(i+n/4)*s);
Wbr2 = round(-C2(i)*s);
%Wbi2 = round( S2(i+n/4)*s);
Wbi2 = round( -S2(i)*s);

Wcr2 = round(C3(i+n/4)*s);
Wci2 = round(S3(i+n/4)*s);

%Calculate Intermediate terms.  This will be pipe stage 1 in the VHDL
%divide by scaling factor and round to simulate fixed point math
Waixb = round((Wai*xb)/s);
Waiyb = round((Wai*yb)/s);
Warxb = round((War*xb)/s);
Waryb = round((War*yb)/s);
Wbixc = round((Wbi*xc)/s);
Wbiyc = round((Wbi*yc)/s);
Wbrxc = round((Wbr*xc)/s);
Wbryc = round((Wbr*yc)/s);
Wcixd = round((Wci*xd)/s);
Wciyd = round((Wci*yd)/s);
Wcrxd = round((Wcr*xd)/s);
Wcryd = round((Wcr*yd)/s);

Wai2xb = round((Wai2*xb)/s);
Wai2yb = round((Wai2*yb)/s);
War2xb = round((War2*xb)/s);
War2yb = round((War2*yb)/s);
Wbi2xc = round((Wbi2*xc)/s);
Wbi2yc = round((Wbi2*yc)/s);
Wbr2xc = round((Wbr2*xc)/s);
Wbr2yc = round((Wbr2*yc)/s);
Wci2xd = round((Wci2*xd)/s);
Wci2yd = round((Wci2*yd)/s);
Wcr2xd = round((Wcr2*xd)/s);
Wcr2yd = round((Wcr2*yd)/s);

Xr = xa + (Warxb - Waiyb) + (Wbrxc - Wbiyc) + (Wcrxd - Wciyd);
%Xi = ya + ((War+Wai)*(xb+yb) - Warxb - Waiyb) +                        ((Wbr+Wbi)*(xc+yc) - Wbrxc - Wbiyc)                   + ((Wcr+Wci)*(xd+yd) - Wcrxd - Wci*yd);

%Xi = ya + (     (War*xb + Wai*xb + War*yb + Wai*yb)  - Warxb - Waiyb) + ((Wbr*xc + Wbi*xc + Wbr*yc + Wbi*yc)  - Wbrxc - Wbiyc) + ((Wcr*xd + Wcr*yd + Wci*xd + Wci*yd ) - Wcrxd - Wci*yd);

%Xi = ya + ( Warxb + Waixb + Waryb + Waiyb  - Warxb - Waiyb + Wbrxc + Wbixc + Wbryc + Wbiyc  - Wbrxc - Wbiyc + Wcrxd + Wcryd + Wcixd + Wciyd - Wcrxd - Wciyd);
Xi = ya + (  Waixb + Waryb + Wbixc + Wbryc  + Wcryd + Wcixd);

%Yr = xa + (War2*xb - Wai2*yb) + (Wbr2*xc - Wbi2*yc) + (Wcr2*xd - Wci2*yd);
Yr = xa + (War2xb - Wai2yb) + (Wbr2xc - Wbi2yc) + (Wcr2xd - Wci2yd);

%Yi = ya + ((War2+Wai2)*(xb+yb) - War2*xb - Wai2*yb) + ((Wbr2+Wbi2)*(xc+yc) - Wbr2*xc - Wbi2*yc) + ((Wcr2+Wci2)*(xd+yd) - Wcr2*xd - Wci2*yd);
%Yi = ya + ( (War2xb + Wai2xb + War2yb + Wai2yb)  - War2xb - Wai2yb) + ((Wbr2xc + Wbi2xc + Wbr2yc + Wbi2yc)  - Wbr2xc - Wbi2yc) + ((Wcr2xd + Wcr2yd + Wci2xd + Wci2yd ) - Wcr2xd - Wci2yd);
Yi = ya + (  Wai2xb + War2yb + Wbi2xc + Wbr2yc  + Wcr2yd + Wci2xd);

fy4b(i) = complex(Xr,Xi);
fy4b(i+n/4) = complex(Yr,Yi);

end

figure;
plot(abs(fy4a(1:n/2)));

figure;
plot(abs(fy4b(1:n/2)));  


• What finds FPGA is this? Considering the data rates, I'd be surprised if no 4096-FFT IP core existed for it.
– mmmm
Aug 16, 2021 at 23:55
• Unfortunately I think you will have to workout the equations. If you receive 4 samples per clock the first question I would ask is if you really need the 4096 point FFT.
– Bob
Aug 17, 2021 at 15:49
• The FPGA is a Xilinx Ultrascale+. Xilinx has a free FFT core that can meet timing at 500Mhz, but it is limited to a single complex input sample per clock cycle. I could keep up with the data rate by time slicing two 4096 pt FFT cores in parallel, but latency is an issue for my application, therefore I need a faster FFT. Aug 17, 2021 at 16:01

Let me try to describe in words how this works:

1. Calculate 4 individual 1k FFTs
2. Repeat each 4 times to make it a 4k length
3. Multiply each result with a vector of the twiddle factors raised to a power: 0 for the first , 1 for the second, etc.
4. Sum them up

I primarily did this way since I was lazy. We can make this a lot more hardware friendly by using circular addressing, e.g. for base 1024, we would count 1021, 1022, 1023, 0, 1 ,2 ... Since all the bases are a power of two, that can be implement with a simple bit-wise and with n-1.

We can also make use of the fact that taking the power of a twiddle factor is the same as multiplying the index, i.e.

$$W(k)^m = W(k\cdot m)$$

So instead of taking the power, we can just use a different step size, provided we use it with circular addressing. Here is a snippet of Matlab code that demonstrates this by unrolling the last summation step

%% unroll the final loop and only use 1k length input vectors
fy4a = 0*fy4;
M = 1; % Matlab array indexing offset
for i = 1:n
realIndex = i-1; % remove matlab indexing offset
k = M + bitand(realIndex,moduloMaskShort); % mod 1024
% add FFTs 0 and 1
fy4a(i) = fp4(k)+W4(i)*fq4(k);
% FFt 2 with a step size of 2
% FFT 3 with a step side of 3
end


This is ugly Matlab code primarily since Matlab starts counting at 1 which makes the whole circular addressing awkward, but should work fine on hardware or C. Please note that this replaces steps 2, 3 & 4 above, so there is no need to replicate the FFT results tp full 4k length. It works with the 1k results "as is".

This can be done as a radix 4 operation. Below is the Matlab. Note that you only need 3 twiddle factors per butterfly.

I implemented the table lookup by have three different pointers with different step sizes (1,2 & 3). You only need to table up 3069 twiddle factors and not the whole 4095.

The three twiddle factors are related as $$W_2 = W_1 \cdot W_1$$ and $$W_3 = W_1 \cdot W_2$$, so if complex multiplication is faster than table lookup, you can do that. If you do, you only need a twiddle factor table up to 1023.

Multiplication with $$j$$ or $$-j$$ doesn't require any actual multiplications, just swapping real and imaginary parts and flipping the proper signs. %% do it as a 4in 4out operation fy4b = 0*fy4; % initialze n4 = n/4;

M = 1; % Matlab array offset
i2 = 0; % index for W^2
i3 = 0; % index for W^3
j = 1i; % imaginary unit

for i1 = 0:n4-1
im =  i1 + M; % index into Matlab arrays starting at 1

% get the tiwddle factors and multipy with inputs
a0 = fp4(im);
a1 = fq4(im)*W4(M+i1);
a2 = fr4(im)*W4(M+i2);
a3 = fs4(im)*W4(M+i3);

% perform the radix as 4 indivdiual operations
fy4b(im)      = a0 +   a1 + a2 +   a3;
fy4b(im+n4)   = a0 - j*a1 - a2 + j*a3;
fy4b(im+2*n4) = a0 -   a1 + a2 -   a3;
fy4b(im+3*n4) = a0 + j*a1 - a2 - j*a3;

% update the twiddle factor indices counters
i2 = i2 + 2;
i3 = i3 + 3;

end

d = (fy4b-fy4);

fprintf('Relative Error = %6.2fdB \n',20*log10(sum(abs(d))./sum(abs(fy4))));
`
• I was not clear in my question, but what I want to do is calculate 4 output points in parallel in an efficient manner so that the FPGA logic that is combining the four 1024 point FFT's can keep up with the data streaming out of the four 1024 point FFT cores. I am using the FFT results to produce power spectra, so I only care about the first 2048 points of the FFT output( if that helps with the problem) Is this calculation fundamentally the same as a radix-4 butterfly? If so, why does the data not need to be re-ordered? Thank you for your help. Aug 17, 2021 at 15:56
• Sorry, I don't understand your question. Nothing needs to re-ordered, the last loop is completely sequential and you can certainly chop it up in chunks of 4 if you like. They way it's currently written, it's a 4-in 1-out operation that's looped 4096 times. You can potentially rewrite as a 4-in 4-out operation that's looped 1024 times, but then your output is not sequential. You would get it as 0,1024,2048,3072,1,1025,2049, ... If that;s what you want I can write this up. Aug 17, 2021 at 16:09
• Yes, I am looking for a 4 in 4 out operation, and also wondering if there are any optimizations that can be made, for instance will I need 12 different W terms (which means 12 lookup tables) in the 4-in 4-out case, or can W terms be shared? I think that this is basically the same as a radix-4 butterfly and there are radix-4 butterfly topologies that only need four W values. Thanks again for your help. My engineer brain is pretty good with implementation details, but terrible on abstract ideas and theory. Aug 17, 2021 at 16:48
• Can the equations derived for the radix-4 butterfly on page 12 of the following TI appnote be used for the "fy4 = fp4 + W4.*fq4 + W4.^2.*fr4 + W4.^3.*fs4;" calculation? ti.com/lit/an/spra152/spra152.pdf Aug 17, 2021 at 18:46
• @robertbristow-johnson: correct. But I my initial recombination algorithm wasn't using a standard butterfly structure (yet), so it actually did wrap around. I fully agree that using a radix 4 butterfly (or whatever insect has four wings on each side) is the better choice and eliminates the need for circular addressing. Aug 19, 2021 at 20:33