typedef enum fft_dir {
  FFT_FORWARD,    /* kernel uses "-1" sign */
  FFT_INVERSE     /* kernel uses "+1" sign */
} fft_dir;

void window (float complex *data, unsigned int log2_N, byte windowType, fft_dir direction) {
  unsigned int N = 1 << log2_N;
  unsigned int Nd2 = N >> 1;
  for (int i = 0; i < Nd2; i++) {
    double indexMinusOne = double(i);
    double ratio = (indexMinusOne / (N - 1));
    double weighingFactor = 1.0;
    // Compute and record weighting factor
    switch (windowType) {
      case RECTANGLE: // rectangle (box car)
        weighingFactor = 1.0;
        break;
      case HAMMING: // hamming
        weighingFactor = 0.54 - (0.46 * cos(2 * PI * ratio));
        break;
      case HANN: // hann
        weighingFactor = 0.54 * (1.0 - cos(2 * PI * ratio));
        break;
      case TRIANGLE: // triangle (Bartlett)
        weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - ((N - 1) / 2.0))) / (N - 1));
        break;
      case NUTTALL: // nuttall
        weighingFactor = 0.355768 - (0.487396 * (cos(2 * PI * ratio))) + (0.144232 * (cos(4 * PI * ratio))) - (0.012604 * (cos(6 * PI * ratio)));
        break;
      case BLACKMAN: // blackman
        weighingFactor = 0.42323 - (0.49755 * (cos(2 * PI * ratio))) + (0.07922 * (cos(4 * PI * ratio)));
        break;
      case BLACKMAN_NUTTALL: // blackman nuttall
        weighingFactor = 0.3635819 - (0.4891775 * (cos(2 * PI * ratio))) + (0.1365995 * (cos(4 * PI * ratio))) - (0.0106411 * (cos(6 * PI * ratio)));
        break;
      case BLACKMAN_HARRIS: // blackman harris
        weighingFactor = 0.35875 - (0.48829 * (cos(2 * PI * ratio))) + (0.14128 * (cos(4 * PI * ratio))) - (0.01168 * (cos(6 * PI * ratio)));
        break;
      case FLT_TOP: // flat top
        weighingFactor = 0.2810639 - (0.5208972 * cos(2 * PI * ratio)) + (0.1980399 * cos(4 * PI * ratio));
        break;
      case WELCH: // welch
        weighingFactor = 1.0 - sq((indexMinusOne - (N - 1) / 2.0) / ((N - 1) / 2.0));
        break;
    }
    if (direction == FFT_FORWARD) {
      data[i] *= weighingFactor;
      data[N - (i + 1)] *= weighingFactor;
    }
    else {
      data[i] /= weighingFactor;
      data[N - (i + 1)] /= weighingFactor;
    }
  }
}

void ffti_shuffle_f(float complex *data, unsigned int log2_N)
{
  /*
     Basic Bit-Reversal Scheme:

     The incrementing pattern operations used here correspond
     to the logic operations of a synchronous counter.

     Incrementing a binary number simply flips a sequence of
     least-significant bits, for example from 0111 to 1000.
     So in order to compute the next bit-reversed index, we
     have to flip a sequence of most-significant bits.
  */

  unsigned int N = 1 << log2_N;   /* N */
  unsigned int Nd2 = N >> 1;      /* N/2 = number range midpoint */
  unsigned int Nm1 = N - 1;       /* N-1 = digit mask */
  unsigned int i;                 /* index for array elements */
  unsigned int j;                 /* index for next element swap location */

  for (i = 0, j = 0; i < N; i++) {
    if (j > i) {
      float complex tmp = data[i];
      data[i] = data[j];
      data[j] = tmp;
    }

    /*
       Find least significant zero bit
    */

    unsigned int lszb = ~i & (i + 1);

    /*
       Use division to bit-reverse the single bit so that we now have
       the most significant zero bit

       N = 2^r = 2^(m+1)
       Nd2 = N/2 = 2^m
       if lszb = 2^k, where k is within the range of 0...m, then
           mszb = Nd2 / lszb
                = 2^m / 2^k
                = 2^(m-k)
                = bit-reversed value of lszb
    */

    unsigned int mszb = Nd2 / lszb;

    /*
       Toggle bits with bit-reverse mask
    */

    unsigned int bits = Nm1 & ~(mszb - 1);
    j ^= bits;
  }
}



