Problem with FFT implementation via api

[Program]

Hi, there! I've decided to write my own program to implement signal acquisition and FFT of the signal using the api and kissfft. I slightly modified some of the functions (inculded in acq_fpga_fft.c below) so that a signal is only acquired on one channel. I'm only really interested in which frequency has the highest magnitude in the FFT, not so much about exact units (dB...) of the magnitudes. Signal acquisition works fine, as far as I can tell, but the FFT is giving me incorrect results. I've based my code on the way the kissfft routine is called in the worker function of the acquirefft program of edgo (see here). After I've called the acq_fpga_fft() function (see code below), I've basically got what I want, right? At a decimation of 1, my frequency bin size is (sampling frequency/2)/(FFT length), so for 125 MHz / 2 / 8192 = 7.6 kHz. However, feeding in a signal of 10MHz and plotting frequency against the output of my FFT program, I get a peak at 45 MHz. Any ideas where the error comes from?

main program:

Code: Select all

//This is a program to acquire samples, compute FFT and output magnitude of Fourier coefficients

//#include <iostream>
#include <stdio.h>
#include <sys/time.h>
#include <unistd.h>
#include <sys/types.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include <sys/param.h>

#include "redpitaya/rp.h"
#include "main.h"
#include "acq_FPGA_FFT.h"
#include "kiss_fftr.h"

double time_diff(struct timeval x , struct timeval y);
double s_mean(double *arr, uint32_t size);

int main()
{
	if(rp_Init() != 0){
		fprintf(stderr, "RP API init failed!\n");
	}
	
                       //BUFFFER_SIZE is number of samples to acquire (for now 16384), FFT length is the length of the output of the kissfft routine, which is
                       //(BUFFER_SIZE/2), both are defined in header files
	uint32_t buff_size = BUFFER_SIZE;
	float *buff = (float *)malloc(buff_size * sizeof(float));
	double *buff_worker = (double *)malloc(buff_size * sizeof(double));
	double *buff_fft = (double *)malloc(sizeof(double) * FFT_LENGTH);
	
	rp_acq_trig_src_t trigger_source;
	
	struct timeval start, now;
	double microseconds;
	
	gettimeofday(&start, NULL);
	rp_AcqReset();
	gettimeofday(&now, NULL);
	microseconds = time_diff(start, now);
	//printf("Time to initialise oscilloscope: %f\n", microseconds);
	rp_AcqSetDecimation(1); //Possible values: 1, 8, 64, 1024, 8192, 65536), corresponding
				//sampling rate = (125 MSps / decimation)
        rp_AcqSetTriggerDelay(buff_size); //There is already a zero-offset delay of 8192 samples
        rp_AcqSetGain(RP_CH_1, RP_LOW);
        
        //Starts acquisition
        rp_AcqStart();
        
        gettimeofday(&start, NULL);
        rp_AcqSetTriggerSrc(RP_TRIG_SRC_NOW);
        
        while(1){
        	rp_AcqGetTriggerSrc(&trigger_source);
        	if(trigger_source == RP_TRIG_SRC_DISABLED){
        		break;
        	}
        	usleep(50);
        }
        gettimeofday(&now, NULL);
	microseconds = time_diff(start, now);
	//printf("Time to acquire data: %f\n", microseconds);
	
        rp_AcqGetLatestDataV(RP_CH_1, &buff_size, buff);
        
        int i;
        //All other functions take array of type double, so convert
        for(i = 0; i < buff_size; i++){
        	buff_worker[i] = (double)buff[i];
        }
        
        //Subtract DC offset as I'm only interested in frequency content
        double dc_offset = s_mean(buff_worker, buff_size);
        for(i = 0; i < buff_size; i++){
        	buff_worker[i] = buff_worker[i] - dc_offset;
        }
        
        /*printf("Voltage samples:\n");
        for(i = 0; i < buff_size; i++){
        	printf("%f\n", buff[i]);
        }*/
        //printf("DC offset: %f\n", dc_offset);
        
        //Hanning filter
        if(acq_fpga_hann_init() < 0){
        	acq_fpga_hann_clean();
        	free(buff);
        	rp_Release();
        	return -1;
        }
        
        if(acq_fpga_fft_init() < 0){
        	acq_fpga_fft_clean();
        	free(buff);
        	rp_Release();
        	return -1;
        }
        //acq_fpga_hann_filter(&buff_worker[0], &buff_worker);
        acq_fpga_fft(&buff_worker[0], (double **)&buff_fft);
        printf("Signal acquired, FFT computed: \n");

        for(i = 0; i < FFT_LENGTH; i++){
        	buff_fft[i] = buff_fft[i] * 1e10; //scaling for plotting, not really necessary
        	printf("%7f\n", buff_fft[i]);
        }
        
        acq_fpga_hann_clean();
        acq_fpga_fft_clean();
        free(buff);
        free(buff_worker);
        free(buff_fft);
        rp_Release();
        
        return 0;
}

