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

Question

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)));  ```
 

Answer 1

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;
% modulo 4096 mask
moduloMask = n-1;
moduloMaskShort = n/4-1; % mod 1024 mask
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
  fy4a(i) = fy4a(i) + W4(M+bitand(2*realIndex,moduloMask))*fr4(k);
  % FFT 3 with a step side of 3
  fy4a(i) = fy4a(i) + W4(M+bitand(3*realIndex,moduloMask))*fs4(k);
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".

EDIT: Radix 4 formulation

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))));