TMS320C6678: How to make DSPLIP DSP_fft32x32 function output same as matlab?

Feng P.

Part Number: TMS320C6678

I have been trying to use DSP_fft32x32 I have looked into the example under C:\ti\dsplib_c66x_3_4_0_2\examples\fft_ex\fft_example.c;

I just can't get the same results as matlab by try an error I noticed that the input samples should change between -1 and 1 and then scale by 2147483648 and then shifting right as below:

The output needs to be scaled back as far as I experienced thought not too sure what I am doing ...

1) I searched the forum before posting the question and all says it's been answered previously but I can't find anything relevant. Please let me know if you any experience.

2) I also found DSPF_dp_fftDPxDP but I have no clue how it would perform better than DSP_fft32x32? I have a hard time understanding "DSPF_sp_fftSPxSP_cn.c is the natural C equivalent of the optimized linear assembly code without restrictions". :-)

I'd appreciate your valuable comments.

Regards

over 5 years ago

0 Rex Chang over 5 years ago

TI__Guru 50170 points

Hi, Feng,

The expert on FFT of the DSPLib is currently out of office. He will get on the issue after he is back in office.Thanks for your patience.

Rex

0 Rahul Prabhu over 5 years ago

TI__Guru** 116180 points

DSPLIB fixed point functions have a long legacy since they were orginally written for C64x+ DSP devices and then maintained in the C66x version. For understanding scaling and overflow handling, you will need to refer to the documentation on C64x+ version of the lib. I have attached the application notes below:

Please refer to section 2.2.5 below and for implementation details refer to page 49-50. Hope this helps.

Regards,

Rahul

0 Rahul Prabhu over 5 years ago in reply to Rahul Prabhu

TI__Guru** 116180 points

Also for fixed point FFT with no data scaling refer to section 2 in the application notes below:

www.ti.com/.../spracn4.pdf

0 Feng P. over 5 years ago in reply to Rahul Prabhu

Intellectual 570 points

Thank you for details I read through and I may post new questions if any.

Appreciate it.

0 Feng P. over 5 years ago in reply to Rahul Prabhu

Intellectual 570 points

Hi Rahul,

Thanks a lot for the tip.

Apparently I succeeded ... so input signal has to be between 1 and -1 plus it has to be scaled down (divided to avoid overflow) by 2^log2(nx) where nx is number of the samples in the input signal; and then the output of the DSP_fft32x32 should be scaled back to (multiplied by) 2^log2(nx).

Below is how I did it .. I hope it helps others.

/*
 * DSP_fft32x32_calc.c
 *
 *  Created on: Aug 17, 2020
 *      Author: Feng-laptop
 */
#include <stdio.h>

#include <math.h>
#include <ti/dsplib/dsplib.h>

#include "DSP_fft32x32_calc.h"
#include "gen_twiddle_fft32x32.h"
#include "utils.h"
#include "ark_memory.h"

int DSP_fft32x32_calc(Complex* img, dt_int nn){
    printf("hellow ...DSP_fft32x32_calc \n");
//    printline
//    ASSERT(round(log2(nn)) == log2(nn), "Number of input samples is power of 2");
    printline
    // allocate mem
    dt_int32 *w_32x32 = (dt_int32*) ark_malloc(2*nn*sizeof(dt_int32));
    dt_int32 *x_32x32 = (dt_int32*) ark_malloc(2*nn*sizeof(dt_int32));
    dt_int32 *y_32x32 = (dt_int32*) ark_malloc(2*nn*sizeof(dt_int32));

    printline
    ASSERT(w_32x32 != NULL && x_32x32 != NULL && y_32x32 != NULL, "DSP_fft32x32_calc: cannot allocte memory" )
    // convert input data to fixed-point
    printline
    dt_double max_x = x_32x32[0];
    for(int i = 0; i < nn; i++){
        double im = abs(CIMAG(img[i]));
        double rl = abs(CREAL(img[i]));

        if (im>max_x){
            max_x = im;
        }
        if (rl>max_x){
            max_x = rl;
        }
    }

    if (max_x <= 1 ){
        max_x = 1;
        printf("max is set to 1.0 \n");
    } else {
//        max_x *= 1.25;
    }

    printf("max_x = %f\n", max_x);
    int scale = log2(nn);
//    int scale = 3;
#define FIX2FLOAT_SCALE_INPUT(x) (((dt_int)(x / max_x * 2147483648)) >> scale)
    for(int i = 0; i < nn; i++){
        x_32x32[2*i] = FIX2FLOAT_SCALE_INPUT(CREAL(img[i]));
        x_32x32[2*i + 1] = FIX2FLOAT_SCALE_INPUT(CIMAG(img[i]));
//        printf("x_32x32[ %d ] = (%d,%d) \n", i, x_32x32[ 2*i ], x_32x32[ 2*i + 1]);
    }

    printline
    // gen twiddles
    gen_twiddle_fft32x32(w_32x32, nn, 2147483647.5);

    printline
    DSP_fft32x32(w_32x32, nn, x_32x32, y_32x32);

//    printf("\n\n");
    printline
    // convert back fft result to ouput as floating point number
#define FIX2FLOAT_SCALE_OUTPUT(x) (pow(2,scale)*(float)(x) * max_x / 2147483648.f)
    for(int i = 0; i < nn; i++){
//        printf("y_32x32[ %d ] = (%d,%d) \n", i, y_32x32[ 2*i ], y_32x32[ 2*i + 1]);

        img[i] = FIX2FLOAT_SCALE_OUTPUT(y_32x32[ 2*i ]) + FIX2FLOAT_SCALE_OUTPUT(y_32x32[ 2*i + 1 ])*I;
    }

//    printf("\n\n");
//    for(int i = 0; i < nn; i++){
//        printf("img[%d] = (%f,%f)\n", i, CREAL(img[i]), CIMAG(img[i]));
//    }

    printline
    // release mem
    ark_free(y_32x32, 2*nn*sizeof(dt_int32));
    ark_free(x_32x32, 2*nn*sizeof(dt_int32));
    ark_free(w_32x32, 2*nn*sizeof(dt_int32));

    printline
    return 0;
}

