1024 point FFT, 30% more efficient (experimental)

Status
Not open for further replies.
Frequency returned different

With the old library the following code returned a frequency of 944 Hz with an input of 941 Hz tuning fork.
With the new library the frequency was 936 Hz.

Code:
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>

int maxbin, k;

const int myInput = AUDIO_INPUT_MIC;

// Create the Audio components.  These should be created in the
// order data flows, inputs/sources -> processing -> outputs
//
AudioInputI2S       audioInput;         // audio shield: mic or line-in
AudioAnalyzeFFT1024  myFFT;
AudioOutputI2S      audioOutput;        // audio shield: headphones & line-out

// Create Audio connections between the components

AudioConnection c1(audioInput, 0, audioOutput, 0);
AudioConnection c2(audioInput, 0, myFFT, 0);
AudioConnection c3(audioInput, 1, audioOutput, 1);

// Create an object to control the audio shield.
 
AudioControlSGTL5000 audioShield;


void setup() {
  // Audio connections require memory to work.  For more
  // detailed information, see the MemoryAndCpuUsage example
  AudioMemory(12);
  myFFT.windowFunction(AudioWindowBlackmanHarris1024);
  
  // Enable the audio shield and set the output volume.
  audioShield.enable();
  audioShield.inputSelect(myInput);
  audioShield.volume(0.6);
}

void loop() {
  if (myFFT.available()) {
    // each time new FFT data is available
    // print it all to the Arduino Serial Monitor
     k=0;
     maxbin=0;
    
    
    for (int i=0; i<512; i++) {
      if (myFFT.output[i]>maxbin) {
        maxbin=myFFT.output[i];
        k=i;
      } 
     
    }
    Serial.println();
    Serial.print("Max Bin ");
    Serial.print(k);
    Serial.println();
    double y1=double(myFFT.output[k-1]);
    double y2=double(myFFT.output[k]);
    double y3=double(myFFT.output[k+1]);
   double d=(y3-y1)/((2.0*(2.0*y2-y1-y3)));
    Serial.print("Diff ");
    Serial.print(d);
    Serial.println();
    double Freq=(AUDIO_SAMPLE_RATE_EXACT/1024.0)*(double(k)+d);
    Serial.print("Freq ");
    Serial.print(Freq);
    Serial.println();
  }
}
 
I'll put up results with that in shortly. I have made edits to sum the magnitude before each line prints, also considering showing Peak. And I have stopped printing trailing '-' so you know when the bins zero.

FFT1024=51,52 all=52.98,53.92 Memory: 4,8 Ready:2256960 Not:401981663 # Samples:120 Sample Time:1392 Time per sample:11

274]FFT: - 13 23 24 15 8 6 5 4 3 3 3 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 3 4 5 12 21 24 18 8 4 3 2 2 1 1 1 1 1
98]FFT:1 13 45 35
97]FFT: - 13 45 35
141]FFT:3 11 21 34 33 21 12 4 1
100]FFT: - - - 5 37 45 12
100]FFT: - - - 5 37 45 12
 
Last edited:
@cartere. I tried to perform a test on both the original and my fft. I created 200 sine waves of exactly 941 Hz, every time with a different phase offset
With a run of 200, I got these results:
old: 942.5619 +/- 0.0050
mine: 942.5618 +/- 0.0057
Anyway, you are looking to a resolution within 10% of the bin width. If you are looking that accurate, better use interpolation algorithms that also incorporate the phase.
Even without phase, for this kind of resolution you might want to include the window effects/scalloping loss in your frequency estimation.
Edit: Just did some calculations, using your equation, based on a Blackman-Harris window, you expect to measure 941.82 Hz

@defragster: Sorry, what am I supposed to see in the graph, is it good, or what is the problem?
If the problem is the amplitudes, then it might be due to the fact that the frequency bins of an fft are uncorrelated by definition. Therefore you should use sqrt(sum(x^2)). If you do that then the amplitude remains 100.
If it is the difference in the peak, you might want to try the AudioWindowFlattop1024, this should have no scalloping loss.
 
