Running FFT on external signal

Hi guys, I'm pretty inexperienced in this but I'm trying to compute FFT on an electrical signal. Could you tell me what I might be doing wrong with my code.

Methodology:
The code creates an array where it stores the samples from the ADC conversion and then uses them for the FFT. Thereafter I do a simple multiplication.

Below are the 2 code files. 1 is the main code, and the other is related to the famous fix_ftt (really long) code for the FFT algorithm done in 1989

#include "fix_fft.h"
#include <math.h> 
#include <arduino.h>

#define LOG2N      10 //base2 log of N [N being number of points and length of spectrum]
#define NPTS       1024  // Number of points or length of spectrum in points
#define FS       2000 // 2kHz is the sampling frequency

int real[NPTS], imag[NPTS];
void sampleAdc(void); 
void startFFT(void); 

void setup() 
{
   Serial.begin(9600);   //baud rate and data transfer speed set
  delay(1000);          // delay to give serial monitor tme to start up 
  Serial.println("Start"); 
}
 
void loop() 
{ 
  sampleAdc();  //calls the sampleADC function below which does the ADC 
  for (int i = 0; i < NPTS; i++) //this loop initializes 1024 samples and prints each as an element of a real array
  {
     Serial.println(real[i]);  
  }

  startFFT(); 

}

void sampleAdc (void)
{
  int i = 0;
  int ADC_CH;
  while (i < NPTS)
  {   //read ADC pin NPTS times at hi
       
          real[i] = analogRead(A0); 
          i++; 
       
  }
   
} //end of void sampleAdc


//do fft to the sampled signal 
//output shows peak frequency 

void startFFT (void) //Computes FFT of sampled signal,identifies peak frequency, calculates approximate top speed and outputs it 
{
  // absolute
  int i;    
 short fix_fft(real, imag, LOG2N, 0);    //performs fft on sampled points
  for(i = 0; i < NPTS/2; i++)       //get the power amplitude in 
    { 
      real[1] = sqrt((int)real[i] *(int)real[i]+ (int) imag[i] * (int)imag[i]);
    }  
    //find the peak 
     int peakHz = 0;
     int peaki = 0;
     for (i = 1; i < NPTS/2; i++)  //bin 0 holds the summation
     {  
        if (real[i] > peakHz) 
           { 
              peakHz = real[i];
              peaki = i;
            } 
     }
    peaki = peaki * FS/NPTS;
    int topSpeed = peaki/70; 
    Serial.print ("Top speed  : ");
    Serial.println (topSpeed);
    Serial.println (""); 

} // end of FFT function
#ifdef fix_fft_h
#define fix_fft_h

#include <arduino.h>

int fix_fft(short fr[], short fi[], short m, short inverse);
int fix_fftr(short f[], int m, int inverse); 

#endif

Could you tell me what I might be doing wrong with my code.

Probably not, until you explain what you expect to be the result, and what the actual result is.

It is ALWAYS best to test unfamiliar code with an example for which you know the answer.

In this case, prepare a sample that consists of a single sine wave of known frequency, and then run the FFT. The result should be a single peak of the correct amplitude at the correct frequency.

I expect it to print out "Top Speed is ____"


So I was initially building this with an MSP430, and so the code was a bit different. But as I was struggling with configuring the ADC I took the unfinished code and started editing it to work for an Arduino.

I am away from the signal generator/electric signal for another few hours so I haven't properly tested the code. I was hoping you could see in my use of analogRead and calling of the functions if there was something amiss.


That being said:

When I compile now it says

Done Compiling

