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