void ffti_evaluate_f(float complex *data, unsigned int log2_N, fft_dir direction)
{
  /*
     In-place FFT butterfly algorithm

     input:
         A[] = array of N shuffled complex values where N is a power of 2
     output:
         A[] = the DFT of input A[]

     for r = 1 to log2(N)
         m = 2^r
         Wm = exp(−j2π/m)
         for n = 0 to N-1 by m
             Wmk = 1
             for k = 0 to m/2 - 1
                 u = A[n + k]
                 t = Wmk * A[n + k + m/2]
                 A[n + k]       = u + t
                 A[n + k + m/2] = u - t
                 Wmk = Wmk * Wm

     For inverse FFT, use Wm = exp(+j2π/m)
  */

  unsigned int N = 1 << log2_N;
  double theta_2pi = (direction == FFT_FORWARD) ? -PI : PI;
  theta_2pi *= 2;

  for (unsigned int r = 1; r <= log2_N; r++)
  {
    unsigned int m = 1 << r;
    unsigned int md2 = m >> 1;
    double theta = theta_2pi / m;
    double re = cos(theta);
    double im = sin(theta);
    double complex Wm = re + I * im;
    for (unsigned int n = 0; n < N; n += m) {
      double complex Wmk = 1.0f;   /* Use double for precision */
      for (unsigned int k = 0; k < md2; k++) {
        unsigned int i_e = n + k;
        unsigned int i_o = i_e + md2;
        double complex t = Wmk * data[i_o];
        data[i_o] = data[i_e] - t;
        data[i_e] = data[i_e] + t;
        t = Wmk * Wm;
        Wmk = t;
      }
    }
  }
}

void ffti_f(float complex data[], unsigned int log2_N, fft_dir direction)
{
  ffti_shuffle_f(data, log2_N);
  ffti_evaluate_f(data, log2_N, direction);
}

void acquireSound () {
  for (int i = 0; i < SAMPLES; i++) {
    unsigned long chrono = micros();
    data[i] = analogRead(MIC);
    while (micros() - chrono < sampling_period_us); // do nothing
  } // 12.8 ms
}

void displaySpectrum () {
  window (data, LOG2SAMPLE, HAMMING, FFT_FORWARD);
  ffti_f(data, LOG2SAMPLE, FFT_FORWARD);
  // window + FFT last 10ms
  int nFreq = SAMPLES / 2;
  int hauteur = display.height();
  display.fillScreen(TFT_BLACK);
  display.setTextColor(TFT_WHITE);
  display.setTextSize(1);
  const float freqs[7] = {.5, 1, 2, 3, 4, 6, 8};
  for (int i = 0; i < 7; i++) {
    int pos =  FREQ2IND * freqs[i] * 1.9 ;
    if (i == 0) display.drawString(".5", pos - 5, hauteur - 13, 2);
    else display.drawString(String(freqs[i], 0), pos, hauteur - 13, 2);
  }
  display.drawString("kHz", display.width() - 22, hauteur - 13, 2);

  const int decal = 15;
  // Affichage du spectre
  int ampmax = 0;
  int imax = 0;
  for (int i = 2; i < nFreq; i++) {
    int amplitude = (int)creal(data[i]) / COEF;
    if (amplitude > ampmax) {
      ampmax = amplitude;
      imax = i;
    }
    amplitude = min(amplitude, hauteur - decal);
    display.drawFastVLine (i * 2, hauteur - amplitude - decal, amplitude, TFT_GREEN);
  }

  // Affichage fréquence de plus forte amplitude
  if (ampmax > 20) {
    display.setTextSize(2);
    display.setTextColor(TFT_BLUE);
    int freqmax = 985.0 * imax * MAX_FREQ / SAMPLES;
    char texte[14];
    sprintf(texte, "%4d Hz: %3d", freqmax, ampmax);
    display.drawString(texte, display.width() - 145, 0, 1);
  }

  // Affichage et décroissance des valeurs peak
  for (int i = 2; i < nFreq - 4; i = i + 4) {
    int amplitude = (int)creal(data[i]);
    amplitude = max(amplitude, (int)creal(data[i + 1]));
    amplitude = max(amplitude, (int)creal(data[i + 2]));
    amplitude = max(amplitude, (int)creal(data[i + 3]));
    amplitude /= COEF;
    if (amplitude > peak[i]) peak[i] = amplitude;
    if (peak[i] > hauteur - decal) peak[i] = hauteur - decal;
    if (peak[i] > 8) {
      peak[i] -= 2; // Decay the peak
      display.drawFastHLine(i * 2, hauteur - peak[i] - decal, 4, TFT_RED);
    }
  }
}

