TI DSP的 FFT与IFFT库函数的使用测试

/**
 * Copyright (C) 2013 Guangzhou Tronlong Electronic Technology Co., Ltd. - www.tronlong.com
 *
 * @file dsplib_fft.c
 *
 * @brief Example application main file.
 * This application will run FFT performance test.
 *
 * @author Tronlong <[email protected]>
 *
 * @version V1.0
 *
 * @date 2019-09-12
 *
 **/

/* C Standard Header files */
#include <stdio.h>
#include <stdlib.h>

/* Compiler Header files */
#include <c6x.h>

/* CSL Header files */
#include <ti/csl/csl_cacheAux.h>

/* DSPLIB Header files */
#include <inc/dsplib.h>

/* MATHLIB Header files */
#include <inc/mathlib.h>

/* Local Header files */
#include "system/platform.h"

/* MAX length of FFT in complex samples */
#define MAXN    4096

/* Sampling Rate */
#define Fs      1000.0

#define PI      3.141592654

/* bit reverse table containing 64 entries  */
unsigned char brev[64] =
{
    0x00, 0x20, 0x10, 0x30, 0x08, 0x28, 0x18, 0x38,
    0x04, 0x24, 0x14, 0x34, 0x0c, 0x2c, 0x1c, 0x3c,
    0x02, 0x22, 0x12, 0x32, 0x0a, 0x2a, 0x1a, 0x3a,
    0x06, 0x26, 0x16, 0x36, 0x0e, 0x2e, 0x1e, 0x3e,
    0x01, 0x21, 0x11, 0x31, 0x09, 0x29, 0x19, 0x39,
    0x05, 0x25, 0x15, 0x35, 0x0d, 0x2d, 0x1d, 0x3d,
    0x03, 0x23, 0x13, 0x33, 0x0b, 0x2b, 0x1b, 0x3b,
    0x07, 0x27, 0x17, 0x37, 0x0f, 0x2f, 0x1f, 0x3f
};

/**
 * @details generate specialized sequence of twiddle factors for FFT
 *
 * @param w pointer to complex twiddle factor
 *
 * @param n length of FFT in complex samples
 *
 * @return NULL
 */
void tw_fft_gen(float *w, int n)
{
    int i, j, k;

    for(j = 1, k = 0; j <= n >> 2; j = j << 2) {
        for(i = 0; i < n >> 2; i += j) {
            w[k]     = (float)sinsp(2 * PI * i / n);
            w[k + 1] = (float)cossp(2 * PI * i / n);
            w[k + 2] = (float)sinsp(4 * PI * i / n);
            w[k + 3] = (float)cossp(4 * PI * i / n);
            w[k + 4] = (float)sinsp(6 * PI * i / n);
            w[k + 5] = (float)cossp(6 * PI * i / n);

            k += 6;
        }
    }

    return;
}

/**
 * @details generate specialized sequence of twiddle factors for IFFT
 *
 * @param w pointer to complex twiddle factor
 *
 * @param n length of FFT in complex samples
 *
 * @return NULL
 */
void tw_ifft_gen(float *w, int n)
{
    int i, j, k;

    for(j = 1, k = 0; j <= n >> 2; j = j << 2) {
        for(i = 0; i < n >> 2; i += j) {
            w[k]     = (float)-sinsp(2 * PI * i / n);
            w[k + 1] = (float) cossp(2 * PI * i / n);
            w[k + 2] = (float)-sinsp(4 * PI * i / n);
            w[k + 3] = (float) cossp(4 * PI * i / n);
            w[k + 4] = (float)-sinsp(6 * PI * i / n);
            w[k + 5] = (float) cossp(6 * PI * i / n);

            k += 6;
        }
    }

    return;
}

