Forum Rule: Always post complete source code & details to reproduce any issue!
Page 2 of 3 FirstFirst 1 2 3 LastLast
Results 26 to 50 of 67

Thread: 1024 point FFT, 30% more efficient (experimental)

  1. #26
    Senior Member
    Join Date
    Oct 2013
    Location
    Rogersville MO
    Posts
    253

    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();
      }
    }

  2. #27
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    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 by defragster; 03-02-2015 at 08:50 PM.

  3. #28
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    @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 by kpc; 03-02-2015 at 09:03 PM.

  4. #29
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    @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 by defragster; 03-02-2015 at 09:16 PM.

  5. #30
    Senior Member
    Join Date
    Oct 2013
    Location
    Rogersville MO
    Posts
    253
    Quote Originally Posted by kpc View Post
    @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?

  6. #31
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    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 by kpc; 03-02-2015 at 11:19 PM.

  7. #32
    Senior Member
    Join Date
    Oct 2013
    Location
    Rogersville MO
    Posts
    253
    Quote Originally Posted by kpc View Post
    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();
      
      }
    }

  8. #33
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    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 by defragster; 03-03-2015 at 11:18 AM.

  9. #34
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    @cartere, at first glance it looks good. For unwindowed fft's see http://www.dspguru.com/dsp/howtos/ho...olate-fft-peak

    @defragster, see attachment I fixed the twiddle array, by actually precalculating the values
    Edit: Oops uploaded wrong file. Proper file will folow.
    Last edited by kpc; 03-04-2015 at 05:37 PM.

  10. #35
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    @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 by kpc; 03-04-2015 at 06:13 PM.

  11. #36
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    @defragster: see attachment
    Attached Files Attached Files

  12. #37
    Senior Member
    Join Date
    Oct 2013
    Location
    Rogersville MO
    Posts
    253
    Quote Originally Posted by kpc View Post
    @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

  13. #38
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    @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.

  14. #39
    Senior Member
    Join Date
    Oct 2013
    Location
    Rogersville MO
    Posts
    253
    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

  15. #40
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    @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 by kpc; 03-05-2015 at 05:42 AM.

  16. #41
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    @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 by kpc; 03-05-2015 at 06:10 AM.

  17. #42
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    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)

  18. #43
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    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 by kpc; 03-05-2015 at 06:16 AM.

  19. #44
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    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);

  20. #45
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    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.

  21. #46
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    Sorry, I am not sure I understand the question in post #44.
    as for post#45, scalling should indeed not matter.
    Last edited by kpc; 03-05-2015 at 06:35 AM.

  22. #47
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    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 by defragster; 03-05-2015 at 06:41 AM.

  23. #48
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    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);

  24. #49
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,973
    Thanks, MY BAD: I misread this line
    Code:
    double f_sample= AUDIO_SAMPLE_RATE_EXACT,  NN = 1024.0, Freq;

  25. #50
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    As for the edit in #47, forget about n, this was a mistake from my side.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •