Update FFT code to include full testbench

LadyPumie · Jul 30, 2018 · a9bf74f · a9bf74f
1 parent 637c025
commit a9bf74f
Show file tree

Hide file tree

Showing 13 changed files with 620 additions and 56 deletions.
diff --git a/examples/fft.h b/examples/fft.h
@@ -0,0 +1,19 @@
+#ifndef FFT_H_
+#define FFT_H_
+
+typedef float DTYPE;
+typedef int INTTYPE;
+#define M 10 			/* Number of Stages = Log2N */
+#define SIZE 1024 		/* SIZE OF FFT */
+#define SIZE2 SIZE>>1	/* SIZE/2 */
+
+void fft(DTYPE XX_R[SIZE], DTYPE XX_I[SIZE]);
+
+//W_real and W_image are twiddle factors for 1024 size FFT.
+//WW_R[i]=cos(e*i/SIZE);
+//WW_I[i]=sin(e*i/SIZE);
+//where i=[0,512) and DTYPE	e = -6.283185307178;
+#include "tw_r.h"
+#include "tw_i.h"
+
+#endif
diff --git a/examples/fft_bit_reverse.c b/examples/fft_bit_reverse.c
diff --git a/examples/fft_stages-top.cpp b/examples/fft_stages-top.cpp
@@ -0,0 +1,103 @@
+/*
+This is traditional 2-radix DIT FFT algorithm implementation.
+It is based on conventional 3-loop structure. 
+INPUT:
+	In_R, In_I[]: Real and Imag parts of Complex signal
+
+OUTPUT:
+	In_R, In_I[]: Real and Imag parts of Complex signal
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <math.h>
+#include "fft_streaming.h"
+
+DTYPE In_R[SIZE], In_I[SIZE];
+DTYPE Out_R[SIZE], Out_I[SIZE];
+DTYPE WW_R[SIZE], WW_I[SIZE];
+
+int main()
+{
+	FILE *fp;
+	FILE *fp_r, *fp_i;
+	printf("GENERATING INPUTS\n");
+	for(int i=0; i<SIZE; i++){
+		In_R[i] = i;
+		In_I[i] = 0.0;
+	}
+
+	//Twiddle factor is calculated here and saved in fft.h to be used in offline.
+		double	e = -6.2831853071795864769;
+		printf("GENERATING %d TWIDDLE FACTORS\n", SIZE);
+		fp_r=fopen("tw_r.h", "w");
+		fp_i=fopen("tw_i.h", "w");
+		fprintf(fp_r, "const DTYPE W_real[]={");
+		fprintf(fp_i, "const DTYPE W_imag[]={");
+		for(int i=0; i<SIZE2; i++)
+		{
+			//COMPLEX W;	// e^(-j 2 pi/ N)
+		  double w = e*double(i)/double(SIZE);
+		  WW_R[i]=cos(w);
+		  WW_I[i]=sin(w);
+		  //printf("%4d\t%f\t%f\n",i,WW_R[i],WW_I[i]);
+			fprintf(fp_r, "%.20f,",WW_R[i]);
+			fprintf(fp_i, "%.20f,",WW_I[i]);
+			if(i%16==0)
+				{
+					fprintf(fp_r, "\n");
+					fprintf(fp_i, "\n");
+				}
+		}
+		fprintf(fp_r, "};\n");	
+		fprintf(fp_i, "};\n");
+		fclose(fp_r);
+		fclose(fp_i);
+
+
+	//Perform FFT
+		fft_streaming(In_R, In_I, Out_R, Out_I);
+	//Print output
+	fp=fopen("out.fft.dat", "w");
+	printf("Printing FFT Output\n");
+	for(int i=0; i<SIZE; i++){
+	  //printf("%4d\t%f\t%f\n",i,Out_R[i],Out_I[i]);
+		fprintf(fp, "%4d\t%f\t%f\n",i,Out_R[i],Out_I[i]);
+	}
+	fclose(fp);
+
+
+	printf ("Comparing against output data \n");
+	std::ifstream golden("out.fft.gold.dat");
+
+	DTYPE error = 0.0;
+	DTYPE maxerror = 0.0;
+	for(int i=0; i<SIZE; i++) {
+	  DTYPE rx, ix;
+	  int j;
+	  golden >> j >> rx >> ix;
+	  DTYPE newerror = fabs(rx-Out_R[i]) + fabs(ix-Out_I[i]);
+	  error += newerror;
+	  if(newerror > maxerror) {
+	    maxerror = newerror; 
+	    fprintf(stdout, "Max Error@%d: %f\n", i, maxerror);
+	  }
+	}
+
+	fprintf(stdout, "Average Error: %f\n",error/SIZE);
+
+	if ((error/SIZE) > .05 || maxerror > 2) { // This is somewhat arbitrary.  Should do proper error analysis.
+	  fprintf(stdout, "*******************************************\n");
+	  fprintf(stdout, "FAIL: Output DOES NOT match the golden output\n");
+	  fprintf(stdout, "*******************************************\n");
+	  return 1;
+	} else {
+	  fprintf(stdout, "*******************************************\n");
+	  fprintf(stdout, "PASS: The output matches the golden output!\n");
+	  fprintf(stdout, "*******************************************\n");
+	  return 0;
+	}
+
+}
diff --git a/examples/fft_stages.cpp b/examples/fft_stages.cpp
@@ -1,21 +1,75 @@
+#include "fft_streaming.h"
+#include "math.h"
+
+unsigned int reverse_bits(unsigned int input) {
+	int i, rev = 0;
+	for (i = 0; i < M; i++) {
+		rev = (rev << 1) | (input & 1);
+		input = input >> 1;
+	}
+	return rev;
+}
+
 void bit_reverse(DTYPE X_R[SIZE], DTYPE X_I[SIZE],
-        DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE]);
-void fft_stage_one(DTYPE X_R[SIZE], DTYPE X_I[SIZE],
-        DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE]);
-void fft_stages_two(DTYPE X_R[SIZE], DTYPE X_I[SIZE],
-        DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE]);
-void fft_stage_three(DTYPE X_R[SIZE], DTYPE X_I[SIZE],
-        DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE]);
+		 DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE]) {
+  unsigned int reversed;
+  unsigned int i;
+  DTYPE temp;
+
+  for (int i = 0; i < SIZE; i++) {
+	  reversed = reverse_bits(i); // Find the bit reversed index
+		if (i <= reversed) {
+			// Swap the real values
+			temp = X_R[i];
+			OUT_R[i] = X_R[reversed];
+			OUT_R[reversed] = temp;
 
-void fft(DTYPE X_R[SIZE], DTYPE X_I[SIZE], DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE])
+			// Swap the imaginary values
+			temp = X_I[i];
+			OUT_I[i] = X_I[reversed];
+			OUT_I[reversed] = temp;
+		}
+	}
+}
+
+void fft_stage(int stage, DTYPE X_R[SIZE], DTYPE X_I[SIZE],
+		     DTYPE Out_R[SIZE], DTYPE Out_I[SIZE]) {
+  int DFTpts = 1 << stage;    // DFT = 2^stage = points in sub DFT
+  int numBF = DFTpts / 2;     // Butterfly WIDTHS in sub-DFT
+  int step = SIZE >> stage;
+  DTYPE k = 0;
+  DTYPE e = -6.283185307178 / DFTpts;
+  DTYPE a = 0.0;
+  // Perform butterflies for j-th stage
+ butterfly_loop:
+  for (int j = 0; j < numBF; j++) {
+    DTYPE c = cos(a);
+    DTYPE s = sin(a);
+    a = a + e;
+    // Compute butterflies that use same W**k
+  dft_loop:
+    for (int i = j; i < SIZE; i += DFTpts) {
+      int i_lower = i + numBF; // index of lower point in butterfly
+      DTYPE temp_R = X_R[i_lower] * c - X_I[i_lower] * s;
+      DTYPE temp_I = X_I[i_lower] * c + X_R[i_lower] * s;
+      Out_R[i_lower] = X_R[i] - temp_R;
+      Out_I[i_lower] = X_I[i] - temp_I;
+      Out_R[i] = X_R[i] + temp_R;
+      Out_I[i] = X_I[i] + temp_I;
+    }
+    k += step;
+  }
+}
+
+void fft_streaming(DTYPE X_R[SIZE], DTYPE X_I[SIZE], DTYPE OUT_R[SIZE], DTYPE OUT_I[SIZE])
 {
   #pragma HLS dataflow
-  DTYPE Stage1_R[SIZE], Stage1_I[SIZE];
-  DTYPE Stage2_R[SIZE], Stage2_I[SIZE];
-  DTYPE Stage3_R[SIZE], Stage3_I[SIZE];
-
+  DTYPE Stage1_R[SIZE], Stage1_I[SIZE],
+    Stage2_R[SIZE], Stage2_I[SIZE],
+    Stage3_R[SIZE], Stage3_I[SIZE];
+  
   bit_reverse(X_R, X_I, Stage1_R, Stage1_I);
-  fft_stage_one(Stage1_R, Stage1_I, Stage2_R, Stage2_I);
-  fft_stages_two(Stage2_R, Stage2_I, Stage3_R, Stage3_R);
-  fft_stage_three(Stage3_R, Stage3_I, OUT_R, OUT_I);
+  fft_stage(1, Stage1_R, Stage1_I, Stage2_R, Stage2_I);
+  fft_stage(2, Stage2_R, Stage2_I, Stage3_R, Stage3_I);
+  fft_stage(3, Stage3_R, Stage3_I, OUT_R, OUT_I);
 }
diff --git a/examples/fft_stages_loop-top.cpp b/examples/fft_stages_loop-top.cpp
@@ -0,0 +1,103 @@
+/*
+This is traditional 2-radix DIT FFT algorithm implementation.
+It is based on conventional 3-loop structure. 
+INPUT:
+	In_R, In_I[]: Real and Imag parts of Complex signal
+
+OUTPUT:
+	In_R, In_I[]: Real and Imag parts of Complex signal
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <math.h>
+#include "fft_streaming.h"
+
+DTYPE In_R[SIZE], In_I[SIZE];
+DTYPE Out_R[SIZE], Out_I[SIZE];
+DTYPE WW_R[SIZE], WW_I[SIZE];
+
+int main()
+{
+	FILE *fp;
+	FILE *fp_r, *fp_i;
+	printf("GENERATING INPUTS\n");
+	for(int i=0; i<SIZE; i++){
+		In_R[i] = i;
+		In_I[i] = 0.0;
+	}
+
+	//Twiddle factor is calculated here and saved in fft.h to be used in offline.
+		double	e = -6.2831853071795864769;
+		printf("GENERATING %d TWIDDLE FACTORS\n", SIZE);
+		fp_r=fopen("tw_r.h", "w");
+		fp_i=fopen("tw_i.h", "w");
+		fprintf(fp_r, "const DTYPE W_real[]={");
+		fprintf(fp_i, "const DTYPE W_imag[]={");
+		for(int i=0; i<SIZE2; i++)
+		{
+			//COMPLEX W;	// e^(-j 2 pi/ N)
+		  double w = e*double(i)/double(SIZE);
+		  WW_R[i]=cos(w);
+		  WW_I[i]=sin(w);
+		  //printf("%4d\t%f\t%f\n",i,WW_R[i],WW_I[i]);
+			fprintf(fp_r, "%.20f,",WW_R[i]);
+			fprintf(fp_i, "%.20f,",WW_I[i]);
+			if(i%16==0)
+				{
+					fprintf(fp_r, "\n");
+					fprintf(fp_i, "\n");
+				}
+		}
+		fprintf(fp_r, "};\n");	
+		fprintf(fp_i, "};\n");
+		fclose(fp_r);
+		fclose(fp_i);
+
+
+	//Perform FFT
+		fft_streaming(In_R, In_I, Out_R, Out_I);
+	//Print output
+	fp=fopen("out.fft.dat", "w");
+	printf("Printing FFT Output\n");
+	for(int i=0; i<SIZE; i++){
+	  //printf("%4d\t%f\t%f\n",i,Out_R[i],Out_I[i]);
+		fprintf(fp, "%4d\t%f\t%f\n",i,Out_R[i],Out_I[i]);
+	}
+	fclose(fp);
+
+
+	printf ("Comparing against output data \n");
+	std::ifstream golden("out.fft.gold.dat");
+
+	DTYPE error = 0.0;
+	DTYPE maxerror = 0.0;
+	for(int i=0; i<SIZE; i++) {
+	  DTYPE rx, ix;
+	  int j;
+	  golden >> j >> rx >> ix;
+	  DTYPE newerror = fabs(rx-Out_R[i]) + fabs(ix-Out_I[i]);
+	  error += newerror;
+	  if(newerror > maxerror) {
+	    maxerror = newerror; 
+	    fprintf(stdout, "Max Error@%d: %f\n", i, maxerror);
+	  }
+	}
+
+	fprintf(stdout, "Average Error: %f\n",error/SIZE);
+
+	if ((error/SIZE) > .05 || maxerror > 2) { // This is somewhat arbitrary.  Should do proper error analysis.
+	  fprintf(stdout, "*******************************************\n");
+	  fprintf(stdout, "FAIL: Output DOES NOT match the golden output\n");
+	  fprintf(stdout, "*******************************************\n");
+	  return 1;
+	} else {
+	  fprintf(stdout, "*******************************************\n");
+	  fprintf(stdout, "PASS: The output matches the golden output!\n");
+	  fprintf(stdout, "*******************************************\n");
+	  return 0;
+	}
+
+}