void rectangleGR (int x0, int w, int h) {
  const uint16_t colors[8] = {TFT_DARKGREEN, TFT_GREEN, TFT_GREENYELLOW, TFT_YELLOW,
                              TFT_GOLD, TFT_ORANGE, TFT_RED, TFT_BROWN
                             };
  const int decal = 15;
  int x = x0 * (w + 2);
  int dh = (display.height() - decal) / 8;
  for (int i = 1; i < 8; i++) {
    if (h > i * dh) {
      int y = display.height() - decal - i * dh;
      display.fillRect(x, y, w, dh, colors[i - 1]);
    } else {
      int y = display.height() - decal - h;
      display.fillRect(x, y, w, h - (i - 1) * dh, colors[i - 1]);
      break;
    }
  }
}

void displaySpectrum2 () {
  window (data, LOG2SAMPLE, HAMMING, FFT_FORWARD);
  ffti_f(data, LOG2SAMPLE, FFT_FORWARD);
  int nFreq = SAMPLES / 2;
  int hauteur = display.height();
  display.fillScreen(TFT_BLACK);
  display.setTextColor(TFT_WHITE);
  display.setTextSize(1);
  const float freqs[6] = {.5, 1, 1.5, 2, 3, 4};
  for (int i = 0; i < 6; i++) {
    int pos =  FREQ2IND * freqs[i] * 3.5 ;
    if (i == 0) display.drawString(".5", pos - 5, hauteur - 13, 2);
    else if (i == 1) display.drawString("1", pos - 10, hauteur - 13, 2);
    else if (i == 2) display.drawString("1.5", pos - 5, hauteur - 13, 2);
    else display.drawString(String(freqs[i], 0), pos, hauteur - 13, 2);
  }
  display.drawString("kHz", display.width() - 22, hauteur - 13, 2);

  const int decal = 15;
  // Affichage du spectre
  for (int i = 0; i < 16; i++) {
    int amplitude = 0;
    for (int j = 0; j < 4; j++) amplitude += abs((int)creal(data[i * 4 + j + 2]) / COEF);
    amplitude = min(amplitude, hauteur - decal);
    rectangleGR (i, 13, amplitude);
  }
}

