Go Down

Topic: Reliable Frequency Measurment Using Autocorrelation (Read 30864 times) previous topic - next topic

jremington

#15
Oct 09, 2014, 09:28 pm Last Edit: Oct 10, 2014, 06:10 pm by jremington Reason: 1
Not my code, but you have a point. It is certainly easy to overlook those traps in someone else's!

Curiously, using Code::Blocks/gcc on the PC as I was, it does not make any difference to the overall result to cast rawData[] as long, i.e.
Code: [Select]
sum += ( (long) rawData[k]) * ( (long) rawData[k+i]);
It appears that on the PC, Code::Blocks promotes the type char to a larger integer for the purpose of the multiplication. I haven't yet figured out how to view the machine code listing to check this supposition, so I printed out the first few values of rawData[]*rawData[] using the following statement:
Code: [Select]
   fprintf(fp,"%d,%d,%d\n",i,rawData[i],rawData[i]*rawData[i]);and got the following values:
Quote

0,115,13225
1,126,15876
2,62,3844
3,57,3249
4,-80,6400
5,77,5929

Clearly, the multiplication is not overflowing the char data type. So, I don't see why Matlab would give a "much different result" for the autocorrelation function (why not show us what it is?).



jremington

I checked the operations performed by Matlab and found that the built-in function autocorr() calculates the same unweighted function that the OP originally described (except that the mean value of the data is properly subtracted). This is fine for short lags but reduces the weight of peaks associated with large lags. From the autocorr help message "the form adopted here follows that of Box, Jenkins, and Reinsel  ... Time Series Analysis: Forecasting and Control. 3rd ed. Upper Saddle River, NJ: Prentice-Hall, 1994". I've attached a sample plot obtained by using the following Matlab script:
Code: [Select]

m=csvread('test.txt');
ACF=autocorr(m(:,2),969);
plot(ACF);


One can see how the peaks diminish as the lag increases, which I consider to be misleading. The Fourier transform of this particular autocorrelation function will not be the same as the power spectrum plot of the zero-filled FFT that I posted earlier.

pantherturm

Hello!

I've tested the code of robtillaart but i'm getting some weird characters on the serial monitor. Any idea why is this happening?

robtillaart

Rob Tillaart

Nederlandse sectie - http://arduino.cc/forum/index.php/board,77.0.html -
(Please do not PM for private consultancy)

pantherturm

I'm using the Uno board. I have uploaded the following code...


Code: [Select]


#include "sample.h"

void setup()
{
  Serial.begin(115200);
  Serial.println("FrequencyDSP");
  Serial.println("------------");

  uint32_t start = millis();
  float fr = autoCorrelateFrequency(rawData, sizeof(rawData), 22050);
  uint32_t stop = millis();

  Serial.print("freq:\t");
  Serial.println(fr);
  Serial.print("time:\t");
  Serial.println(stop - start);
}

void loop()
{
}

float autoCorrelateFrequency(char * sample, int len, float sampleFreq)
{
  long sum = 0;
  long sum_old = 0;
  int thresh = 0;
  byte pd_state = 0;
  int period = 0;  // 0 results in inf

  // Autocorrelation
  for (int i=0; i < len && (pd_state != 3); i++)
  {
    sum_old = sum;
    sum = 0;

    for (int k=0; k <len-i; k++)
    {
      sum += (rawData[k]) * (rawData[k+i]);
    }
    sum /= 256;

    // Peak Detect State Machine
    // 0: initial, set threshold
    // 1: detect positive slope
    // 2: detect negative slope
    // 3: peak found
    if (pd_state == 0)
    {
      thresh = sum / 2;
      pd_state = 1;
    }
    else if (pd_state == 1 && (sum > thresh) && (sum - sum_old) > 0)
    {
      pd_state = 2;
    }
    else if (pd_state == 2 && (sum - sum_old) <= 0)
    {
      period = i;
      pd_state = 3;
    }
  }

  return sampleFreq/period;
}


but instead of a frequency what i'm getting on the serial monitor is this

Quote

üÕ&
K±μ
Yþ)¤

Sorry i'm not an experienced arduino user not sure what else i should provide.

wildbill

Check that you have set the serial monitor's baud rate to match that in your sketch (115200).

pantherturm


michinyon

Quote
When you perform the multiplication in the step which multiplies the two rawData values you are overflowing the type char.  The multiplication is done on type char before it is assigned to variable sum.
I am unconvinced by this.   The multiplication of 2 single-byte numbers on the Arduino results in a 16-bit result.

michinyon

The other problem with the autoCorrelateFrequency( )  function in Reply # 1,  is that this function gets passed a pointer to the data array,  as a function argument,   and then doesn't use it -  it goes straight to the global data.

