Go Down

Topic: Reliable Frequency Measurment Using Autocorrelation (Read 9580 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.

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)

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).


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.






Go Up