void displayAmplitudeBars() {
  const uint16_t colors[8] = {TFT_DARKGREEN, TFT_GREEN, TFT_GREENYELLOW, TFT_YELLOW,
                              TFT_GOLD, TFT_ORANGE, TFT_RED, TFT_BROWN
                             };
  int ampmax = -5000;
  int ampmin =  5000;
  int N = display.width();
  int decal = 15;
  int hauteur = display.height() - decal;
  for (int i = 0; i < N; i++) {
    sound[i] = analogRead(MIC);
    delayMicroseconds(50);
  }
  display.fillScreen(TFT_BLACK);
  int prevamp = 1;
  for (int i = 0; i < N; i++) {
    uint16_t colour;
    //    colour = TFT_GREEN;
    //    colour = colors[(i / 10) % 8];
    uint16_t red  = 0b11111;
    uint16_t red2 = 0b00001;
    if (i < N / 2) colour = (red - i * (red - red2) * 2 / N) << 11;
    else colour = (i + N / 2) * 32 / N;
    int amplitude = (sound[i] - SILENCE) / COEF2;
    constrain(amplitude, -hauteur / 2, hauteur / 2);
    if (amplitude > ampmax) ampmax = amplitude;
    if (amplitude < ampmin) ampmin = amplitude;
    display.drawFastVLine (i, hauteur / 2 - amplitude + decal, abs(amplitude), colour);
    if (amplitude * prevamp < 0 && i > 0)
      display.drawLine(i - 1, hauteur / 2 - prevamp + decal, i, hauteur / 2 - amplitude + decal, colour);
  }
  // Affichage plus forte amplitude
  display.setTextColor(TFT_BLUE);
  display.setTextSize(2);
  char texte[25];
  sprintf(texte, "Min: %3d, Max: %3d", ampmin, ampmax);
  display.drawString(texte, 10, 0, 1);
  delay(50);
}

void displayEnvelope(bool erase) {
  static int i = 0;
  if (i == 1 || erase) display.fillScreen(TFT_BLACK);
  const int decal = 15;
  const int hauteur = display.height();
  const int N = display.width();
  static unsigned int maxP2P = 0;
  static unsigned int minP2P = 5000;
  static unsigned int maxP2Pold = 0;
  static unsigned int minP2Pold = 5000;
  static int amplitudeold = 0;

  if (i >= N || erase) {
    i = 0;
    maxP2P = 0;
    minP2P = 5000;
    maxP2Pold = 0;
    minP2Pold = 5000;
  }

  unsigned int signalMax = 0;
  unsigned int signalMin = 4096;
  unsigned long chrono = micros(); // Sample window 10ms
  while (micros() - chrono < 10000ul) {
    int sample = analogRead(MIC) / COEF3;
    if (sample > signalMax) signalMax = sample;
    else if (sample < signalMin) signalMin = sample;
  }
  unsigned int peakToPeak = signalMax - signalMin;
  int amplitude = peakToPeak - SENSITIVITY;
  if (amplitude < 0) amplitude = 0;
  else if (amplitude > 500) amplitude = 500;
  int ampPrev = map(amplitudeold, 0, 500, hauteur - 1, decal);
  int amp     = map(amplitude,    0, 500, hauteur - 1, decal);
  if (amplitude > maxP2P) maxP2P = amplitude;
  if (amplitude < minP2P) minP2P = amplitude;
  if (i > 0) display.drawLine (i - 1, ampPrev, i, amp, TFT_GREEN);
  amplitudeold = amplitude;
  // Affichage plus forte amplitude
  if ((maxP2Pold != maxP2P) || (minP2Pold != minP2P) || (i == 0)) {
    display.fillRect(0, 0, N, 15, TFT_BLACK);
    display.setTextColor(TFT_BLUE);
    char texte[25];
    sprintf(texte, "Min : %3d, Max: %3d", minP2P, maxP2P);
    display.drawString(texte, 0, 0, 1);
  }
  maxP2Pold = maxP2P;
  minP2Pold = minP2P;
  ++i;
}