Last edited:
@kpc - just getting a handle on this - enough for practical use. This sample snippet it the PJRC algorithm - I have altered the 'stats line that will shows when KPC is used'

Good Feedback - Thanks! I was showing I was summing the bins for a quick glance at cumulative value from the sample I started with. Before they seemed to add to the specified amplitude, I will use what you added "sqrt(sum(x^2))". I was looking for ways to compare the PJRC and KPC versions [edit: as well as the various algorithms] at a glance for a utility/sanity test since I had it at hand and want to learn enough that I can use it going forward. If you have other feedback [or a pointer to where I could learn] on what I can do to differentiate algorithms, please let me know. I want to keep the output simple and informative - but if I can add workload that simulates normal processing needed to use the data - that will also be good for judging the available CPU power after the FFT is done.
 
Last edited:
@cartere. I tried to perform a test on both the original and my fft. I created 200 sine waves of exactly 941 Hz, every time with a different phase offset
With a run of 200, I got these results:
old: 942.5619 +/- 0.0050
mine: 942.5618 +/- 0.0057
Anyway, you are looking to a resolution within 10% of the bin width. If you are looking that accurate, better use interpolation algorithms that also incorporate the phase.
Even without phase, for this kind of resolution you might want to include the window effects/scalloping loss in your frequency estimation.
Edit: Just did some calculations, using your equation, based on a Blackman-Harris window, you expect to measure 941.82 Hz

Interesting, running off internal sine wave at 941 Hz your library returns 942.0 Hz with streaming input about every 12 millisecond.

Any thoughts on improving the algorithm? Using a different window?
 
For the Blackman-Harris window:
k = index of max
d=(y(k+1)-y(k-1))/(y(k-1)+y(k)+y(k+1));
f=(f_sample/N)*(k+d*11/8)*(N+1)/(N);
 
Last edited:
for the blackman-harris window:
K = index of max
d=(y(k+1)-y(k-1))/(y(k-1)+y(k)+y(k+1));
f=(f_sample/n)*(k+d*11/8)*(n+1)/(n);


Did I cipher this correctly? Original code 942.0 using the above 948.77

Code:
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>

double f_sample= AUDIO_SAMPLE_RATE_EXACT,  NN = 1024.0, Freq;
int maxbin, k;

const int myInput = AUDIO_INPUT_MIC;

// Create the Audio components.  These should be created in the
// order data flows, inputs/sources -> processing -> outputs
//
AudioInputI2S       audioInput;         // audio shield: mic or line-in
AudioSynthWaveformSine  T941;
AudioAnalyzeFFT1024  myFFT;
AudioOutputI2S      audioOutput;        // audio shield: headphones & line-out

// Create Audio connections between the components


AudioConnection c2(T941, 0, myFFT, 0);


// Create an object to control the audio shield.
 
AudioControlSGTL5000 audioShield;


void setup() {
  // Audio connections require memory to work.  For more
  // detailed information, see the MemoryAndCpuUsage example
  AudioMemory(12);
  myFFT.windowFunction(AudioWindowBlackmanHarris1024);
  
  // Enable the audio shield and set the output volume.
  audioShield.enable();
  audioShield.inputSelect(myInput);
  audioShield.volume(0.6);
  T941.amplitude(0.9);
   T941.frequency(941.0);

   
}