Don't forget to multiply input by 2147483647.5 as 32 bit fixed point calculation takes up 4 bytes. and then divide the output by 2147483647.5.

Regards,

0 Feng P. over 5 years ago in reply to Rahul Prabhu

Intellectual 570 points

hi Rahul,

Thank so much for the documents provided ... I was successful calculating DSP_fft32x32 but the inverse function DSP_ifft32x32 doesn't work using the below mentioned scaling by the document provided. I did exact same scaling as I did for DSP_fft32x32 but the same scaling doesn't seem work for the inverse function ifft.

Below is the code I use please let me know your valuable comments.

Much appreciate it.

/*
 * DSP_fft32x32_calc.c
 *
 *  Created on: Aug 17, 2020
 *      Author: Feng-laptop
 */
#include <stdio.h>

#include <math.h>
#include <ti/dsplib/dsplib.h>

#include "DSP_fft32x32_calc.h"
#include "gen_twiddle_fft32x32.h"
#include "utils.h"
#include "ark_memory.h"

int DSP_fft32x32_calc(Complex* img, dt_int nn){
    printf("hellow ...DSP_fft32x32_calc \n");
//    printline
//    ASSERT(round(log2(nn)) == log2(nn), "Number of input samples is power of 2");
    printline
    // allocate mem
    dt_int32 *w_32x32 = (dt_int32*) ark_malloc(2*nn*sizeof(dt_int32));
    dt_int32 *x_32x32 = (dt_int32*) ark_malloc(2*nn*sizeof(dt_int32));
    dt_int32 *y_32x32 = (dt_int32*) ark_malloc(2*nn*sizeof(dt_int32));

    printline
    ASSERT(w_32x32 != NULL && x_32x32 != NULL && y_32x32 != NULL, "DSP_fft32x32_calc: cannot allocte memory" )
    // convert input data to fixed-point
    printline
    dt_double max_x = x_32x32[0];
    for(int i = 0; i < nn; i++){
        double im = abs(CIMAG(img[i]));
        double rl = abs(CREAL(img[i]));

        if (im>max_x){
            max_x = im;
        }
        if (rl>max_x){
            max_x = rl;
        }
    }

    if (max_x <= 1 ){
        max_x = 1;
        printf("max is set to 1.0 \n");
    } else {
//        max_x *= 1.25;
    }

    printf("max_x = %f\n", max_x);
    int scale = log2(nn);
//    int scale = 3;
#define FIX2FLOAT_SCALE_INPUT(x) (((dt_int)(x / max_x * 2147483648)) >> scale)
    for(int i = 0; i < nn; i++){
        x_32x32[2*i] = FIX2FLOAT_SCALE_INPUT(CREAL(img[i]));
        x_32x32[2*i + 1] = FIX2FLOAT_SCALE_INPUT(CIMAG(img[i]));
//        printf("x_32x32[ %d ] = (%d,%d) \n", i, x_32x32[ 2*i ], x_32x32[ 2*i + 1]);
    }

    printline
    // gen twiddles
    gen_twiddle_fft32x32(w_32x32, nn, 2147483647.5);

    printline
    DSP_fft32x32(w_32x32, nn, x_32x32, y_32x32);

//    printf("\n\n");
    printline
    // convert back fft result to ouput as floating point number
#define FIX2FLOAT_SCALE_OUTPUT(x) (pow(2,scale)*(float)(x) * max_x / 2147483648.f)
    for(int i = 0; i < nn; i++){
//        printf("y_32x32[ %d ] = (%d,%d) \n", i, y_32x32[ 2*i ], y_32x32[ 2*i + 1]);

        img[i] = FIX2FLOAT_SCALE_OUTPUT(y_32x32[ 2*i ]) + FIX2FLOAT_SCALE_OUTPUT(y_32x32[ 2*i + 1 ])*I;
    }

//    printf("\n\n");
//    for(int i = 0; i < nn; i++){
//        printf("img[%d] = (%f,%f)\n", i, CREAL(img[i]), CIMAG(img[i]));
//    }

    printline
    // release mem
    ark_free(y_32x32, 2*nn*sizeof(dt_int32));
    ark_free(x_32x32, 2*nn*sizeof(dt_int32));
    ark_free(w_32x32, 2*nn*sizeof(dt_int32));

    printline
    return 0;
}

Processors

Processors forum

TMS320C6678: How to make DSPLIP DSP_fft32x32 function output same as matlab?