void displayRunningEnvelope(bool erase) {
  if (erase) display.fillScreen(TFT_BLACK);
  static int ampMaxOld = 0;
  static int ampMinOld = 1024;
  const int decal = 15;
  const int hauteur = display.height();
  const int N = display.width();
  static int sound[240] = {0};
  int save = sound[0];
  for (int i = 0; i < N - 1 ; i++) sound[i] = sound[i + 1];

  unsigned int signalMax = 0;
  unsigned int signalMin = 4096;
  unsigned long chrono = micros(); // Sample window 10ms
  while (micros() - chrono < 10000ul) {
    int sample = analogRead(MIC) / COEF3;
    if (sample > signalMax) signalMax = sample;
    else if (sample < signalMin) signalMin = sample;
  }
  unsigned int peakToPeak = signalMax - signalMin;
  int x = peakToPeak - SENSITIVITY;
  if (x < 0)   x = 0;
  if (x > 500) x = 500;
  sound[N - 1] = map(x, 0, 500, 1, hauteur - decal);

  int ampMin = 500;
  int ampMax = 0;
  display.drawLine (0, hauteur - save, 1, hauteur - sound[0], TFT_BLACK);
  for (int j = 1; j < N; j++) {
    if (sound[j] > ampMax) ampMax = sound[j];
    if (sound[j] < ampMin) ampMin = sound[j];
    if (ampMin == 1) ampMin = 0;
    display.drawLine (j, hauteur - sound[j - 1], j + 1, hauteur - sound[j], TFT_BLACK);
    display.drawLine (j - 1, hauteur - sound[j - 1], j, hauteur - sound[j], TFT_GREEN);
  }
  // Affichage bornes d'amplitude
  if ((ampMaxOld != ampMax) || (ampMinOld != ampMin) || erase) {
    display.fillRect(0, 0, N, 15, TFT_BLACK);
    char texte[25];
    sprintf(texte, "Min : %3d, Max: %3d", ampMin, ampMax);
    display.drawString(texte, 0, 0, 1);
    ampMaxOld = ampMax;
    ampMinOld = ampMin;
  }
}

/*
    Get "ambient" volume
*/
void calibrate() {
  for (int i = 0; i < 200; i++) backgroundSound += analogRead(MIC);
  backgroundSound /= 200;
  //  Serial.print("Threshold set at ");
  //  Serial.println(backgroundSound);
}

bool soundDetected() {
  return abs(analogRead(MIC) - backgroundSound) >= SOUND_THRESHOLD;
}

