Loop not running when using timer interrupts

This is my first project using lower level functionality like timers and interrupts and I am struggling to understand the problem with my code. I hope someone can help.

I am making a remote controlled car timing system. Each car will carry an infrared (IR) emitter flashing at a unique frequency (modulated at 38khz) in order to identify itself. By the side of the track will be an IR receiver which will sample any incoming data at a steady rate and run it through a fast fourier transform (FFT) to identify frequencies present in the signal.

The emitters seem to be working perfectly, I can monitor the frequency on an oscilloscope and I have also identified the frequency on my receiver with code that simply samples some data in some set amount of time and then runs the FFT calculation before looping.

So far so good.

I now want to do this whole process at a certain rate so that it is useful as a race timer. I can perform the FFT calculation in 7 or 8 milliseconds so I believe I can run the timer at a resolution of 0.01 seconds. That is to say that I can gather samples over a 0.1 second period and then in the following 0.1 second period perform the FFT calculation whilst using interrupts to collect the next set of samples.

This is where I am having problems. I thought I understood how interrupts worked but I cannot explain the output of my code. In the following example I have slowed the timer down by a factor of 1024 to ensure that I am not running out of time to perform the FFT. I would expect to see the following over serial output.

Processing
(Number)
Processing
(Number)
Processing
(Number)
Processing
(Number)
...

However I see the following

Processing
(Number)
Processing