int main(void)
{
    uint32_t N = 0, i = 0;
    int32_t status = 0;
    float *fft_input, *fft_output, *ifft_output;
    float *signal_data, *fft_twiddle, *ifft_twiddle;
    uint32_t main_pll_freq, t_start, t_cost;
    uint8_t   radix;

    TSCH = 0;
    TSCL = 0;

    /* Get the cpu freq */
    main_pll_freq = platform_get_main_pll_freq();

    /* malloc buffer */
    fft_input = (float *)malloc(MAXN * sizeof(float) * 2);
    if(fft_input == NULL) {
        status = -1;
        printf("Failed to alloc fft_input memory !\r\n");
        goto err_return;
    }

    fft_output = (float *)malloc(MAXN * sizeof(float) * 2);
    if(fft_output == NULL) {
        status = -1;
        printf("Failed to alloc fft_output memory !\r\n");
        goto err_return;
    }

    ifft_output = (float *)malloc(MAXN * sizeof(float) * 2);
    if(ifft_output == NULL) {
        status = -1;
        printf("Failed to alloc ifft_output memory !\r\n");
        goto err_return;
    }

    fft_twiddle = (float *)malloc(MAXN * sizeof(float) * 2);
    if(fft_twiddle == NULL) {
        status = -1;
        printf("Failed to alloc fft_twiddle memory !\r\n");
        goto err_return;
    }

    ifft_twiddle = (float *)malloc(MAXN * sizeof(float) * 2);
    if(ifft_twiddle == NULL) {
        status = -1;
        printf("Failed to alloc ifft_twiddle memory !\r\n");
        goto err_return;
    }

    signal_data = (float *)malloc(MAXN * sizeof(float) * 2);
    if(signal_data == NULL) {
        status = -1;
        printf("Failed to alloc signal_data memory !\r\n");
        goto err_return;
    }

    /*
     * FFT input signal = signalA + signalB
     * signalA : sine wave signal with amplitude of 5 and frequency of 50 Hz
     * signalB : sine wave signal with amplitude of 15 and frequency of 150 Hz
     */
    for(i = 0; i < MAXN; i++) {
        signal_data[2 * i] = (float)5 * sinsp(2 * PI * 50 * (i / Fs)) + \
                           (float)15 * sinsp(2 * PI * 150 * (i / Fs));
        signal_data[2 * i + 1] = (float)0.0;
    }

    /* generate FFT twiddle factors */
    tw_fft_gen(fft_twiddle, MAXN);

    /* generate IFFT twiddle factors */
    tw_ifft_gen(ifft_twiddle, MAXN);

    /* set all L1 cache is 0, disable all cache */
    CACHE_setL1PSize(CACHE_L1_0KCACHE);
    CACHE_setL1DSize(CACHE_L1_0KCACHE);
    /* Stall CPU while memory system is busy */
    _mfence();
    _mfence();

    printf("--------------------------- Cache Disabled ---------------------------\n");
    /* between the 8 ~ MAXN, increase the FFT samples length by a multiple of 2 */
    for(N = 8; N <= MAXN; N = N * 2) {
        /* if N can be represented as Power of 4, radix=4, else radix=2 */
        if(N & 0x55555555) {
            radix = 4;
        } else {
            radix = 2;
        }

        /* Copy signal data to fft_input buffer */
        memcpy(fft_input, signal_data, (N * sizeof(float) * 2));
        /* start FFT conversion and calculate cost time */
        /* fast fourier transform - input data: fft_input, output data: fft_output */
        t_start = _itoll(TSCH, TSCL);
        DSPF_sp_fftSPxSP(N, fft_input, fft_twiddle, fft_output, brev, radix, 0, N);
        t_cost = _itoll(TSCH, TSCL) - t_start;

        printf("DSPF_sp_fftSPxSP   N = %5d| radix = %5d| cycles: %5d| time: %5.7f us\n", \
            N, radix, t_cost, t_cost/(main_pll_freq/1000000.0));

        /* start IFFT conversion and calculate cost time */
        /* inverse fast fourier transform - input data: fft_output, output data: ifft_output */
        t_start = _itoll(TSCH, TSCL);
        DSPF_sp_ifftSPxSP(N, fft_output, ifft_twiddle, ifft_output, brev, radix, 0, N);
        t_cost = _itoll(TSCH, TSCL) - t_start;

        printf("DSPF_sp_ifftSPxSP  N = %5d| radix = %5d| cycles: %5d| time: %5.7f us\n", \
            N, radix, t_cost, t_cost/(main_pll_freq/1000000.0));
    }

    /* enable all L1 cache */
    CACHE_setL1PSize(CACHE_L1_32KCACHE);
    CACHE_setL1DSize(CACHE_L1_32KCACHE);
    /* Stall CPU while memory system is busy */
    _mfence();
    _mfence();

    printf("--------------------------- Cache Enabled ----------------------------\n");
    /* between the 8 ~ MAXN, increase the FFT samples length by a multiple of 2 */
    for(N = 8; N <= MAXN; N = N * 2) {
        /* if N can be represented as Power of 4, radix=4, else radix=2 */
        if(N & 0x55555555) {
            radix = 4;
        } else {
            radix = 2;
        }

        /* Copy signal data to fft_input buffer */
        memcpy(fft_input, signal_data, (N * sizeof(float) * 2));
        /* start FFT conversion and calculate cost time */
        /* fast fourier transform - input data: fft_input, output data: fft_output */
        t_start = _itoll(TSCH, TSCL);
        DSPF_sp_fftSPxSP(N, fft_input, fft_twiddle, fft_output, brev, radix, 0, N);
        t_cost = _itoll(TSCH, TSCL) - t_start;

        printf("DSPF_sp_fftSPxSP   N = %5d| radix = %5d| cycles: %5d| time: %5.7f us\n", \
            N, radix, t_cost, t_cost/(main_pll_freq/1000000.0));

        /* start IFFT conversion and calculate cost time */
        /* inverse fast fourier transform - input data: fft_output, output data: ifft_output */
        t_start = _itoll(TSCH, TSCL);
        DSPF_sp_ifftSPxSP(N, fft_output, ifft_twiddle, ifft_output, brev, radix, 0, N);
        t_cost = _itoll(TSCH, TSCL) - t_start;

        printf("DSPF_sp_ifftSPxSP  N = %5d| radix = %5d| cycles: %5d| time: %5.7f us\n", \
            N, radix, t_cost, t_cost/(main_pll_freq/1000000.0));
    }

    /*
     * Dont care about the correctness of the calculation results,
     * loop FFT&IFFT, make the CPU and algorithm unit in a high load state
     */
    while(1) {
        DSPF_sp_fftSPxSP(N, fft_input, fft_twiddle, fft_output, brev, radix, 0, N);
        DSPF_sp_ifftSPxSP(N, fft_output, ifft_twiddle, ifft_output, brev, radix, 0, N);
    }

err_return:
    return status;
}