[C:\directory: In function 'void startFFT()':

[C:\directory: warning: expression list treated as compound expression in initializer [-fpermissive]

short fix_fft(real, imag, LOG2N, 0);    //performs fft on sampled points

I'm not seeing where the "imag" array gets assigned anything sensible (zeroed?) before computing the FFT.

I don't know what that warning means. Because you haven't posted all the code, I can't reproduce the message.

I am away from the signal generator/electric signal

Again: DO NOT use a signal generator, until you have tested the code with a set of calculated input values.

DO NOT expose the ADC to negative voltages, or you may destroy it.

Also you must not allow external signals of frequency greater than 1/2 the ADC sample frequency to reach the input. Otherwise "aliasing artifacts" can make the output incomprehensible.

Example test function (untested):

void sampleAdc (void)
{
  int i = 0;
  float twopi = 2.0*3.14159;
  float freq = 16.; //example signal frequency

  int ADC_CH;
  while (i < NPTS)
  {   //create fake ADC signal NPTS
       
          real[i] = 100.0*sin(twopi*freq*i/NPTS)+200; //positive only sine wave signal
          imag[i] = 0; //do not forget to set imaginary part to zero
          i++;
       
  }
   
} //end of void sampleAdc

Which Arduino is this running on? I presume (hope) it isn't UNO/NANO because they don't have enough ram for the 1000 samples.
Sampling frequency is going to be closer to 9kHz (assuming this is a 16MHz Arduino) because the sampleAdc function samples as fast as it can.
When you get a sample it will be a positive number in the range 0 - 1023. The FFT expects signed integers so you have subtract the DC offset - usually by doing this:

real[i] = analogRead(A0) - 512;

Have you offset your signal so that it is in the range 0 to 5V?

As @jremington has requested, post all your code.

Pete

(Thank you for replies so far)

The rest of the code ; a third code file named fix_fft.cpp

[I have to edit the int "Sinewave[N_WAVE-N_WAVE/4]" function because there were lots of numbers there and I couldn't post the full code because of the 9000 character limit.]

/* fix_fft.c - Fixed-point in-place Fast Fourier Transform  */
*/

#define N_WAVE      1024    /* full length of Sinewave[] */
#define LOG2_N_WAVE 10      /* log2(N_WAVE) */

/*
  Henceforth "short" implies 16-bit word. If this is not
  the case in your architecture, please replace "short"
  with a type definition which *is* a 16-bit word.
*/

/*
  Since we only use 3/4 of N_WAVE, we define only
  this many samples, in order to conserve data space.
*/
int Sinewave[N_WAVE-N_WAVE/4] = {
      0,    201,    402,  to 32767, then back to 0, and then steps down to -32766,
};

/*
  FIX_MPY() - fixed-point multiplication & scaling.
  Substitute inline assembly for hardware-specific
  optimization suited to a particluar DSP processor.
  Scaling ensures that result remains 16-bit.
*/
inline int FIX_MPY(int a, int b)
{
  /* shift right one less bit (i.e. 15-1) */
  int c = ((int)a * (int)b) >> 14;
  /* last bit shifted out = rounding-bit */
b = c & 0x01;
  /* last shift + rounding bit */
  a = (c >> 1) + b;
  return a;
}

/*
  fix_fft() - perform forward/inverse fast Fourier transform.
  fr[n],fi[n] are real and imaginary arrays, both INPUT AND
  RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to
  0 for forward transform (FFT), or 1 for iFFT.
*/
int fix_fft(int fr[], int fi[], int m, int inverse)
{
  int mr, nn, i, j, l, k, istep, n, scale, shift;
  int qr, qi, tr, ti, wr, wi;

  n = 1 << m;

  /* max FFT size = N_WAVE */
  if (n > N_WAVE)
    return -1;

  mr = 0;
  nn = n - 1;
  scale = 0;

  /* decimation in time - re-order data */
  for (m=1; m<=nn; ++m) {
    l = n;
    do {
      l >>= 1;
    } while (mr+l > nn);
    mr = (mr & (l-1)) + l;

    if (mr <= m)
      continue;
    tr = fr[m];
    fr[m] = fr[mr];
    fr[mr] = tr;
    ti = fi[m];
    fi[m] = fi[mr];
    fi[mr] = ti;
  }

  l = 1;
  k = LOG2_N_WAVE-1;
  while (l < n) {
    if (inverse) {
      /* variable scaling, depending upon data */
      shift = 0;
      for (i=0; i<n; ++i) {
        j = fr[i];
        if (j < 0)
          j = -j;
        m = fi[i];
        if (m < 0)
          m = -m;
        if (j > 16383 || m > 16383) {
          shift = 1;
          break;
        }
      }
      if (shift)
        ++scale;
    } else {
      /*
        fixed scaling, for proper normalization --
        there will be log2(n) passes, so this results
        in an overall factor of 1/n, distributed to
        maximize arithmetic accuracy.
      */
      shift = 1;
    }
    /*
      it may not be obvious, but the shift will be
      performed on each data point exactly once,
      during this pass.
    */
    istep = l << 1;
    for (m=0; m<l; ++m) {
      j = m << k;
      /* 0 <= j < N_WAVE/2 */
      wr =  Sinewave[j+N_WAVE/4];
      wi = -Sinewave[j];
      if (inverse)
        wi = -wi;
      if (shift) {
        wr >>= 1;
        wi >>= 1;
      }
      for (i=m; i<n; i+=istep) {
        j = i + l;
        tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
        ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
        qr = fr[i];
        qi = fi[i];
        if (shift) {
          qr >>= 1;
          qi >>= 1;
        }
        fr[j] = qr - tr;
        fi[j] = qi - ti;
        fr[i] = qr + tr;
        fi[i] = qi + ti;
      }
    }
    --k;
    l = istep;
  }
  return scale;
}


int fix_fftr(int f[], int m, int inverse)
{
  int i, N = 1<<(m-1), scale = 0;
  int tt, *fr=f, *fi=&f[N];

  if (inverse)
    scale = fix_fft(fi, fr, m-1, inverse);
  for (i=1; i<N; i+=2) {
    tt = f[N+i-1];
    f[N+i-1] = f[i];
    f[i] = tt;
  }
  if (! inverse)
    scale = fix_fft(fi, fr, m-1, inverse);
  return scale;
}

===========

@jremington :

I put your sampleADC() function into the code and tried it. The top speed printed is always "0" even when I change the frequency to higher numbers (e.g 210).

Also, what are you suggesting by saying I shouldn't expose Arduino to negative voltages? Do I have to rectify the signal first? Or is that something I should be able to do on the signal generator?

========

@el_supremo:

I'm using an Arduino Mega 2560.

Offset the signal? I haven't returned the to the signal generator, but thanks for pointing that out.

Also, what are you suggesting by saying I shouldn't expose Arduino to negative voltages?

We are suggesting that you will destroy the Arduino if you do so. I thought that was stated rather clearly.

The top speed printed is always "0" even when I change the frequency to higher numbers (e.g 210).

Print out the entire result of the FFT instead.

This project is too difficult for a beginner. We recommend that you start with and work your way through the simple examples that come with the Arduino development software, in order to learn the programming language and the special features of the Arduino.

I understood that. The follow up question explained what I was trying to say: Are you suggesting that I simply put a DC offset at the signal generator or AC rectify the signal first? [I assume the former but just wanted to see]

[Sorry I can't really back down from this project now]

The result of the startFFT function appears to be a bin number where the highest magnitude is stored. However it seems this bin's values were assigned to it just by simply dividing the sine wave spectrum. i.e. 0 would be 0, but the peak magnitude at 1 Hz is the amplitude of the test sine wave at number/bin 256 [1024/4, which makes sense to me].

It doesn't seem any FFT is going on.

I realized it was giving 0 previously because it's set to 'int' and I was getting a decimal less than 0 because I was dividing by 70 (I've taken this out for the moment)

Are you suggesting that I simply put a DC offset at the signal generator or AC rectify the signal first

You definitely should not rectify the signal. That changes its spectrum.
You need to DC offset the signal so that it is entirely within the range 0 - 5V (it must not be negative) and it is centred around 2.5V. This puts the whole signal within the range of the ADC. It is also why you then have to subtract 512 from each sample to remove the DC bias that you've added.
If it isn't centred around 2.5V (e.g. maybe you have a 1Vp-p signal centred on 2V) you have to figure out the value that's required to get the correct offset value to subtract (about 410 for a 2V centre, 1024*2/5).

Pete