void loop() {
  if (myFFT.available()) {
    // each time new FFT data is available
    // print it all to the Arduino Serial Monitor
     k=0;
     maxbin=0;
    
    
    for (int i=0; i<512; i++) {
      if (myFFT.output[i]>maxbin) {
        maxbin=myFFT.output[i];
        k=i;
      } 
     
    }
    Serial.println();
    Serial.print("Max Bin ");
    Serial.print(k);
    Serial.println();
    double y1=double(myFFT.output[k-1]);
    double y2=double(myFFT.output[k]);
    double y3=double(myFFT.output[k+1]);
   //double d=(y3-y1)/((2.0*(2.0*y2-y1-y3)));
    double d=(myFFT.output[k+1]-myFFT.output[k-1])/(myFFT.output[k-1]+myFFT.output[k]+myFFT.output[k+1]);
   
    Serial.print("Diff ");
    Serial.print(d);
    Serial.println();
  // double Freq=(AUDIO_SAMPLE_RATE_EXACT/1024.0)*(double(k)+d);
   
    Freq=(f_sample/NN)*(k+d*11/8)*(NN+1)/(NN);
    Serial.print("Freq ");
    Serial.print(Freq);
    Serial.println();
  
  }
}
 
kpc: Have you gone to the new TeensyD 1.21-b7? Can you post a copy of your current .cpp/.h? In hand editing then backing out and cloning the tree it would make sure I'm testing what you've got.

edit: Since I put in your un-un-rolled code I don't have the cases after 7. A fresh copy would be best - leaving the twiddle init in sketch setup()

edit: I calculated a correlated Magnitude [sqrt(sum(x^2))] like this:

int kk;
for ( ii ...) {
kk = 100 * myFFT.read(ii);
bSumXsq += kk * kk;
}
kk = sqrt( bSumXsq);

edit: I collected this summary data on the first peak [some index errors still]: Peak_Val, Peak_Bin; Start_bin, AscDuration, (Peak), DecDuration, SecondPeak_Seen:
FOR:: 59 ]FFT:15 21 25 20 11 4 2 2 1 1 1
bMax=25 bMaxId=2 bStart=1 bAscDur=1 bDecDur=7 bSecond= 32
 
Last edited:
@cartere: The problem was twofold. First I accidentally used another window, that I did not restore. The other is that in your code the division should be forced to be double.
double a = (int)/(int). will first do an integer division.
The other problem is that due to the window you need some other correction coefficients

Code:
double d=((double) (myFFT.output[k+1]-myFFT.output[k-1]))/(myFFT.output[k-1]+myFFT.output[k]+myFFT.output[k+1]);
Freq=(f_sample/NN)*(k+d*2.251);
This should result in an error of at maximum 0.1 Hz for a 1024 point fft, given a single frequency as input.
 
Last edited:
@cartere: The problem was twofold. First I accidentally used another window, that I did not restore. The other is that in your code the division should be forced to be double.
double a = (int)/(int). will first do an integer division.
The other problem is that due to the window you need some other correction coefficients

Code:
double d=((double) (myFFT.output[k+1]-myFFT.output[k-1]))/(myFFT.output[k-1]+myFFT.output[k]+myFFT.output[k+1]);
Freq=(f_sample/NN)*(k+d*2.251);
This should result in an error of at maximum 0.1 Hz for a 1024 point fft, given a single frequency as input.

Wonderful, right on the mark.

Code:
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>

double f_sample= AUDIO_SAMPLE_RATE_EXACT,  NN = 1024.0, Freq;
int maxbin, k;

const int myInput = AUDIO_INPUT_MIC;

// Create the Audio components.  These should be created in the
// order data flows, inputs/sources -> processing -> outputs
//
AudioInputI2S        audioInput;
AudioSynthWaveformSine   sine1;          //xy=236,270
AudioAnalyzeFFT1024      myFFT;      //xy=650,312
AudioConnection          patchCord1(sine1, 0, myFFT, 0);
// Create an object to control the audio shield.
AudioControlSGTL5000 audioShield;

void setup() {
  // Audio connections require memory to work.  For more
  // detailed information, see the MemoryAndCpuUsage example
  AudioMemory(12);
  myFFT.windowFunction(AudioWindowBlackmanHarris1024);
  audioShield.enable();
  audioShield.inputSelect(myInput);
  audioShield.volume(0.6);
  sine1.amplitude(0.7);
  sine1.frequency(941.0);
 }