double time_diff(struct timeval x , struct timeval y)
{    
	double x_ms , y_ms , diff;
	
	x_ms = (double)x.tv_sec*1000000 + (double)x.tv_usec;
	y_ms = (double)y.tv_sec*1000000 + (double)y.tv_usec;
	diff = (double)y_ms - (double)x_ms;

	return diff;
}

double s_mean(double *arr, uint32_t size)
{
	double sum = 0.0;
	double mean;
	
	int i;
	for(i = 0; i < size; i++){
		sum += arr[i];
	}
	
	mean = sum / size;
	
	return mean;
}

acq_FPGA_FFT.c:

Code: Select all

//This file contains the functions used in main.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#include "acq_FPGA_FFT.h"
#include "kiss_fftr.h"
#include "main.h"

kiss_fft_cpx *rp_kiss_fft_out = NULL;
kiss_fftr_cfg rp_kiss_fft_cfg  = NULL;
double *rp_hann_window   = NULL;

int acq_fpga_hann_init()
{
    int i;

    acq_fpga_hann_clean(rp_hann_window);

    rp_hann_window = (double *)malloc(BUFFER_SIZE * sizeof(double));
    if(rp_hann_window == NULL) {
        fprintf(stderr, "rp_spectr_hann_create() can not allocate mem");
        return -1;
    }
    
    for(i = 0; i < BUFFER_SIZE; i++) {
        rp_hann_window[i] = RP_SPECTR_HANN_AMP * 
            (1 - cos(2*M_PI*i / (double)(BUFFER_SIZE-1)));
    }

    return 0;
}

int acq_fpga_hann_clean()
{
    if(rp_hann_window){
        free(rp_hann_window);
        rp_hann_window = NULL;
    }
    return 0;
}

int acq_fpga_hann_filter(double *ch_in, double **ch_out)
{
    int i;
    double *ch_o = *ch_out;

    if(!ch_in || !*ch_out){
        return -1;
    }
    
    for(i = 0; i < BUFFER_SIZE; i++) {
        ch_o[i] = ch_in[i] * rp_hann_window[i];
    }

    return 0;
}

int acq_fpga_fft_init()
{
    if(rp_kiss_fft_out || rp_kiss_fft_cfg){
        acq_fpga_fft_clean();
    }

    rp_kiss_fft_out = 
        (kiss_fft_cpx *)malloc(BUFFER_SIZE * sizeof(kiss_fft_cpx));

    rp_kiss_fft_cfg = kiss_fftr_alloc(BUFFER_SIZE, 0, NULL, NULL);

    return 0;
}

int acq_fpga_fft_clean()
{
    kiss_fft_cleanup();
    if(rp_kiss_fft_out){
        free(rp_kiss_fft_out);
        rp_kiss_fft_out = NULL;
    }
    if(rp_kiss_fft_cfg){
        free(rp_kiss_fft_cfg);
        rp_kiss_fft_cfg = NULL;
    }
    
    return 0;
}

int acq_fpga_fft(double *ch_in, double **ch_out)
{
	double *ch_o = *ch_out;
	int i;
	
	if(!ch_in || !*ch_out){
		return -1;
	}
	
	if(!rp_kiss_fft_out || !rp_kiss_fft_cfg){
	 	 fprintf(stderr, "acq_fpga_fft not initialized");
	 	 return -1;
        }
	
	kiss_fftr(rp_kiss_fft_cfg, (kiss_fft_scalar *)ch_in, rp_kiss_fft_out);
	
	for(i = 0; i < FFT_LENGTH; i++){                     // FFT limited to fs/2, specter of amplitudes
		ch_o[i] = sqrt(pow(rp_kiss_fft_out[i].r, 2) + 
                               pow(rp_kiss_fft_out[i].i, 2));
        }
        
        return 0;
}

[Program]

Never mind, I found the error: In my code above, the function call to set the decimation has to read AcqSetDecimation(RP_DEC_1) instead of AcqSetDecimation(1). I copied this from the acquire_trigger_software.c file that is contained in C/Examples on the redpitaya. This needs to be corrected in the redpitaya github repository, too.
By the way, why did the compiler not throw an error, because 1 is not of type rp_acq_decimation_t?

[Nils Roos]

By the way, why did the compiler not throw an error, because 1 is not of type rp_acq_decimation_t?

In C, enum types are considered compatible with int, even when you give a number that does not correspond to a defined enum value (which you didn't, in this case).