Compiler/TMDSEVM6678: how to build fftlib 2 unit tests/examples ? ... here is how I did

I have one question too ... why twiddles coefficients are computed as below in the fft_dp_1d_c2c_d.c code and also all other places. 

/* Function for generating Specialized sequence of twiddle factors */
void tw_gen_cn (double *w, int n)
{
    int i, j, k;
    const double PI = 3.141592654;

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

but the matlab code looks very different !!?? see below

% Program for generating n-length FFT's twiddle factor
% By: Denny Hermawanto
% Puslit Metrologi LIPI, INDONESIA
% Copyright 2015
function [real_twiddle, im_twiddle] = tw_gen_cn(fft_length)
    if nargin < 1
        [real_twiddle, im_twiddle] = tw_gen_cn(8);
        return;
    end   
% fft_length = input('Enter FFT length:');
for mm = 0:1:(fft_length-1)
    theta = (-2*pi*mm*1/fft_length);
%   Twiddle factor equation
%   twiddle = exp(1i*theta);
%   Euler equation for complex exponential 
%   e^(j*theta) = cos(theta) + j(sin(theta)) 
    
    twiddle(mm+1) = cos(theta) + (1i*(sin(theta)));
    real_twiddle = real(twiddle);
    real_twiddle = real_twiddle';
    im_twiddle = imag(twiddle);
    im_twiddle = im_twiddle';
    
end;
end

Mike,

Thanks for your contribution on how to compile the FFTLIB. This is helpful info and will benefit forum users.

Rex

ok, I found that the codes are mathematically equivalent ... look here

But now I have another question. From the code, I could see there are two ways to take FFT (1.Direct 2.ECPY); however the later doesn't seem to work:

1. What ECPY do exactly?

2. How to get it to work ?

Below you can see the code ... let me know if you any idea ...

Regards,

int main ()
{
    int     i, j, N, rad_cn;
    int     rad3 = 1, rad5 = 1, rad15 = 1;
    clock_t t_start, t_stop, t_overhead, t_opt;
    double  diff, max_diff = 0;
    fft_plan_t p;
    fft_callout_t plan_fxns;

    /* --------------------------------------------------------------------- */
    /* intialize hardware timers                                             */
    /* --------------------------------------------------------------------- */
    TSCL=0;TSCH=0;

    /* initalize callout functions */
    plan_fxns.memoryRequest   = fft_memory_request;
    plan_fxns.memoryRelease   = fft_memory_release;

    /* initialize ECPY */
    fft_assert( (EdmaMgr_init(DNUM, NULL) == EdmaMgr_SUCCESS), DNUM, "EdmaMgr_init() return error!");
    fftEdmaState = FFT_EDMA_STATE_INIT;

    /* radix 2&4 testing */
    for (N = 10; N <= 10; N = N*2)
    {
        memset (x_i,  0x55, sizeof (x_i) );
        memset (x_cn, 0x55, sizeof (x_cn));
        memset (w_i,  0x00, sizeof (w_i) );

        /* ---------------------------------------------------------------- */
        /* Initialize input vector temporarily.                             */
        /* ---------------------------------------------------------------- */

        for (i = 0; i < N; i++) {
          x_cn[PAD + 2*i  ] = i+1; // sin (2 * 3.1415 * 50 * i / (double) N);
          x_cn[PAD + 2*i+1] = 0; // sin (2 * 3.1415 * 60 * i / (double) N);
        }
        for (j = 0; j < 2*N; j++) {
          x_i[PAD + j] = x_cn[PAD + j];
        }

        /* ---------------------------------------------------------------- */
        /* Force uninitialized arrays to fixed values.                      */
        /* ---------------------------------------------------------------- */
        memset (y_i,  0xA5, sizeof (y_i) );
        memset (y_cn, 0xA5, sizeof (y_cn));

        /* ---------------------------------------------------------------- */
        /* Generate twiddle factors.                                        */
        /* ---------------------------------------------------------------- */
        j = 0;
        for (i = 0; i <= 31; i++)
            if ((N & (1 << i)) == 0)
                j++;
            else
                break;

        if (j % 2 == 0) {
            rad_cn = 4;
        }
        else {
            rad_cn = 2;
        }

        tw_gen_cn (ptr_w_cn, N);

        DSPF_dp_fftDPxDP (N, ptr_x_cn, ptr_w_cn, ptr_y_cn, rad_cn, 0, N);

        /* ---------------------------------------------------------------- */
        /* Compute the overhead of allocating and freeing EDMA              */
        /* ---------------------------------------------------------------- */
        p.edmaState = fft_assign_edma_resources();
        fft_free_edma_resources(p.edmaState);
        t_start = _itoll(TSCH, TSCL);
        p.edmaState = fft_assign_edma_resources();
        fft_free_edma_resources(p.edmaState);
        t_stop  = _itoll(TSCH, TSCL);
        t_overhead = t_stop - t_start;

        /***************************************
         * ecpy fft test
         ***************************************/
        plan_fxns.ecpyRequest = fft_assign_edma_resources;
        plan_fxns.ecpyRelease = fft_free_edma_resources;

        p = fft_dp_plan_1d_c2c (N, FFT_ECPY, plan_fxns);
        t_start = _itoll(TSCH, TSCL);
        fft_execute (p); // this doesn't work
        t_stop = _itoll(TSCH, TSCL);
        fft_destroy_plan (p);
        t_opt  = (t_stop - t_start) - t_overhead;

        /* ---------------------------------------------------------------- */
        /* compute difference and track max difference                      */
        /* ---------------------------------------------------------------- */
        diff = 0; max_diff = 0;
        for(i=0; i<2*N; i++) {
            printf("FFT_ECPY ptr_y_i[%d] = %f\n", i, ptr_y_i[i]);
          diff = _fabs(ptr_y_cn[i] - ptr_x_i[i]); // 
          if (diff > max_diff) max_diff = diff;
        }
        printf("fft_dp_1d_c2c_ecpy\tsize= %d\n", N);
        printf("max_diff = %f", max_diff);
        printf("\tN = %d\tCycle: %d\n\n", N, t_opt);

//        /********************************************************
//         * direct form fft test, up to 2048 due to twiddle size
//         *******************************************************/
        if (N <= 2048) {

          /* ---------------------------------------------------------------- */
          /* Initialize input vector temporarily.                             */
          /* ---------------------------------------------------------------- */

          for (i = 0; i < N; i++) {
            x_i[PAD + 2*i  ] = i+1; //sin (2 * 3.1415 * 50 * i / (double) N);
            x_i[PAD + 2*i+1] = 0; sin (2 * 3.1415 * 60 * i / (double) N);
          }

          /* ---------------------------------------------------------------- */
          /* Compute the overhead of EDMA assigment                           */
          /* ---------------------------------------------------------------- */
          t_start = _itoll(TSCH, TSCL);
          t_stop  = _itoll(TSCH, TSCL);
          t_overhead = t_stop - t_start;

          plan_fxns.ecpyRequest = NULL;
          plan_fxns.ecpyRelease = NULL;


          p = fft_dp_plan_1d_c2c (N, FFT_DIRECT, plan_fxns);
          t_start = _itoll(TSCH, TSCL);
          fft_execute (p); // this works
          t_stop = _itoll(TSCH, TSCL);
          fft_destroy_plan (p);
          t_opt  = (t_stop - t_start) - t_overhead;


          /* ---------------------------------------------------------------- */
          /* compute difference and track max difference                      */
          /* ---------------------------------------------------------------- */
          diff = 0; max_diff = 0;
          for(i=0; i<2*N; i++) {
              printf("FFT_DIRECT ptr_y_i[%d] = %f\n", i, ptr_y_i[i]);
            diff = _fabs(ptr_y_cn[i] - ptr_y_i[i]);
            if (diff > max_diff)
              max_diff = diff;
          }

          printf("fft_dp_1d_c2c_direct\tsize= %d\n", N);
          printf("max_diff = %f", max_diff);
          printf("\tN = %d\tCycle: %d\n\n", N, t_opt);

          for(i=0; i<2*N; i+=2) {
                printf("FFT_DIRECT fft[%d] = %f + i(%f) \n", i+1, ptr_y_i[i], ptr_y_i[i+1]);
          }
        }
    }

}

Rex Chang said:

Mike,

Thanks for your contribution on how to compile the FFTLIB. This is helpful info and will benefit forum users.

Rex

No problem ... That's all we do hear :)