void loop() {
  if (myFFT.available()) {
    // each time new FFT data is available
     k=0;
     maxbin=0;
     for (int i=0; i<512; i++) {
      if (myFFT.output[i]>maxbin) {
        maxbin=myFFT.output[i];
        k=i;
      } 
     
    }
    double d=((double) (myFFT.output[k+1]-myFFT.output[k-1]))/(myFFT.output[k-1]+myFFT.output[k]+myFFT.output[k+1]);
    Freq=(f_sample/NN)*(k+d*2.251);
    Serial.print("Freq ");
    Serial.print((Freq));
    Serial.println();
  }
}

Produced this output:
Freq 940.90
Freq 940.90
Freq 940.91
Freq 940.92
Freq 940.91
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.91
Freq 940.90
Freq 940.92
Freq 940.90
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.92
Freq 940.91
Freq 940.91
Freq 940.92
Freq 940.91
Freq 940.92
Freq 940.91
Freq 940.90
Freq 940.93
Freq 940.91
 
@kpc - I'll move to that and take out the sketch twiddle init. Also add the peak bin FREQ calc too to my loop to give it more work and data to compare the old and new since it should ID the sin FREQ as it is used.
 
We are now back to the original problem, but with much more stability.
With an internal 941 Hz sine both your latest library and Paul's original library return 941 Hz.
With a tuning fork counted with an external freq counter at 942.9 Hz using the mic input your latest library returns 939 Hz while Paul's original library returns 943 Hz.
To remove any possible DC interaction I ran the input through a biquad configured for highpass Fc at 800Hz Q=5 , with the same results.
All this using the latest frequency estimation code. Again thanks for your help with that :)
 
@carrera: This piece of code constantly dumps blocks of 1024 samples.
Could you record one block for me, with the tuning fork data in it?
Each block contains 64 lines with 16 samples each and is separated by whitelines from the other blocks.

Code:
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>

AudioInputI2S i2sin;
//AudioSynthWaveformSine sine1;
AudioRecordQueue queue1;
AudioConnection patchCord1(i2sin, 0, queue1, 0);
AudioControlSGTL5000 audioShield;

void setup() {
  AudioMemory(12);
  audioShield.enable();
  audioShield.inputSelect(AUDIO_INPUT_MIC);
  audioShield.volume(0.6);
//  sine1.amplitude(0.7);
//  sine1.frequency(941.0);
  queue1.begin();
 }

void loop() {
  if (queue1.available() > 1024 / AUDIO_BLOCK_SAMPLES) {
    for (int i = 0; i < 1024 / AUDIO_BLOCK_SAMPLES; i++) {
      int16_t *buf = queue1.readBuffer();
      for (int s = 0; s < AUDIO_BLOCK_SAMPLES; s++) {
        Serial.print(buf[s]);
        Serial.print((s % 16 == 15) ? "\n" : " ");
      }
      queue1.freeBuffer();
    }
    Serial.println();
    queue1.clear();
  }
}
 
Last edited:
@defragster: Frequency detection will only work for every third fft in your code. Since the fft is overlapping, you will always have two ffts which have spectral widening due to the frequency transition. Only one in three will have a pure sine wave. In your code if you set the count per frequnecy higher you will see that the number of "good" fft's will increase, with always two fft's "bad" at the transition. I hope the following diagram will make things clear. So to obtain a good spectrum and therefore a good frequency estimation the fft should not cross a border between your input frequencies.
Code:
      0000001111111111112222222222
fft0: xxxxxxxx          |
fft1:     xxxxxxxx      |
fft2:       | oooooooo  |
fft3:       |     xxxxxxxx
fft4:       |         xxxxxxxx
fft5:       |           | oooooooo
The diagram show horizontally a time axis with each charactr representing a block of 128 samples, vertically the fft's. The 0,1,2 at the top there indicates that the block contains the first sine wave, second etc. The xxx and ooo show which eight blocks are used for the fft. the xxx resulting in a 'bad' spectrum and the ooo in a 'good' spectrum.
 