This problem becomes more obvious if you put each function in a separate file.

michinyon

Quote
I get the freq = 5512.5Hz result,
So do I.

I tried to run this program using Netbeans using cygwin/gcc on a windows computer,  and had all sorts of grief.

I had to define the data as unsigned char to get it to compile ,    and then I got 5512 Hz

But if I defined the data as char,  I got 259 Hz.

And if I defined the data as signed char,   I also got 259 Hz.

I believe there is an "issue" with C, on different platforms,  concerning whether char numbers are signed or unsigned if not explicitly stated.   That is probably where your 5512 Hz comes from.






nenea_dani

Hi Akellyirl,
In the version that collects data in A0 why just at low frequencies the result is correct? For 1000hz real they get 891hz and for 2000 they get 1783hz. Is the problem in autocorrelation or in machine state? I used http://www.akellyirl.com/arduino-guitar-tuner/

MrMark

Hi Akellyirl,
In the version that collects data in A0 why just at low frequencies the result is correct? For 1000hz real they get 891hz and for 2000 they get 1783hz. Is the problem in autocorrelation or in machine state? I used http://www.akellyirl.com/arduino-guitar-tuner/
The original post is several years old, so I'll be surprised if the original poster responds.

The issue is that the measured frequency resolution decreases as the frequency increases.  This is because the autocorrelation peaks are measured at discrete steps.  Resolution is approximately f2/fs where "f" is the frequency being measured and "fs is the sample rate.  So at 2000 Hz with a sample rate of 8919 Hz, the resolution is only about 448 Hz.

A brute force way of getting better resolution is to use a higher sampling rate, but RAM to store the sample then becomes an issue.  A more sophisticated approach would be to interpolate between correlation points by using some sort of curve fit around the discrete peak.

nenea_dani

I migrated from a FHT frequency meter for 150-4000Hz that has an error of less than 10 percent in the idea that I can measure it more precisely.  The author said he can get even decimal but not specified as just for two hundred or three hundred Hz. Looks like I'm out of math. The formula " f square / fs"  " is very weird for me and I still want to know where it comes from unless I ask too much.  Do I understand that autocorrelation is convenient to implement in general only for low frequencies ?  Anyway thanks for the information. For now I understand that I am in a direction that surpasses me and I return to FHT that even if I do not understand how it works, at least it gives very good results.

tim77777

I found two major problems with the changes made to akellyirl's original code.

1. Although the changes work OK on the supplied data file, when applied to a changing input such as from an A/D it becomes obvious that often a frequency exactly half of the actual, is being outputted.

2. Changing the raw data from 'byte' to 'signed char', increases the overall processing time for each outputted frequency from 100mS to about 160mS.

The modifications involved a lot of formatting/structural changes from the original code (added function etc), so maybe I messed it up somewhere.  This is how I combined the changes with Akellyirl's original code.

Code: [Select]

#define LENGTH 256

signed char rawData[LENGTH];


// Sample Frequency in kHz
const float sample_freq = 8919;

int len = sizeof(rawData);
int i,k;
long sum, sum_old;
int thresh = 0;
float freq_per = 0;
byte pd_state = 0;
int count;


void setup()
{
  analogReference(DEFAULT);   // Connect to 3.3V
  analogRead(A0);
  Serial.begin(115200);
  count = 0;
}


void loop(){

    if (count < LENGTH)
    {   
      rawData[count] = analogRead(A0)>>2; 
      count++;
    }
 
    else
    {
       sum = 0;
       pd_state = 0;
       int period = 0;
       
       for(i=0; i < len; i++)
       {       
            sum_old = sum;
            sum = 0;
           
            for (k=0; k <len-i; k++)
            {
               sum += (rawData[k]) * (rawData[k+i]);
            }
            sum /= k ;
               
            if (pd_state == 0)
            {
              thresh = sum / 2;
              pd_state = 1;
            }
            else if (pd_state == 1 && (sum > thresh) && (sum - sum_old) > 0)
            {
              pd_state = 2;
            }
            else if (pd_state == 2 && (sum - sum_old) <= 0)
            {
              period = i;
              pd_state = 3;
            }
 
       }
 
        if (thresh >100 && period > 0)
        {
          freq_per = sample_freq/period;
          Serial.println(freq_per);
        }
     
        count = 0;
   
    }
 
}


MarkT

The accuracy of frequency samples without interpolation should be fs/N where N is the FFT size.

However in the time domain the time error is 1/fs, which is f/fs wave periods, so the frequency error is about
f * f/fs, or f^2/fs.

Once you have an estimate of the frequency you can use other methods to refine the estimate, perhaps as
simple as crude band-pass filtering and zero-crossing counting.
[ I will NOT respond to personal messages, I WILL delete them, use the forum please ]

Go Up