And then it stops outputting. The first time it processes, the FFT calculation is done immediately (single digit milliseconds as mentioned above) but it is as though the second one never completes. I cannot come up with a reason for this. If I log the counter inside the interrupt (I know this isn't something you should do but I believe that it is ok for debugging this problem as I have slowed down the process immensely) then I can see it continuing to count and to reset the canProcess flag but the loop never seems to read it again.

Have I missed something to do with interrupts? Or something else?

The FFT library is this one, I have copied the code from the first snippet block near the bottom of that page. I have left it out of my code here for simplicity. The only change I made was to lower the accuracy of the fast_sine funtion from 5 to 2 to ensure it is performing comfortably under 10ms.

My code:

int IRPin = 9;

volatile int dataIndex = 0;
volatile int dataOne[64] = {};
volatile int dataTwo[64] = {};
volatile bool fillingDataOne = true;
volatile bool canProcess = false;

void setup() {
  pinMode(IRPin, INPUT);

  cli();
  TCCR1A = 0;
  TCCR1B = 0;
  // TCCR1B |= _BV(CS10); // No prescaler
  TCCR1B |= _BV(CS10) | _BV(CS12); // Prescaler of 1024 for debugging purposes only
  TCCR1B |= _BV(WGM12); // Count up to OCR1A
  OCR1A = 2499; // 6400Hz (64 samples every 10 milliseconds)
  TCNT1 = 0; // Reset counter
  TIMSK1 |= _BV(OCIE1A); // Interrupt
  sei();

  Serial.begin(115200);
}

ISR (TIMER1_COMPA_vect) {
  if (fillingDataOne) {
    dataOne[dataIndex++] = digitalRead(IRPin);
  } else {
    dataTwo[dataIndex++] = digitalRead(IRPin);
  }
  if (dataIndex == 64) {
    dataIndex = 0;
    canProcess = true;
    fillingDataOne = !fillingDataOne;
  }
}

void loop() {
  if (canProcess) {
    Serial.println("Processing");
    canProcess = false;
    float f = Approx_FFT(fillingDataOne ? dataTwo : dataOne, 64, 6400);
    Serial.println(f);
  }
}

Any help is greatly appreciated :pray:

I'll leave it to the heavies, but it seems complicated.

Do you expect to be able to detect more than one car, more than one frequency modulated carrier, at a time?

I'd ditch the FFT and just measure the frequency by looking at the edges of the IR receiver output.

TV remote control reception and decoding does not usually involve FFT, I'm not sure I see the point.

Either you will get a pulse train from the demodulation or you won't.

Maybe you are using a different kind of receiver that is more analog-y, tell ~us~ me.

a7

You code doesn't include the library header, it just won't compile...
Please post the whole code.

PS
Also I didn't see where you started the timer in your code.

I think this code is non-functional... please show exact the code that you tested.

I'm so impressed by how quickly people reply on here :+1:

Yes, I expect to detect more than one frequency at the same time and that is the only reason for the FFT. I have a (mostly) working timer based on car color and analysing the input from a camera so I know its pretty common to have more than one car passing the line at the same time (especially at the start of a race).

As I mentioned, the FFT functionality seems to work really well. In fact, here is the code that works flawlessly (apart from the fact that I'm not always collecting data or indeed at the right sample rate therefore it's not going to work as a race timer)...

int IRPin = 9;

volatile int dataIndex = 0;
volatile int data[64] = {};

void setup() {
  pinMode(IRPin, INPUT);
  Serial.begin(115200);
}

void loop() {
  if (dataIndex == 64) {
    Serial.println("Processing");
    float f = Approx_FFT(data, 64, 100);
    Serial.println(f);
    dataIndex = 0;
  } else {
    data[dataIndex++] = digitalRead(IRPin);
    delay(1);
  }
}

OK, cool.

Can you say what transmitters are on the cars and what receiver is stationed by the track?

Make and model, or if you are rolling your own, schematics.

I'm only familiar with parts for that what just would not work.

TIA

a7

Hi b707, thanks for the reply.

The full code is below, the library was just pasted straight into my code.

Also I didn't see where you started the timer in your code.

I'm not sure what you mean by this. Everything between cli() and sei() in setup is timer code. It triggers the ISR function directly below the setup function.

I can confirm the code compiles, uploads and runs - I see serial output.

int IRPin = 9;

volatile int dataIndex = 0;
volatile int dataOne[64] = {};
volatile int dataTwo[64] = {};
volatile bool fillingDataOne = true;
volatile bool canProcess = false;

void setup() {
  pinMode(IRPin, INPUT);

  cli();
  TCCR1A = 0;
  TCCR1B = 0;
  // TCCR1B |= _BV(CS10); // No prescaler
  TCCR1B |= _BV(CS10) | _BV(CS12); // Prescaler of 1024 for debugging purposes only
  TCCR1B |= _BV(WGM12); // Count up to OCR1A
  OCR1A = 2499; // 6400Hz (64 samples every 10 milliseconds)
  TCNT1 = 0; // Reset counter
  TIMSK1 |= _BV(OCIE1A); // Interrupt
  sei();

  Serial.begin(115200);
}

ISR (TIMER1_COMPA_vect) {
  if (fillingDataOne) {
    dataOne[dataIndex++] = digitalRead(IRPin);
  } else {
    dataTwo[dataIndex++] = digitalRead(IRPin);
  }
  if (dataIndex == 64) {
    dataIndex = 0;
    canProcess = true;
    fillingDataOne = !fillingDataOne;
  }
}

void loop() {
  if (canProcess) {
    Serial.println("Processing");
    canProcess = false;
    float f = Approx_FFT(fillingDataOne ? dataTwo : dataOne, 64, 6400);
    Serial.println(f);
  }
}

// Everything below here is from Approx_FTT. I have passed it through a formatter to fix indentations

//---------------------------------lookup data------------------------------------//
byte isin_data[128] = {
  0,
  1,
  3,
  4,
  5,
  6,
  8,
  9,
  10,
  11,
  13,
  14,
  15,
  17,
  18,
  19,
  20,
  22,
  23,
  24,
  26,
  27,
  28,
  29,
  31,
  32,
  33,
  35,
  36,
  37,
  39,
  40,
  41,
  42,
  44,
  45,
  46,
  48,
  49,
  50,
  52,
  53,
  54,
  56,
  57,
  59,
  60,
  61,
  63,
  64,
  65,
  67,
  68,
  70,
  71,
  72,
  74,
  75,
  77,
  78,
  80,
  81,
  82,
  84,
  85,
  87,
  88,
  90,
  91,
  93,
  94,
  96,
  97,
  99,
  100,
  102,
  104,
  105,
  107,
  108,
  110,
  112,
  113,
  115,
  117,
  118,
  120,
  122,
  124,
  125,
  127,
  129,
  131,
  133,
  134,
  136,
  138,
  140,
  142,
  144,
  146,
  148,
  150,
  152,
  155,
  157,
  159,
  161,
  164,
  166,
  169,
  171,
  174,
  176,
  179,
  182,
  185,
  188,
  191,
  195,
  198,
  202,
  206,
  210,
  215,
  221,
  227,
  236
};
unsigned int Pow2[14] = {
  1,
  2,
  4,
  8,
  16,
  32,
  64,
  128,
  256,
  512,
  1024,
  2048,
  4096
};
byte RSSdata[20] = {
  7,
  6,
  6,
  5,
  5,
  5,
  4,
  4,
  4,
  4,
  3,
  3,
  3,
  3,
  3,
  3,
  3,
  2,
  2,
  2
};
//---------------------------------------------------------------------------------//

//-----------------------------FFT Function----------------------------------------------//
/*
Code   to perform High speed and Accurate FFT on arduino,
setup:

1. in[]     :   Data array, 
2. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
3.   Frequency: sampling frequency required as input (Hz)

It will by default return   frequency with max aplitude,
if you need complex output or magnitudes uncomment   required sections

If sample size is not in power of 2 it will be clipped   to lower side of number. 
i.e, for 150 number of samples, code will consider   first 128 sample, remaining sample  will be omitted.
For Arduino nano, FFT of   more than 256 sample not possible due to mamory limitation 
Code by ABHILASH
Contact:   abhilashpatel121@gmail.com
Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/

Update(06/05/21):   Correction made for support on Arduino Due
*/

float Approx_FFT(int in [], int N, float Frequency) {
  int a, c1, f, o, x, data_max, data_min = 0;
  long data_avg, data_mag, temp11;
  byte scale, check = 0;

  data_max = 0;
  data_avg = 0;
  data_min = 0;

  for (int i = 0; i < 12; i++) //calculating the levels
  {
    if (Pow2[i] <= N) {
      o = i;
    }
  }
  a = Pow2[o];
  int out_r[a]; //real part of transform
  int out_im[a]; //imaginory part of transform

  for (int i = 0; i < a; i++) //getting   min max and average for scalling
  {
    out_r[i] = 0;
    out_im[i] = 0;
    data_avg = data_avg + in [i];
    if (in [i] > data_max) {
      data_max = in [i];
    }
    if (in [i] < data_min) {
      data_min = in [i];
    }
  }

  data_avg = data_avg >> o;
  scale = 0;
  data_mag = data_max - data_min;
  temp11 = data_mag;

  //scalling data  from +512 to -512

  if (data_mag > 1024) {
    while (temp11 > 1024) {
      temp11 = temp11 >> 1;
      scale = scale + 1;
    }
  }

  if (data_mag < 1024) {
    while (temp11 < 1024) {
      temp11 = temp11 << 1;
      scale = scale + 1;
    }
  }

  if (data_mag > 1024) {
    for (int i = 0; i < a; i++) {
      in [i] = in [i] - data_avg;
      in [i] = in [i] >> scale;
    }
    scale = 128 - scale;
  }

  if (data_mag < 1024) {
    scale = scale - 1;
    for (int i = 0; i < a; i++) {
      in [i] = in [i] - data_avg;
      in [i] = in [i] << scale;
    }

    scale = 128 + scale;
  }

  x = 0;
  for (int b = 0; b < o; b++) // bit reversal order stored   in im_out array
  {
    c1 = Pow2[b];
    f = Pow2[o] / (c1 + c1);
    for (int j = 0; j < c1; j++) {
      x = x + 1;
      out_im[x] = out_im[j] + f;
    }
  }

  for (int i = 0; i < a; i++) // update input array as per bit reverse order
  {
    out_r[i] = in [out_im[i]];
    out_im[i] = 0;
  }

  int i10, i11, n1, tr, ti;
  float e;
  int c, s, temp4;
  for (int i = 0; i < o; i++) //fft
  {
    i10 = Pow2[i]; // overall values of sine/cosine  
    i11 = Pow2[o] / Pow2[i + 1]; // loop with similar sine cosine
    e = 1024 / Pow2[i + 1]; //1024 is equivalent   to 360 deg
    e = 0 - e;
    n1 = 0;

    for (int j = 0; j < i10; j++) {
      c = e * j; //c is angle as where 1024 unit is 360 deg
      while (c < 0) {
        c = c + 1024;
      }
      while (c > 1024) {
        c = c - 1024;
      }

      n1 = j;

      for (int k = 0; k < i11; k++) {
        temp4 = i10 + n1;
        if (c == 0) {
          tr = out_r[temp4];
          ti = out_im[temp4];
        } else if (c == 256) {
          tr = -out_im[temp4];
          ti = out_r[temp4];
        } else if (c == 512) {
          tr = -out_r[temp4];
          ti = -out_im[temp4];
        } else if (c == 768) {
          tr = out_im[temp4];
          ti = -out_r[temp4];
        } else if (c == 1024) {
          tr = out_r[temp4];
          ti = out_im[temp4];
        } else {
          tr = fast_cosine(out_r[temp4], c) - fast_sine(out_im[temp4], c); //the   fast sine/cosine function gives direct (approx) output for A*sinx
          ti = fast_sine(out_r[temp4], c) + fast_cosine(out_im[temp4], c);
        }

        out_r[n1 + i10] = out_r[n1] - tr;
        out_r[n1] = out_r[n1] + tr;
        if (out_r[n1] > 15000 || out_r[n1] < -15000) {
          check = 1;
        } //check for int size, it can handle only +31000   to -31000,

        out_im[n1 + i10] = out_im[n1] - ti;
        out_im[n1] = out_im[n1] + ti;
        if (out_im[n1] > 15000 || out_im[n1] < -15000) {
          check = 1;
        }

        n1 = n1 + i10 + i10;
      }
    }

    if (check == 1) { //   scalling the matrics if value higher than 15000 to prevent varible from overflowing
      for (int i = 0; i < a; i++) {
        out_r[i] = out_r[i] >> 1;
        out_im[i] = out_im[i] >> 1;
      }
      check = 0;
      scale = scale - 1; //   tracking overall scalling of input data
    }

  }

  if (scale > 128) {
    scale = scale - 128;
    for (int i = 0; i < a; i++) {
      out_r[i] = out_r[i] >> scale;
      out_im[i] = out_im[i] >> scale;
    }
    scale = 0;
  } // revers   all scalling we done till here,
  else {
    scale = 128 - scale;
  } //   in case of nnumber getting higher than 32000, we will represent in as multiple of   2^scale

  /*
  for(int i=0;i<a;i++)
  {
  Serial.print(out_r[i]);Serial.print("\  ");                     // un comment to print RAW o/p    
  Serial.print(out_im[i]);   
  Serial.print("i");Serial.print("\ "); 
  Serial.print("*2^");Serial.println(scale);   
  }
  */

  //---> here onward out_r contains amplitude and our_in conntains   frequency (Hz)
  int fout, fm, fstp;
  float fstep;
  fstep = Frequency / N;
  fstp = fstep;
  fout = 0;
  fm = 0;

  for (int i = 1; i < Pow2[o - 1]; i++) // getting amplitude from compex   number
  {
    out_r[i] = fastRSS(out_r[i], out_im[i]);
    // Approx RSS function used to calculated magnitude quickly

    out_im[i] = out_im[i - 1] + fstp;
    if (fout < out_r[i]) {
      fm = i;
      fout = out_r[i];
    }

    // un comment to   print Amplitudes (1st value (offset) is not printed)
    //Serial.print(out_r[i]);   Serial.print("\  "); 
    //Serial.print("*2^");Serial.println(scale);   
  }

  float fa, fb, fc;
  fa = out_r[fm - 1];
  fb = out_r[fm];
  fc = out_r[fm + 1];
  fstep = (fa * (fm - 1) + fb * fm + fc * (fm + 1)) / (fa + fb + fc);

  return (fstep * Frequency / N);
}

//---------------------------------fast   sine/cosine---------------------------------------//

int fast_sine(int Amp, int th) {
  int temp3, m1, m2;
  byte temp1, temp2, test, quad, accuracy;
  accuracy = 2; // set it value from 1 to 7, where 7 being most accurate but slowest
  //   accuracy value of 5 recommended for typical applicaiton
  while (th > 1024) {
    th = th - 1024;
  } // here 1024 = 2*pi or 360 deg
  while (th < 0) {
    th = th + 1024;
  }
  quad = th >> 8;

  if (quad == 1) {
    th = 512 - th;
  } else if (quad == 2) {
    th = th - 512;
  } else if (quad == 3) {
    th = 1024 - th;
  }

  temp1 = 0;
  temp2 = 128; //2 multiple
  m1 = 0;
  m2 = Amp;

  temp3 = (m1 + m2) >> 1;
  Amp = temp3;
  for (int i = 0; i < accuracy; i++) {
    test = (temp1 + temp2) >> 1;
    temp3 = temp3 >> 1;
    if (th > isin_data[test]) {
      temp1 = test;
      Amp = Amp + temp3;
      m1 = Amp;
    } else if (th < isin_data[test]) {
      temp2 = test;
      Amp = Amp - temp3;
      m2 = Amp;
    }
  }

  if (quad == 2) {
    Amp = 0 - Amp;
  } else if (quad == 3) {
    Amp = 0 - Amp;
  }
  return (Amp);
}

int fast_cosine(int Amp, int th) {
  th = 256 - th; //cos th = sin (90-th) formula
  return (fast_sine(Amp, th));
}

//--------------------------------------------------------------------------------//

//--------------------------------Fast   RSS----------------------------------------//
int fastRSS(int a, int b) {
  if (a == 0 && b == 0) {
    return (0);
  }
  int min, max, temp1, temp2;
  byte clevel;
  if (a < 0) {
    a = -a;
  }
  if (b < 0) {
    b = -b;
  }
  clevel = 0;
  if (a > b) {
    max = a;
    min = b;
  } else {
    max = b;
    min = a;
  }

  if (max > (min + min + min)) {
    return max;
  } else {
    temp1 = min >> 3;
    if (temp1 == 0) {
      temp1 = 1;
    }
    temp2 = min;
    while (temp2 < max) {
      temp2 = temp2 + temp1;
      clevel = clevel + 1;
    }
    temp2 = RSSdata[clevel];
    temp1 = temp1 >> 1;
    for (int i = 0; i < temp2; i++) {
      max = max + temp1;
    }
    return (max);
  }
}
//--------------------------------------------------------------------------------//

Hi @alto777. Apologies, I was blocked by the spam filter.

It's all self built using two Arduino Pro Micros (5v). Here is the schematic, please excuse my doodling...

Emitter and receiver hardware...

Thanks

Yes, the code between cli() and sei() configures timer registers. But after configure the timer you should to start it by setting the lowest bit of TCC1B to one. Where are you doing it in your code?

With this sketch your timer just will never start.

OK, I understand what you're saying now.

I'm still confused though - assuming you're referring to TCCR1B (unless TCC1B is something I have missed?) then the lowest bit is CS10 I believe. In the below guide it is used for the prescaler. In fact I am setting it in my debug code for the prescaler of 1024. Have I misunderstood something?

I am certain that the timer runs - I can serial log the count.

https://medesign.seas.upenn.edu/index.php/Guides/MaEvArM-timer1

Edit:

I am just realising that I didn't mention which arduino I was using up front - it's a Pro Micro 5v (ATMega32u4). This could lead to confusion.

NP, that's actually helpful.

You state that FFT is working, have you tried placing two cars, each modulating its 38 kHz carrier but at a distinct frequency in front of the receiver, and applying FFT to the signal you get?

Forget for the moment the fancy footwork so you can be gathering the next samples whilst processing a batch you have already.

Does that show energy at the two frequencies?

I'm going for a proof of concept test here. I remain to be convinced that two digital signals can be received by that kind of setup and become total trash no matter how you are processing the one result digital signal.

I get that if you had a microphone and two whistlers you'd expect success. But analog, superposition and linearity and stuff.

Like I said differently, above my pay grade.

a7

Sorry, I didn't notice that you set this bit in the prescaler. Then the objection is removed.
Although, in general, the datasheet recommends setting all the timer registers when the counter is stopped, that is, the start of the timer should be the last setting item.
Changing registers while the timer is running can lead to undefined behavior

@b707 Ok cool, that's good to know. Does the cli() / sei() process not start and stop the timers essentially? Anyway you've given me something to learn about, thanks.

@alto777 No I haven't tried two cars yet and I appreciate that it could all fall down there. I don't believe it will be a problem though, I think this is exactly what FFT is for.

However, for the purposes of this post I don't think it matters at all. If I end up getting junk data that's OK and I can go in a different direction. I posted here because I can't understand why the if statement in my loop only ever runs once. The FFT code could be doing something else here, I just want to understand the interaction between the interrupt and the canProcess flag and the loop.

Hope that makes sense.

No, it just for interrupts. Starting/stopping timer has nothing to do with interrupts.

@b707 Sure, makes sense

OK I've fixed it, or at least I've gone around it. I still don't understand the problem however.

I commented parts of the FFT function bit by bit and realised that the problem was somewhere in the scaling section - when the input data is normalised to a range of [-512, 512]. As my data is sourced from digitalRead I know it is all 1s and 0s. So I removed all of the scaling code and simply set all my data points to -512 or 512 at the point of reading them.

This worked perfectly and I am now identifying input frequencies 100 times per second :tada: The last obstacle is to check whether I can identify multiple frequencies at the same time. I'll have to wait for a delivery of more Arduinos with my fingers crossed.

Thank you @b707 @alto777.

The offending code:

if (data_mag < 1024) {
  while (temp11 < 1024) {
    temp11 = temp11 << 1;
    scale = scale + 1;
  }
}

The fully updated code:

int IRPin = 9;

volatile int dataIndex = 0;
volatile int dataOne[64] = {};
volatile int dataTwo[64] = {};
volatile bool fillingDataOne = true;
volatile bool canProcess = false;

void setup() {
  pinMode(IRPin, INPUT);

  cli();
  TCCR1A = 0;
  TCCR1B = 0;
  TCCR1B |= _BV(CS10); // No prescaler
  TCCR1B |= _BV(WGM12); // Count up to OCR1A
  OCR1A = 2499; // 6400Hz (64 samples every 10 milliseconds)
  TCNT1 = 0; // Reset counter
  TIMSK1 |= _BV(OCIE1A); // Interrupt
  sei();

  Serial.begin(115200);
}

ISR(TIMER1_COMPA_vect) {
  int level = (digitalRead(IRPin) == 1) ? 512 : -512;
  if (fillingDataOne) {
    dataOne[dataIndex++] = level;
  } else {
    dataTwo[dataIndex++] = level;
  }
  if (dataIndex == 64) {
    dataIndex = 0;
    canProcess = true;
    fillingDataOne = !fillingDataOne;
  }
}

void loop() {
  if (canProcess) {
    canProcess = false;
    float f = Approx_FFT(fillingDataOne ? dataTwo : dataOne, 64, 6400);
    Serial.println(f);
  }
}

// Everything below here is from Approx_FTT

byte isin_data[128] = {
  0, 1, 3, 4, 5, 6, 8, 9, 10, 11,
  13, 14, 15, 17, 18, 19, 20, 22,
  23, 24, 26, 27, 28, 29, 31, 32,
  33, 35, 36, 37, 39, 40, 41, 42,
  44, 45, 46, 48, 49, 50, 52, 53,
  54, 56, 57, 59, 60, 61, 63, 64,
  65, 67, 68, 70, 71, 72, 74, 75,
  77, 78, 80, 81, 82, 84, 85, 87,
  88, 90, 91, 93, 94, 96, 97, 99,
  100, 102, 104, 105, 107, 108,
  110, 112, 113, 115, 117, 118,
  120, 122, 124, 125, 127, 129,
  131, 133, 134, 136, 138, 140,
  142, 144, 146, 148, 150, 152,
  155, 157, 159, 161, 164, 166,
  169, 171, 174, 176, 179, 182,
  185, 188, 191, 195, 198, 202,
  206, 210, 215, 221, 227, 236
};

unsigned int Pow2[14] = {
  1, 2, 4, 8, 16, 32, 64, 128,
  256, 512, 1024, 2048, 4096
};

byte RSSdata[20] = {
  7, 6, 6, 5, 5, 5, 4, 4, 4, 4,
  3, 3, 3, 3, 3, 3, 3, 2, 2, 2
};

/*
Code   to perform High speed and Accurate FFT on arduino,
setup:

1. in[]     :   Data array,
2. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
3.   Frequency: sampling frequency required as input (Hz)

It will by default return   frequency with max aplitude,
if you need complex output or magnitudes uncomment   required sections

If sample size is not in power of 2 it will be clipped   to lower side of number. 
i.e, for 150 number of samples, code will consider   first 128 sample, remaining sample  will be omitted.
For Arduino nano, FFT of   more than 256 sample not possible due to mamory limitation 
Code by ABHILASH
Contact:   abhilashpatel121@gmail.com
Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/

Update(06/05/21):   Correction made for support on Arduino Due
*/

float Approx_FFT(int in[], int N, float Frequency) {
  int a, c1, f, o, x = 0;
  byte scale, check = 0;

  for (int i = 0; i < 12; i++)  //calculating the levels
  {
    if (Pow2[i] <= N) {
      o = i;
    }
  }
  a = Pow2[o];
  int out_r[a];   //real part of transform
  int out_im[a];  //imaginory part of transform

  scale = 1024;

  x = 0;
  for (int b = 0; b < o; b++)  // bit reversal order stored   in im_out array
  {
    c1 = Pow2[b];
    f = Pow2[o] / (c1 + c1);
    for (int j = 0; j < c1; j++) {
      x = x + 1;
      out_im[x] = out_im[j] + f;
    }
  }

  for (int i = 0; i < a; i++)  // update input array as per bit reverse order
  {
    out_r[i] = in[out_im[i]];
    out_im[i] = 0;
  }

  int i10, i11, n1, tr, ti;
  float e;
  int c, s, temp4;
  for (int i = 0; i < o; i++)  //fft
  {
    i10 = Pow2[i];                // overall values of sine/cosine
    i11 = Pow2[o] / Pow2[i + 1];  // loop with similar sine cosine
    e = 1024 / Pow2[i + 1];       //1024 is equivalent   to 360 deg
    e = 0 - e;
    n1 = 0;

    for (int j = 0; j < i10; j++) {
      c = e * j;  //c is angle as where 1024 unit is 360 deg
      while (c < 0) {
        c = c + 1024;
      }
      while (c > 1024) {
        c = c - 1024;
      }

      n1 = j;

      for (int k = 0; k < i11; k++) {
        temp4 = i10 + n1;
        if (c == 0) {
          tr = out_r[temp4];
          ti = out_im[temp4];
        } else if (c == 256) {
          tr = -out_im[temp4];
          ti = out_r[temp4];
        } else if (c == 512) {
          tr = -out_r[temp4];
          ti = -out_im[temp4];
        } else if (c == 768) {
          tr = out_im[temp4];
          ti = -out_r[temp4];
        } else if (c == 1024) {
          tr = out_r[temp4];
          ti = out_im[temp4];
        } else {
          tr = fast_cosine(out_r[temp4], c) - fast_sine(out_im[temp4], c);  //the   fast sine/cosine function gives direct (approx) output for A*sinx
          ti = fast_sine(out_r[temp4], c) + fast_cosine(out_im[temp4], c);
        }

        out_r[n1 + i10] = out_r[n1] - tr;
        out_r[n1] = out_r[n1] + tr;
        if (out_r[n1] > 15000 || out_r[n1] < -15000) {
          check = 1;
        }  //check for int size, it can handle only +31000   to -31000,

        out_im[n1 + i10] = out_im[n1] - ti;
        out_im[n1] = out_im[n1] + ti;
        if (out_im[n1] > 15000 || out_im[n1] < -15000) {
          check = 1;
        }

        n1 = n1 + i10 + i10;
      }
    }

    if (check == 1) {  //   scalling the matrics if value higher than 15000 to prevent varible from overflowing
      for (int i = 0; i < a; i++) {
        out_r[i] = out_r[i] >> 1;
        out_im[i] = out_im[i] >> 1;
      }
      check = 0;
      scale = scale - 1;  //   tracking overall scalling of input data
    }
  }

  if (scale > 128) {
    scale = scale - 128;
    for (int i = 0; i < a; i++) {
      out_r[i] = out_r[i] >> scale;
      out_im[i] = out_im[i] >> scale;
    }
    scale = 0;
  }  // revers   all scalling we done till here,
  else {
    scale = 128 - scale;
  }  //   in case of nnumber getting higher than 32000, we will represent in as multiple of   2^scale

  /*
  for(int i=0;i<a;i++)
  {
  Serial.print(out_r[i]);Serial.print("\  ");                     // un comment to print RAW o/p    
  Serial.print(out_im[i]);   
  Serial.print("i");Serial.print("\ "); 
  Serial.print("*2^");Serial.println(scale);   
  }
  */

  //---> here onward out_r contains amplitude and our_in conntains   frequency (Hz)
  int fout, fm, fstp;
  float fstep;
  fstep = Frequency / N;
  fstp = fstep;
  fout = 0;
  fm = 0;

  for (int i = 1; i < Pow2[o - 1]; i++)  // getting amplitude from compex   number
  {
    out_r[i] = fastRSS(out_r[i], out_im[i]);
    // Approx RSS function used to calculated magnitude quickly

    out_im[i] = out_im[i - 1] + fstp;
    if (fout < out_r[i]) {
      fm = i;
      fout = out_r[i];
    }

    // un comment to   print Amplitudes (1st value (offset) is not printed)
    //Serial.print(out_r[i]);   Serial.print("\  ");
    //Serial.print("*2^");Serial.println(scale);
  }

  float fa, fb, fc;
  fa = out_r[fm - 1];
  fb = out_r[fm];
  fc = out_r[fm + 1];
  fstep = (fa * (fm - 1) + fb * fm + fc * (fm + 1)) / (fa + fb + fc);

  return (fstep * Frequency / N);
}

//---------------------------------fast   sine/cosine---------------------------------------//

int fast_sine(int Amp, int th) {
  int temp3, m1, m2;
  byte temp1, temp2, test, quad, accuracy;
  accuracy = 2;  // set it value from 1 to 7, where 7 being most accurate but slowest
  //   accuracy value of 5 recommended for typical applicaiton
  while (th > 1024) {
    th = th - 1024;
  }  // here 1024 = 2*pi or 360 deg
  while (th < 0) {
    th = th + 1024;
  }
  quad = th >> 8;

  if (quad == 1) {
    th = 512 - th;
  } else if (quad == 2) {
    th = th - 512;
  } else if (quad == 3) {
    th = 1024 - th;
  }

  temp1 = 0;
  temp2 = 128;  //2 multiple
  m1 = 0;
  m2 = Amp;

  temp3 = (m1 + m2) >> 1;
  Amp = temp3;
  for (int i = 0; i < accuracy; i++) {
    test = (temp1 + temp2) >> 1;
    temp3 = temp3 >> 1;
    if (th > isin_data[test]) {
      temp1 = test;
      Amp = Amp + temp3;
      m1 = Amp;
    } else if (th < isin_data[test]) {
      temp2 = test;
      Amp = Amp - temp3;
      m2 = Amp;
    }
  }

  if (quad == 2) {
    Amp = 0 - Amp;
  } else if (quad == 3) {
    Amp = 0 - Amp;
  }
  return (Amp);
}

int fast_cosine(int Amp, int th) {
  th = 256 - th;  //cos th = sin (90-th) formula
  return (fast_sine(Amp, th));
}

//--------------------------------------------------------------------------------//

//--------------------------------Fast   RSS----------------------------------------//
int fastRSS(int a, int b) {
  if (a == 0 && b == 0) {
    return (0);
  }
  int min, max, temp1, temp2;
  byte clevel;
  if (a < 0) {
    a = -a;
  }
  if (b < 0) {
    b = -b;
  }
  clevel = 0;
  if (a > b) {
    max = a;
    min = b;
  } else {
    max = b;
    min = a;
  }

  if (max > (min + min + min)) {
    return max;
  } else {
    temp1 = min >> 3;
    if (temp1 == 0) {
      temp1 = 1;
    }
    temp2 = min;
    while (temp2 < max) {
      temp2 = temp2 + temp1;
      clevel = clevel + 1;
    }
    temp2 = RSSdata[clevel];
    temp1 = temp1 >> 1;
    for (int i = 0; i < temp2; i++) {
      max = max + temp1;
    }
    return (max);
  }
}

This topic was automatically closed 180 days after the last reply. New replies are no longer allowed.