Non-synthesisable VHDL code for 8 point FFT algorithm

A Fast Fourier Transform(FFT) is an efficient algorithm for calculating the discrete Fourier transform of a set of data. A DFT basically decomposes a set of data in time domain into different frequency components. DFT is defined by the following equation:

   A FFT algorithm uses some interesting properties of the above formula to simply the calculations. You can read more about these FFT algorithms here.
   Many students have been asking doubts regarding vhdl implementation of FFT, so I decided to write a sample code. I have selected 8 point decimation in time(DIT) FFT algorithm for this purpose. In short I have wrote the code for this flow diagram of FFT.
   This is just  a sample code, which means it is not synthesisable. I have used real data type for the inputs and outputs and all calculations are done using the math_real library. The inputs can be complex numbers too.
   To define the basic arithmetic operations between two complex numbers I have defined some new functions which are available in the package named fft_pkg. The component named, butterfly , contains the basic butterfly calculations for FFT as shown in this flow diagram.
    I wont be going any deep into the theory behind FFT here. Please visit the link given above or Google  for in depth theory. There are 4 vhdl codes in the design,including the testbench code, and are given below.

The package file - fft_pkg.vhd:

library IEEE;
use IEEE.std_logic_1164.all;
use IEEE.MATH_REAL.ALL;

package fft_pkg is

type complex is
    record
        r : real;
        i : real;
    end record;

type comp_array is array (0 to 7) of complex;
type comp_array2 is array (0 to 3) of complex;

function add (n1,n2 : complex) return complex;
function sub (n1,n2 : complex) return complex;
function mult (n1,n2 : complex) return complex;

end fft_pkg;

package body fft_pkg is

--addition of complex numbers
function add (n1,n2 : complex) return complex is

variable sum : complex;

begin
sum.r:=n1.r + n2.r;
sum.i:=n1.i + n2.i;
return sum;
end add;

--subtraction of complex numbers.
function sub (n1,n2 : complex) return complex is

variable diff : complex;

begin
diff.r:=n1.r - n2.r;
diff.i:=n1.i - n2.i;
return diff;
end sub;

--multiplication of complex numbers.
function mult (n1,n2 : complex) return complex is

variable prod : complex;

begin
prod.r:=(n1.r * n2.r) - (n1.i * n2.i);
prod.i:=(n1.r * n2.i) + (n1.i * n2.r);
return prod;
end mult;

end fft_pkg;

The top level entity - fft8.vhd:

library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
use IEEE.MATH_REAL.ALL;
library work;
use work.fft_pkg.ALL;

entity fft8 is
port(   s : in comp_array; --input signals in time domain
        y : out comp_array  --output signals in frequency domain
        );
end fft8;

architecture Behavioral of fft8 is

component butterfly is
   port(
      s1,s2 : in complex;      --inputs
      w :in complex;      -- phase factor
      g1,g2 :out complex      -- outputs
   );
end component; 
   
signal g1,g2 : comp_array := (others => (0.0,0.0));
--phase factor, W_N = e^(-j*2*pi/N)  and N=8 here.
--W_N^i = cos(2*pi*i/N) - j*sin(2*pi*i/N);  and i has range from 0 to 7.
constant w : comp_array2 := ( (1.0,0.0), (0.7071,-0.7071), (0.0,-1.0), (-0.7071,-0.7071) );

begin

--first stage of butterfly's.
bf11 : butterfly port map(s(0),s(4),w(0),g1(0),g1(1));
bf12 : butterfly port map(s(2),s(6),w(0),g1(2),g1(3));
bf13 : butterfly port map(s(1),s(5),w(0),g1(4),g1(5));
bf14 : butterfly port map(s(3),s(7),w(0),g1(6),g1(7));

--second stage of butterfly's.
bf21 : butterfly port map(g1(0),g1(2),w(0),g2(0),g2(2));
bf22 : butterfly port map(g1(1),g1(3),w(2),g2(1),g2(3));
bf23 : butterfly port map(g1(4),g1(6),w(0),g2(4),g2(6));
bf24 : butterfly port map(g1(5),g1(7),w(2),g2(5),g2(7));

--third stage of butterfly's.
bf31 : butterfly port map(g2(0),g2(4),w(0),y(0),y(4));
bf32 : butterfly port map(g2(1),g2(5),w(1),y(1),y(5));
bf33 : butterfly port map(g2(2),g2(6),w(2),y(2),y(6));
bf34 : butterfly port map(g2(3),g2(7),w(3),y(3),y(7));
   
end Behavioral;

Butterfly component - butterfly.vhd:

library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
library work;
use work.fft_pkg.ALL;

entity butterfly is
   port(
      s1,s2 : in complex;      --inputs
      w :in complex;      -- phase factor
      g1,g2 :out complex      -- outputs
   );
end butterfly;

architecture Behavioral of butterfly is

begin

--butterfly equations.
g1 <= add(s1,mult(s2,w));
g2 <= sub(s1,mult(s2,w));

end Behavioral;

Testbench code - tb_fft8.vhd:

LIBRARY ieee;
USE ieee.std_logic_1164.ALL;
library work;
use work.fft_pkg.all;

ENTITY tb IS
END tb;

ARCHITECTURE behavior OF tb IS
   signal s,y : comp_array;

BEGIN

    -- Instantiate the Unit Under Test (UUT)
   uut: entity work.fft8 PORT MAP (
          s => s,
             y => y
        );
   
   -- Stimulus process
   stim_proc: process
   begin       
    --sample inputs in time domain.
      s(0) <= (-2.0,1.2);  
        s(1) <= (-2.2,1.7);
        s(2) <= (1.0,-2.0);
        s(3) <= (-3.0,-3.2);   
        s(4) <= (4.5,-2.5);
        s(5) <= (-1.6,0.2);
        s(6) <= (0.5,1.5); 
        s(7) <= (-2.8,-4.2);   
      wait;
   end process;

END;

Copy and paste the above codes into their respective files and simulate. You will get the output in the output signal 'y'. 

NOTE: I have used real data types in this code. 'real' type is not synthesisable. You cannot convert the real type into a hardware equivalent circuit. so if we have to synthesis the code we have to use fixed point arithmetic for the data types..