void displaySpectrogram (bool erase, bool displayMax) {
#define nCoul 20
  if (erase) {
    display.fillScreen(TFT_BLACK);
    display.setTextColor(TFT_BLUE);
    display.setTextSize(1);
    display.drawString("WAITING", 65, 40, 4);
    display.drawString("FOR SOUND", 40, 70, 4);
  }
  const int rectW = display.width() / nTimes;
  const int nbFreqD = 60;
  const int decal = 15;
  const int rectH = (display.height() - decal) / nbFreqD;
  const uint16_t colors[nCoul] = {TFT_BLACK, TFT_BLACK, TFT_BLACK, TFT_BLACK,
                                  TFT_BLACK, TFT_BLACK, TFT_BLACK, TFT_BLACK,
                                  TFT_BLACK, TFT_BLACK, TFT_DARKGREEN, TFT_VIOLET,
                                  TFT_BLUE, TFT_SKYBLUE, TFT_GREEN, TFT_GREENYELLOW,
                                  TFT_YELLOW, TFT_GOLD, TFT_ORANGE, TFT_RED //, TFT_BROWN
                                 };
  // Wait until a sound arrives
  if (erase) delay(300); // to avoid the sound of the button
  if (!soundDetected()) return;
  // Acquire samples over the 700ms range
  for (int j = 0; j < nTimes; j++) {
    for (int i = 0; i < SAMPLES; i++) {
      unsigned long chrono = micros();
      data2D[i][j] = analogRead(MIC);
      while (micros() - chrono < sampling_period_us); // do nothing
    } // 12.8 ms
  }
  // Fourier transform and display each sound
  display.fillScreen(TFT_BLACK);
  for (int k = 1; k < 5; k++) display.drawString(String(k), 0, display.height() - k * 25 - decal - 8, 2);
  for (int k = 0; k < 8; k++) // '8' is for 700ms
    display.drawString(String((float)(k / 10.0), 1), k * 30, display.height() - 13, 2);
  display.drawString("s", display.width() - 7, display.height() - 13, 2);

  //Histogram
  // int nb[22] = {0};

  int indMax0[nTimes] = {0}; int thresh0 =  400;
  int indMax1[nTimes] = {0}; int thresh1 =  800;
  int indMax2[nTimes] = {0}; int thresh2 = 2000;
  int indMax3[nTimes] = {0}; int thresh3 = 4000;

  for (int j = 0; j < nTimes; j++) {
    for (int i = 0; i < SAMPLES; i++) data[i] = (float)data2D[i][j] + I * 0.0f;
    window (data, LOG2SAMPLE, HAMMING, FFT_FORWARD);
    ffti_f(data, LOG2SAMPLE, FFT_FORWARD);

    int amplitude[nbFreqD] = {0};
    for (int i = 3; i < nbFreqD; i++) {
      amplitude[i] = abs((int)creal(data[i]));
      if (!displayMax) {
        // draw spectrogram
        int iCoul = map(amplitude[i], 0, 160, 0, nCoul - 1);
        if (iCoul < 0) iCoul = 0;
        if (iCoul > nCoul - 1) iCoul = nCoul - 1;
        //    ++nb[iCoul];
        int y = display.height() - decal - i * rectH;
        display.fillRect((j + 2) * rectW + 3, y, rectW, rectH, colors[iCoul]);
      }
    }

    // Search the envelopes in time - frequency
    if (displayMax) {
      for (int i = nbFreqD; i > 2; i--) {
        if (amplitude[i] > thresh0 && indMax0[j] == 0) indMax0[j] = i;
        if (amplitude[i] > thresh1 && indMax1[j] == 0) indMax1[j] = i;
        if (amplitude[i] > thresh2 && indMax2[j] == 0) indMax2[j] = i;
        if (amplitude[i] > thresh3 && indMax3[j] == 0) indMax3[j] = i;
      }
    }

  }
  // Print histogram
  //  for (int h = 0; h < nCoul; h++) Serial.printf("Colour %2d : %d\n", h, nb[h]);
  //  Serial.println();
  display.drawString("kHz", 15, 0, 2);

  if (!displayMax) {
    // Display scale
    int dx = 13;
    int dy = (display.height() - 2 * decal) / 11;
    for (int i = 9; i < nCoul; i++) {
      int y = display.height() - (i - 9) * dy - decal * 2 - 3;
      display.fillRect(display.width() - dx, y, dx, dy - 3, colors[i]);
      display.drawRect(display.width() - dx, y, dx, dy - 3, TFT_WHITE);
    }    
  } else {
    // Display envelopes
    int xp = 2 * rectW + 3;
    int yp = display.height() - decal - indMax0[0] * rectH;
    for (int j = 1; j < nTimes; j++) {
      int x = (j + 2) * rectW + 3;
      int y = display.height() - decal - indMax0[j] * rectH;
      display.drawLine(xp, yp, x, y, TFT_RED);
      xp = x;
      yp = y;
    }
    xp = 2 * rectW + 3;
    yp = display.height() - decal - indMax1[0] * rectH;
    for (int j = 1; j < nTimes; j++) {
      int x = (j + 2) * rectW + 3;
      int y = display.height() - decal - indMax1[j] * rectH;
      display.drawLine(xp, yp, x, y, TFT_CYAN);
      xp = x;
      yp = y;
    }
    xp = 2 * rectW + 3;
    yp = display.height() - decal - indMax2[0] * rectH;
    for (int j = 1; j < nTimes; j++) {
      int x = (j + 2) * rectW + 3;
      int y = display.height() - decal - indMax2[j] * rectH;
      display.drawLine(xp, yp, x, y, TFT_BLUE);
      xp = x;
      yp = y;
    }
    xp = 2 * rectW + 3;
    yp = display.height() - decal - indMax3[0] * rectH;
    for (int j = 1; j < nTimes; j++) {
      int x = (j + 2) * rectW + 3;
      int y = display.height() - decal - indMax3[j] * rectH;
      display.drawLine(xp, yp, x, y, TFT_GREEN);
      xp = x;
      yp = y;
    }
  }

}