Discrete Fast Fourier Transform

In Verilog

- Decompose frequency from sound

1. Extracting frequencies from sound

2. Uncertainty principle

3. Riemann Zeta function and primes

4. Differential equations

Part 1. Sampling

    /* 
    Part 1: Volume Sampling
        Sampling frequency; 20kHz
        Number of samples; 10000
        Period: 0.5 seconds
    */
    // initialisation of the array to 12'b0
    initial begin
        freq = 12'd0;
        for (c = 0; c < 14'd10000;c = c + 1) begin
            SAMPLE[c] = 12'd0;
        end
    end
    //sample volume at 20kHz frequency
    reg onecycle = 1'b0;
    always @(posedge MCLK) begin
        if (counter == 14'd9999) begin
        onecycle = 1'b1;
            if(onecycle == 1'b1) begin
                FREQ = (RESULT_1000>RESULT_2000 && RESULT_1000>RESULT_3000) ? 3'b001 : (RESULT_2000>RESULT_3000)? 3'b011:3'b111;
                counter <= 14'b00000000000000;
                onecycle <= 1'b0;
            end
        //reset counter value after the calculations are done
        end
        else begin
            counter <= counter + 1;    
            SAMPLE[counter] = mic_in;   //g(t)
            if(vart_1000==5'd20) begin
                    vart_1000<= 5'b00001;
                end
            else vart_1000 <= vart_1000 + 1'b1;
            if(vart_2000==5'd10) begin
                    vart_2000<= 5'b00001;
            end
            else vart_2000 <= vart_2000 + 1'b1;
            if(vart_3000==5'd7) begin
                    vart_3000<= 5'b00001;
            end
            else vart_3000 <= vart_3000 + 1'b1;
        end
        
    end
    

Part 2. Calculation

Functions Used

• Look-Up Table Method

module FFT_sin (input [90:0] angle, output [12:0] SIN_value);//Divide by 4096
    reg [12:0] SIN_TABLE[90:0];
    always @(*) begin
        SIN_TABLE = {13'd0,13'd71,13'd142,13'd214,13'd285,13'd356,13'd428,13'd499,13'd570,13'd640,
                    13'd711,13'd781,13'd851,13'd921,13'd990,13'd1060,13'd1128,13'd1197,13'd1265,13'd1333,
                    13'd1400,13'd1467,13'd1534,13'd1600,13'd1665,13'd1731,13'd1795,13'd1859,13'd1922,13'd1985,
                    13'd2047,13'd2109,13'd2170,13'd2230,13'd2290,13'd2349,13'd2407,13'd2465,13'd2521,13'd2577,
                    13'd2632,13'd2687,13'd2740,13'd2793,13'd2845,13'd2896,13'd2946,13'd2995,13'd3043,13'd3091,
                    13'd3137,13'd3183,13'd3227,13'd3271,13'd3313,13'd3355,13'd3395,13'd3435,13'd3473,13'd3510,
                    13'd3547,13'd3582,13'd3616,13'd3649,13'd3681,13'd3712,13'd3741,13'd3770,13'd3797,13'd3823,
                    13'd3848,13'd3872,13'd3895,13'd3917,13'd3937,13'd3956,13'd3974,13'd3991,13'd4006,13'd4020,
                    13'd4033,13'd4045,13'd4056,13'd4065,13'd4073,13'd4080,13'd4086,13'd4090,13'd4093,13'd4095,13'd4096};
    end
    assign SIN_value = SIN_TABLE [angle];
endmodule

module FFT_cos (input [90:0] angle, output [12:0] COS_value);
    reg [12:0] COS_TABLE[90:0];
    always @(*) begin
        COS_TABLE = {13'd4096,13'd4095,13'd4093,13'd4090,13'd4086,13'd4080,13'd4073,13'd4065,13'd4056,13'd4045,
                    13'd4033,13'd4020,13'd4006,13'd3991,13'd3974,13'd3956,13'd3937,13'd3917,13'd3895,13'd3872,
                    13'd3848,13'd3823,13'd3797,13'd3770,13'd3741,13'd3712,13'd3681,13'd3649,13'd3616,13'd3582,
                    13'd3547,13'd3510,13'd3473,13'd3435,13'd3395,13'd3355,13'd3313,13'd3271,13'd3227,13'd3183,
                    13'd3137,13'd3091,13'd3043,13'd2995,13'd2946,13'd2896,13'd2845,13'd2793,13'd2740,13'd2687,
                    13'd2632,13'd2577,13'd2521,13'd2465,13'd2407,13'd2349,13'd2290,13'd2230,13'd2170,13'd2109,
                    13'd2048,13'd1985,13'd1922,13'd1859,13'd1795,13'd1731,13'd1666,13'd1600,13'd1534,13'd1467,
                    13'd1400,13'd1333,13'd1265,13'd1197,13'd1129,13'd1060,13'd990,13'd921,13'd851,13'd781,
                    13'd711,13'd640,13'd570,13'd499,13'd428,13'd357,13'd285,13'd214,13'd143,13'd71,13'd0};
    end
    assign COS_value = COS_TABLE [angle];
endmodule

Example: Calculation for 1000Hz

// Section 1: 1000Hz
    FFT_sin sin_1000(angle_90_1000, SIN_value_1000);
    FFT_cos cos_1000(angle_90_1000, COS_value_1000);
    always @(posedge MCLK) begin
        //angle = (13'd6280 * vart_1000 / 15'd20000) * 9'd360 / 6.28;  // (0~6.28) need to improve the calculation
        angle_1000 = 19'd360000 * vart_1000; //(0~360)
        
        if (angle_1000 < 8'd180) begin   //final value is positive
            if (angle_1000 <= 7'd90) begin
                angle_90_1000 = angle_1000;
                IMAG_pos_1000 <=  IMAG_pos_1000 + mic_in * SIN_value_1000;
            end
            else begin
                angle_90_1000 = angle_1000 - 7'd90; //sin£¨¦Ð/2£«¦Á£©£½cos¦Á
                IMAG_pos_1000 <= IMAG_pos_1000 + mic_in * COS_value_1000;
            end
        end
        else begin //final value is negative  //sin£¨¦Ð£«¦Á£©£½£­sin¦Á
            angle_180_1000 = angle_1000 - 9'd180;
            if (angle_180_1000 <= 7'd90) begin
                angle_90_1000 = angle_1000;
                IMAG_neg_1000 <=  IMAG_neg_1000 + mic_in * SIN_value_1000;
            end
            else begin
                angle_90_1000 = angle_1000 - 7'd90; //sin£¨¦Ð/2£«¦Á£©£½cos¦Á
                IMAG_neg_1000 <= IMAG_neg_1000 + mic_in * COS_value_1000;
          end
        end
        
        if(vart_1000==5'd20) begin
            final_cal_1000 = final_cal_1000 + 1;
            REAL_1000 = REAL_pos_1000 - REAL_neg_1000;
            IMAG_1000 = IMAG_pos_1000 - IMAG_neg_1000;
        end
        if (final_cal_1000 == 2'b01) begin
            RESULT_1000 <= RESULT_1000 + IMAG_1000*IMAG_1000 + REAL_1000*REAL_1000;
            final_cal_1000 = final_cal_1000 + 1;
        end
        if (final_cal_1000 -- 2'b10) final_cal_1000 <= 2'b00;
    end