diff options
Diffstat (limited to 'FFT_Test/fft.hcc')
| -rw-r--r-- | FFT_Test/fft.hcc | 409 |
1 files changed, 409 insertions, 0 deletions
diff --git a/FFT_Test/fft.hcc b/FFT_Test/fft.hcc new file mode 100644 index 0000000..82a1a29 --- /dev/null +++ b/FFT_Test/fft.hcc @@ -0,0 +1,409 @@ +#include <pal_master.hch>
+#include <stdlib.hch>
+#include "weights256.hch"
+#include "config.hch"
+#include "debug.hch"
+#include "xilinxmult.hch"
+
+#define PERFORM_FFT_CALCULATION 1
+#define USE_UNSIGNED_AUDIO 0
+#define HARDWARE_MULTIPLY 1
+#define PRINT_DEBUG 0
+/* Define two multi-port RAMs for FFT calculation; one for real and one for imaginary values
+ * Extra block RAM settings are defined to make sure read and write actions can be performed
+ * within one clock-cycle.
+ * Left out extra settings on new board the clock changes TODO !!!!
+ */
+mpram
+{
+ ram signed 24 rwrite[256];
+ rom signed 24 read[256];
+} real with {block = "BlockRAM"/*, westart=2.5, welength=1, rclkpos={1.5}, wclkpos={3}, clkpulselen=0.5*/};
+
+mpram
+{
+ ram signed 24 rwrite[256];
+ rom signed 24 read[256];
+} imaginary with {block = "BlockRAM"/*, westart=2.5, welength=1, rclkpos={1.5}, wclkpos={3}, clkpulselen=0.5*/};
+
+// multiplication factors for equalizer function
+ram signed 7 eq_settings[16] = {0,2,4,7,10,13,16,19,22,26,30,35,41,48,55,63};
+
+
+/****************************************************************
+* Function: multiply *
+* *
+* Arguments *
+* x,y signed variables *
+* *
+* Description *
+* Just a multiplier. But by doing this in a function the *
+* FPGA space needed is reduced. *
+* *
+* Return Values *
+* The result after multiplication *
+****************************************************************/
+macro proc multiply(result, op_a, op_b)
+{
+#if HARDWARE_MULTIPLY
+ xilinxmult(result, adjs(op_a,18), adjs(op_b,18));
+#else
+ macro expr mul_width = (width(op_a) + width(op_b));
+ result = adjs(op_a, mul_width) * adjs(op_b, mul_width);
+#endif
+/* signed result;
+
+ fixed_mul(result,x,y);
+ return result;*/
+}
+/*macro proc multiply(result, x,y)
+{
+#if HARDWARE_MULTIPLY
+ xilinxmult(result, adjs(x,18), adjs(y,18));
+#else
+ result = ((adjs(x,36))*(adjs(y,36)));
+#endif
+
+}
+*/
+
+
+
+/*******************************************************************
+* Function: calculate_fft *
+* *
+* Arguments *
+* select_inverse Boolean that indicates FFT or iFFT calculation *
+* *
+* Description *
+* This routine performs the Fast Fourier Transform for *
+* calculation of the frequency spectrum *
+* *
+*******************************************************************/
+void calculate_fft(unsigned 1 select_inverse)
+{
+ unsigned 4 level;
+ unsigned 8 point1,point2,j,f,k;
+ unsigned 9 e,i;
+ signed 16 weight1,weight2;
+ signed 24 p,q,r,t;
+ signed a,b;
+
+ macro expr rescale (x) = (x[35] @ x[28:14]);
+
+ for(level=1;level<=NUMBER_OF_COLUMNS;level++) // count all the columns
+ {
+ e=1<<(NUMBER_OF_COLUMNS-level+1); // number of points in each block in this column
+ f=(e>>1)<-8; // number of butterflies in each block in this column
+
+ for(j=1;j<=f;j++) // count all the butterflies in each block
+ {
+ par
+ {
+ // Weight factors for real (the same for FFT and iFFT)
+ weight1 = weight_re[((j-1)<<(level-1))<-7];
+
+
+ // Weight factors for imaginary (opposite for FFT and iFFT)
+ weight2 = (!select_inverse) ? (weight_im[((j-1)<<(level-1))<-7]) : -(weight_im[((j-1)<<(level-1))<-7]);
+
+ /* ORIGINAL CODE BELOW, MODIFIED BECAUSE OF MISMATCHING OUTPUT WITH BORLAND TESTAPP
+ weight2 = (!select_inverse) ? -(weight_im[((j-1)<<(level-1))<-7]) : weight_im[((j-1)<<(level-1))<-7];
+ */
+
+
+
+ for(i=0@j;i<=NUMBER_OF_POINTS;i+=e) // count all the blocks in this column
+ { // Butterfly calculation
+ par
+ {
+ point1 = ((i<-8)-1);
+ point2 = (((i<-8)+f)-1);
+ }
+
+ par
+ {
+ p=real.read[point1]+real.rwrite[point2];
+ q=imaginary.read[point1]+imaginary.rwrite[point2];
+ }
+
+ par
+ {
+ r=real.read[point1]-real.rwrite[point2];
+ t=imaginary.read[point1]-imaginary.rwrite[point2];
+ }
+
+ multiply(a,r,weight1);
+ multiply(b,t,weight2);
+
+ par
+ {
+ real.rwrite[point2] = ((a-b)>>FRACBITS)<-24;
+ imaginary.rwrite[point1] = q;
+ }
+
+ multiply(a,t,weight1);
+ multiply(b,r,weight2);
+
+ par
+ {
+ real.rwrite[point1] = p;
+ imaginary.rwrite[point2] = ((a+b)>>FRACBITS)<-24;
+ }
+
+ }
+ }
+ }
+ }
+
+ j=1;
+ for(i=1;i<NUMBER_OF_POINTS;i++)
+ {
+ if(i<(0@j))
+ {
+ par
+ {
+ point1=j-1;
+ point2=(i-1)<-8;
+ }
+ /*
+ COPYING ARRAY VALUES FROM ONE PLACE TO ANOTHER IN THE ARRAT MUST BE DONE IN
+ 2 STEPS. FIRSTLY THE VALUES ARE COPIED TO SEPARATE VARIABLES AFTER THAT THEY
+ ARE COPIED BACK TO THEIR NEW POSITION IN THE ARRAY. THIS MUST BE DONE TO
+ PREVENT TIMING ISSUES FROM OCCURING.
+ */
+ par
+ {
+ p = real.read[point1];
+ q = imaginary.read[point1];
+ }
+ par
+ {
+ r = real.read[point2];
+ t = imaginary.read[point2];
+ }
+ par
+ {
+ real.rwrite[point1] = r;
+ imaginary.rwrite[point1] = t;
+ }
+ par
+ {
+ real.rwrite[point2] = p;
+ imaginary.rwrite[point2] = q;
+ }
+ }
+
+ k = NUMBER_OF_POINTS>>1;
+
+
+ while(k<j)
+ {
+ j = j-k;
+ k = k>>1;
+ }
+
+ j+=k;
+ }
+
+}
+
+/********************************************************************/
+
+void perform_fft(signed 16 *pcm_audio)
+{
+ unsigned 8 k;
+ signed 24 sample;
+#if USE_UNSIGNED_AUDIO
+ unsigned 24 utemp;
+ signed 24 stemp;
+#endif
+
+ //initialize variables for the copying pipeline
+ k=0;
+ sample = adjs(pcm_audio[k],24);
+ // copy audio data to real-array before starting FFT calculation
+ // and set imaginary values to zero
+ do
+ {
+#if USE_UNSIGNED_AUDIO
+ stemp = adjs(pcm_audio[k],24);
+ stemp >>= 1;
+ utemp = 0@((unsigned)((stemp + 32768)<-16));
+ real.rwrite[k] = (signed)(utemp);
+#else
+ //Copying the array values has been pipelined to prevent parallel access to the
+ //pcm_audio array. This copying procedure must be finished before another
+ //sample is read from the audio input. The time available for this loop is
+ //determined by the sampling rate of 44,1 Khz
+ par
+ {
+ //COPYING NEEDS TO BE DONE IN 2 STEPS, BECAUSE THE VALUE THAT NEEDS TO WRITTEN
+ //TO THE REAL-RAM NEEDS TO BE AVAILABLE ON THE START OFF THE CLOCKCYCLE.
+ sample = adjs(pcm_audio[k+1],24);
+ real.rwrite[k] = sample;
+ imaginary.rwrite[k] = 0;
+ k++;
+ }
+
+#endif
+ } while (k);
+
+
+
+#if PERFORM_FFT_CALCULATION
+ calculate_fft(0);
+#endif
+
+
+}
+
+/********************************************************************/
+
+void perform_ifft(signed 16 *modified_audio/*, unsigned 6 *ifft_info*/)
+{
+ unsigned 6 k;
+ signed 24 p;
+
+#if PERFORM_FFT_CALCULATION
+ calculate_fft(1);
+#endif
+
+ k=0;
+#if PRINT_DEBUG
+ do
+ {
+ print_string("real[");
+ print_hex_value(0@k);
+ print_string("]: ");
+ print_hex_value((unsigned)real.read[k]);
+ print_string(" imaginary[");
+ print_hex_value(0@k);
+ print_string("]: ");
+ print_hex_value((unsigned)imaginary.read[k]);
+ print_eol();
+ k++;
+ }
+ while(k!=0);
+#endif
+//initialize variables for the copying pipeline
+#if PERFORM_FFT_CALCULATION
+ p = (real.read[(0@k)+95] >> NUMBER_OF_COLUMNS);
+#else
+ p = (real.read[(0@k)+95]);
+#endif
+
+ do
+ {
+ //Copying the array values has been pipelined to prevent parallel access to the
+ //pcm_audio array. This copying procedure must be finished before another
+ //sample is read from the audio input. The time available for this loop is
+ //determined by the sampling rate of 44,1 Khz
+ par
+ {
+#if PERFORM_FFT_CALCULATION
+ // divide samples by number of points
+ p = (real.read[(0@k+1)+95] >> NUMBER_OF_COLUMNS);
+#else
+ p = (real.read[(0@k+1)+95]);
+#endif
+#if USE_UNSIGNED_AUDIO
+ p -= 32768;
+ p <<= 1;
+#endif
+ modified_audio[k] = (p <-16);
+ //ifft_info[k] = (unsigned)(p[15:10]);
+ k++;
+ }
+ } while(k);
+}
+
+/********************************************************************/
+void equalize_audio(unsigned 4 *eq_level, unsigned 7 *fft_info)
+{
+#if 0
+ signed 24 p,q;
+ signed 16 a;
+ unsigned 8 i, mirror_i, bit, m, n;
+ unsigned 7 old_value;
+ unsigned 9 tmp;
+
+ //macro expr equalize_bar = adjs(q,24)-adjs(a,24);
+
+ macro proc equalize_bar()
+ {
+ signed result;
+ multiply(result, q,a);
+ return result[29:6];
+ }
+
+ p = real.read[0] - DC_COMPONENT; // remove DC component for calculations
+ real.rwrite[0] = p;
+// imaginary.rwrite[0] = 0; // remove DC component
+
+
+ for(i=0;i<=NUMBER_OF_FREQUENCIES;i++)
+ {
+ // set multiplication factor (0..64) for current frequency bar
+ a = adjs(eq_settings[eq_level[i<-7]],16);
+
+ // multiply frequency with this factor and divide by 64 (drop 6 LSB's)
+ q = real.read[i];
+ p = equalize_bar;
+ real.rwrite[i] = p;
+
+ q = imaginary.read[i];
+ p = equalize_bar;
+ imaginary.rwrite[i] = p;
+
+ // the upper part(128..255) of the spectrum is mirrored to the lower part;
+ // these values need to be adjusted too
+ if ((i<-7)!=0) // if not in DC component bar
+ {
+ mirror_i = (NUMBER_OF_POINTS-1)-i+1;
+ q = real.read[mirror_i];
+ p = equalize_bar;
+ real.rwrite[mirror_i] = p;
+
+ q = imaginary.read[mirror_i];
+ p = equalize_bar;
+ imaginary.rwrite[mirror_i] = p;
+ }
+ }
+
+ //write data to fft_info for display purposes
+ for(i=0;i<NUMBER_OF_FREQUENCIES;i++)
+ {
+ p = real.read[i];
+ q = imaginary.read[i];
+
+ if (p[23] == 1) p = -p; else delay;
+ if (q[23] == 1) q = -q; else delay;
+ p = (p<q) ? q : p; // This is done to get the best visual frequency result
+
+ if (display_log)
+ {
+
+ bit=126;
+ while ((p[21] == 0) && (bit != 0))
+ par
+ {
+ p = p<<1;
+ bit = bit - 18;
+ }
+ old_value = fft_info[i<-7];
+ tmp = ((0@old_value) + (0@bit))>>1;
+ fft_info[i<-7] = (old_value <= (tmp<-7)) ? (tmp<-7) : old_value-1;
+ }
+ else
+ {
+ old_value = fft_info[i<-7];
+ fft_info[i<-7] = (old_value<=(unsigned)(p[21:15])) ? (unsigned)(p[21:15]) : old_value-1;
+ }
+ }
+
+ // add DC component again before inverse FFT calculation is performed
+ p = real.read[0] + DC_COMPONENT;
+ real.rwrite[0] = p;
+#endif
+}
|