Last edited:
Yes, I expected/assumed as much. While the freq is transitioning over sample groups it will have no true reference.

I'm doing my best to emulate the posted detect code and still looking for my error in adjusting the altered variables names and data source points.

Not knowing the Audio structures underlying the sampling I'm seeing that my FFT INO from the Audio Examples used .read() and the detect calc as shown uses .output[] and wondering if :
myFFT.output[bMaxId - 1] is the same as myFFT.read(bMaxId - 1)
 
Read also tries to integrate and scale
If you ask for one bin at the time, apart from the scaling, both should be equal.
The read code also has the problem of adding voltages, whereas it should be adding powers (sqrt(sum(x^2))
Code:
       float read(unsigned int binFirst, unsigned int binLast) {
                if (binFirst > binLast) {
                        unsigned int tmp = binLast;
                        binLast = binFirst;
                        binFirst = tmp;
                }
                if (binFirst > 511) return 0.0;
                if (binLast > 511) binLast = 511;
                uint32_t sum = 0;
                do {
                        sum += output[binFirst++];
                } while (binFirst < binLast);
                return (float)sum * (1.0 / 16384.0);
        }
 
Last edited:
Not sure if this leaves out any unclear variables:

Code:
    for (ii = 0; ii < 200; ii++) {
      kk = 100 * myFFT.read(ii);
      if ( kk ) {
        bSumXsq += kk * kk;
        if (kk > bMax) {
          bMax = kk;
          bMaxId = ii;
          if ( -1 == bStart )  bStart = ii;  // Start on first non-zero
        }
        if (0 != bDecDur )  bSecond += 1;  // detect multiple peaks
        iLastVal = ii;
      }
      else if (0 != bMax && 0 == bDecDur) {
        bAscDur = bMaxId - bStart;  // Duration before Max
        bDecDur = ii - bMaxId - 2;   // Duration after Max
      }
    }

    // ----------------------- CalcFreq
    double dd = ((double) (myFFT.output[bMaxId + 1] - myFFT.output[bMaxId - 1])) / (myFFT.output[bMaxId - 1] + myFFT.output[bMaxId] + myFFT.output[bMaxId + 1]);
    CalcFreq = (AUDIO_SAMPLE_RATE_EXACT / Freq) * (bMaxId + dd * 2.251);
 
My code following the FFT example looked at .read() - and I tried both .read() and .output[] on the on the freq detect.

Just reading one bin at a time - so I should use all of one or all of another? I suppose the 'units' (scaling) cancels as I get the same trend in the provided code. My freq is always near 1,000.

I just adjusted my freq jump count so I'd have three settled sample sets before changing and all three, as well as the mixed sample set always come in near 1,000 - I must have a bad reference value in detect.
 
Sorry, I am not sure I understand the question in post #44.
as for post#45, scalling should indeed not matter.
 
Last edited:
Is there an obvious error in the code snipped above (msg #44) - or a question about it's value? My Freq Detect Calc is not working. I'm still on PJRC FFT code until I see this working right, so it isn't a problem in your code - just getting a baseline.

[edit] - I need to switch to blackmanharris. And I don't see where the 'n' is used in code sample above code from in msg #35:
for the blackman-harris window:
K = index of max
d=(y(k+1)-y(k-1))/(y(k-1)+y(k)+y(k+1));
f=(f_sample/n)*(k+d*11/8)*(n+1)/(n);
 
Last edited:
The problem is this line, the first division should be by the number of samples.
CalcFreq = (AUDIO_SAMPLE_RATE_EXACT / 1024.) * (bMaxId + dd * 2.251);
 
Status
Not open for further replies.
Back
Top