Pages: 1 [2]   Go Down
Author Topic: kurtosis calculations on an analog signal  (Read 985 times)
0 Members and 1 Guest are viewing this topic.
Offline Offline
Newbie
*
Karma: 0
Posts: 12
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Yes. It measures the peakedness of a signal and can be defined as either the fourth cumulant divided by the second cumulant or the fourth moment around the mean divided by the square of the variance.

We are using and EMG/EOG test to score REM Sleep and possibly be a prediagnostic for early neurological disorders.

http://petrsu.ru/Chairs/Physiology/JEK.pdf

Kurtosis is one of the needed extracted features.

Thanks.
Logged

Offline Offline
Faraday Member
**
Karma: 62
Posts: 3077
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

OK,  so they are not sinusoidal signals they are measuring.     And,  kurtosis does not measure the peakiness of the signal,   it measures the peakiness of a histogram showing the distribution of values in the sample.    Not the same thing at all.

So,  you have your sampled signal.   You have your equation defining the kurtosis of the sample numbers.  Go for it.    Whether that result is at all meaningful,   would be for others to decide.

From the arduino point of view,  your biggest obstacle is going to be shortage of memory to contain all the data.
Logged

Oregon, USA
Offline Offline
Edison Member
*
Karma: 69
Posts: 2381
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

For those interested in the theory, a concise and well written summary of statistical analyses of a time series of measurements can be found in the "Handbook of Human Vibration".  The kurtosis is defined and a couple of key results are that the kurtosis of a sine wave is 1.5, while the kurtosis of random, i.e. Gaussian-distributed, signals is about 3. The technique can be used for discriminating signal-containing transmissions from noise in radio communications. The formulas would be trivial to implement if the memory to contain the time series were on hand.

Courtesy of Google Books, here is the relevant page: http://books.google.com/books?id=-_JV63L6RMcC&pg=PA468&dq=handbook+human+vibration+kurtosis&hl=en&sa=X&ei=dkRtUsncLqmajALd3oHICQ&ved=0CC4Q6AEwAA#v=onepage&q=handbook%20human%20vibration%20kurtosis&f=false

An Arduino would be able to handle a few hundred data points, but the question is: would that be enough to give statistically significant results for low-frequency signals like those from electrophysiology measurements?


Logged

"It seems to run on some form of electricity"

Offline Offline
Newbie
*
Karma: 0
Posts: 12
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

This is the code I wrote to find the kurtosis of an analog input at pin A0. Does anyone see any problems with this code? I want it to output a value with two decimal places on the serial monitor but for some reason it is not, do you see why?

const int numReadings = 500.00;
float inputPin = A0;
int kurtindex= 0.00;
int kurtVoltages[numReadings];
int total = 0.00;
int average = 0.00;
int sqDevSum = 0.00;
int sum= 0.00;
int kurtosis=0.00;
int stDev4=0.00;
int denominator = 0.00;
int numerator = 0.00;

void setup() {
  Serial.begin(9600);   
}

void loop() {
  kurtVoltages[kurtindex] = analogRead(inputPin);
 
  if (kurtindex >= numReadings);
    kurtindex =0.00; 
   
for (int i=0; i < numReadings; i++) {
    total += kurtVoltages;
  }
 
  average = total / numReadings;
 
 for (int i=0; i < numReadings; i++) {
    kurtVoltages = kurtVoltages - average;
  }
 
 for (int i=0; i < numReadings; i++) {
    sqDevSum += pow(kurtVoltages, 2);
  }
 
int stDev = sqrt(sqDevSum/numReadings);
int stDev4= pow(stDev, 4);
 
  for (int i=0; i < numReadings; i++) {
    kurtVoltages = pow(kurtVoltages, 4);
  }
 
  for (int i=0; i < numReadings; i++) {
   sum += kurtVoltages;
  }

denominator = numReadings*stDev4;
numerator = sum;
kurtosis = numerator/denominator;
kurtosis = kurtosis -3;

Serial.println(kurtosis);

}
Logged

USA
Offline Offline
Sr. Member
****
Karma: 17
Posts: 387
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

It's good to see the original poster put up some code.  It suggests that he's actually interested in a solution.

The posted code doesn't compile.  There are repeated references to kurtvoltages, as if it were a variable. But kurtvoltages is an array, and, if it's going to be treated like a single value, it needs an index.

There are logic errors, too, but they can wait.  If it doesn't compile, you'll never get any results at all.

[Edit: Add this]
And, until you know what you're doing, please don't try to hook up an analog source - like a signal generator or a biological monitor - to your Arduino.  It can be quite easy to damage an analog input when the input voltage goes out of bounds.  I suspect that your signal source is centered at zero volts, and that it's negative for at least some of the time.  A negative voltage on an analog pin can damage it.

This sounds like an academic assignment.  Is it?
« Last Edit: October 28, 2013, 05:42:48 pm by tmd3 » Logged

Offline Offline
Edison Member
*
Karma: 49
Posts: 1670
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

@blwill428
Read this and fix up your posted code: How to post code properly

Pete
Logged

Offline Offline
Newbie
*
Karma: 0
Posts: 12
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

I am sorry that I am new to this forum and to Arduino, so excuse me for my errors.

It is not for an academic assignment, but just a project some classmates and I are working on. With my old code, it did compile and upload to my Arduino Uno, but it was using 68% of the Dynamic Memory.

I wrote a new code using a Moving Averager Approach:

Code:
const int numberSamples = 5000;
float sensorPin = A3;
float alpha = 0.0002;
float beta = 0.9998;
float runningaverage= 0.000;
float subtract = 0.000;
float sqsubtract = 0.000;
float qudsubtract = 0.000;
float variance = 0.000;
float sigma4 = 0.000;
float fourthmoment = 0.000;
float denominator = 0.000;
float kurtosis = 0.000;
float xi;
int i;
long previousMillis = 0;
long interval = 1;   



void setup() {
  Serial.begin(9600);

}

void loop() {
  unsigned long currentMillis = millis();
  if(currentMillis - previousMillis > interval) {
    xi = analogRead(sensorPin);
    previousMillis = currentMillis; 
  }
  i=i+1;
  runningaverage = beta*runningaverage + alpha*xi;
  subtract = xi - runningaverage;
  sqsubtract = pow(subtract, 2);
  qudsubtract = pow(subtract, 4);
  variance = beta*variance + alpha*sqsubtract;
  sigma4 = pow(variance, 2);
  fourthmoment = beta*fourthmoment + alpha*qudsubtract;
  denominator = sigma4*numberSamples;
  kurtosis = fourthmoment/denominator;
 
 
  if (i>=numberSamples)
    Serial.println(kurtosis);
  if (i>=numberSamples)
    i=0;
 
 
}

This only uses 12% of the memory since it is not holding so much in an array.

I understand that negative voltage will harm the Arduino. I am using a DC Offset on a function generator to remedy this. When I move it to a biosignal, I will use a summing amplifier to make sure the bipolar signal is converted to unipolar. (http://masteringelectronicsdesign.com/measure-a-bipolar-signal-with-an-arduino-board/)

Logged

Seattle, WA USA
Offline Offline
Brattain Member
*****
Karma: 638
Posts: 50304
Seattle, WA USA
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Code:
float sensorPin = A3;
So you can read from pin 5.7 later, right? And then pin 0.001 and pin -14.7 and 3.14159.
Logged

USA
Offline Offline
Sr. Member
****
Karma: 17
Posts: 387
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Ingenious.  That's a lot of reasonable-looking code to come from a guy who didn't know where to start a few days ago.

Here are some suggestions:
  • [Edit:  Add this]
    Something's wrong with the timing and program flow.  The code decides when to take a sample, but it processes a sample regardless of whether or not it takes one.  Examine it:  loops starts; the test passes, it takes a sample, and processes it; loop terminates, and starts over; the test fails, and it doesn't take a sample, but it processes the last sample again anyway.  Most of that code needs to move into the "if(currentMillis ..." section, and the rest should only be executed when it's time to do final calculations and print.
  • You're trying to calculate a function of a series of values, and you don't yet know whether the calculation is correct.  You also don't know what the values are, since they're analog readings, and they aren't displayed.  Troubleshooting the calculation is difficult enough without having to wonder what the input data was.  While you're troubleshooting the calculation, use some convenient calculated data whose kurtosis is well-known, instead of analog readings.  Shutting off sample timing during troubleshooting will make the results show up faster.  Here are two suggestions:
Code:
xi = i  % 10
          // uniform distribution, expected excess kurtosis = -1.2
xi = sin(i*0.001*2*3.1415927*100) + 1
         // sine function, 100 Hz, sampled at 1/millisecond, expected excess kurtosis = -1.5

  • Decide whether you want to calculate kurtosis the old-fashioned way, or whether you want to calculate excess kurtosis.
  • Print some debugging information, like the value of the running average, fourth moment, and variance.  Be prepared to get deeper into the calculations if those values don't help.
  • Your "running average" technique, where newValue = alpha*sampleValue + (1-alpha)*oldValue, is a single-pole digital low-pass filter.  Its output value tends toward the mean of the input signal.  When you filter variance, you get something that looks like the statistical variance - the mean of the squares of the differences between the samples and the mean of the signal.  When you filter the fourth moment, you get something that looks like the fourth moment - the mean of the fourth power of those same differences.  Division by the number of samples is essentially accomplished in filtering those values.   Consequently, I'll recommend that you reconsider whether this line needs the multiplication by the number of samples:
Code:
 denominator = sigma4*numberSamples;
  • When you start getting output data, experiment with larger values of alpha.  I don't think it's necessary that alpha be the same as 1.0/numberSamples.  You might want to define beta  = 1.0 - alpha.
  • This looks to me to be an ingenious solution.  My experience with ingenious solutions to statistical calculations is that they often yield unpextected results, usually because there's some nugget of probability theory that I've overlooked that makes them invalid.  I can't see the flaw here, though - the level of filtering is pretty intense, and the signal mean will get pretty close to the actual value; once it stabilizes, I think that the variance and fourth moment calculations will be valid.  Nevertheless, I'd recommend testing it with several kinds of distributions to see that it gives the expected results, and be prepared for disappointment.
  • After you're confident of the kurtosis calculation, think about ways to mange the calculation of the signal mean.  When I run this code with known input data and print the value of runningaverage, it takes a while to stabilize.  Until it does, the output data will not be valid.  If the DC offset of the signal won't change during the measurement, consider instead taking some time to determine the signal mean while outputting no values, and then using that mean to calculate the second and fourth moments.  There's no loss of performance - there's a delay between startup and valid data either way.  With the mean known, the calculations are pretty simple.  Once the signal mean is established, I think that there's no necessity to filter the second and fourth moments - the calculation can be straightforward from the formula.
  • Consider AC-coupling the input signal.  The signal mean will then be the value of the analog offset voltage.  That'll give control of the signal mean to you, rather than to some setting on a sensor amplifier.  There are a lot of options for how to manage getting that voltage into the calculation, some involving hardware, some not.  A simple one would be to use a resistor divider to set it at half scale, and use that value as the initial value of runningaverage, rather than starting with zero.  That would reduce the time it takes for data to become valid, though possibly not by all that much - the second and fourth moments are still filtered.
  • As PaulS suggests, don't make sensorPin a float.  It seems to work, but it's just too weird.  sensorPin doesn't change, so it could easily be const int.
« Last Edit: October 29, 2013, 09:32:40 am by tmd3 » Logged

Offline Offline
Newbie
*
Karma: 0
Posts: 12
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Thank you tmd3 for being so understanding to my beginner's approach to Arduino. It is very intimidating in these forums and its nice to get a kind and constructive response. I appreciate your feedback and will work on it immediately.  I will let you know of my results.
Logged

USA
Offline Offline
Sr. Member
****
Karma: 17
Posts: 387
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Don't forget:  the first goal is to get the calculation working.  It's a whole lot easier to modify a working program than to start over with something new.  Fix this one, and worry about optimization and and improvements later.
Logged

Offline Offline
Newbie
*
Karma: 0
Posts: 12
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Here's the code I am using. It was working, outputting the right kurtosis value every 5 seconds, but now its not. Its outputting at a rapid rate that is definitely not 1 every 5 seconds. What went wrong?

Code:
#include <SoftwareSerial.h>
SoftwareSerial xbee(2, 3);

 

const int numberSamples = 5000;
float sensorPin = A2;
float alpha = 0.0002;
float beta = 1-alpha;
float runningaverage= 0.000;
float subtract = 0.000;
float sqsubtract = 0.000;
float qudsubtract = 0.000;
float variance = 0.000;
float sigma4 = 0.000;
float fourthmoment = 0.000;
float denominator = 0.000;
float kurtosis = 0.000;
float xi;
int i=0;

long previousMillis = 0;
long interval = 1;


void setup() {
Serial.begin(19200);
xbee.begin(19200);

}

void loop() {

i=i+1;

unsigned long currentMillis = millis();
if(currentMillis - previousMillis >= interval) {
previousMillis = currentMillis;
xi = i;
runningaverage = beta*runningaverage + alpha*xi;
subtract = xi - runningaverage;
sqsubtract = pow(subtract, 2);
qudsubtract = pow(subtract, 4);
variance = beta*variance + alpha*sqsubtract;
sigma4 = pow(variance, 2);
fourthmoment = beta*fourthmoment + alpha*qudsubtract;
denominator = sigma4*numberSamples;
kurtosis = fourthmoment/denominator;
}


if (i>=numberSamples)
Serial.println(kurtosis);
if (i>=numberSamples)
xbee.print(kurtosis);
if (i>=numberSamples)
i=0;


}
Logged

USA
Offline Offline
Sr. Member
****
Karma: 17
Posts: 387
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Every time loop() executes, the counter variable - i - is incremented, whether it's time to take a sample or not.  loop() executes many times between samples, so the counter reaches its limit quite quickly, and the output data prints.

Since the counter is supposed to be counting samples, it needs to be incremented only when a sample is taken.  It needs to be inside the brackets of the "if (currentMillis ... " test.

[Edit: Add this]  After a closer look:
  • Proper indentation will make your code easier to read and understand, for others and for you.  You can manage that yourself, or you can use "Auto Format" under the "Tools" menu in the Arduino IDE.
  • I recommend that you refrain from adding complexity to this code until you get the kurtosis calculation working.  Leave xbee stuff out of your program until you get kurtosis working.
  • Patience, Grasshopper.
« Last Edit: October 29, 2013, 10:27:47 pm by tmd3 » Logged

Seattle, WA USA
Offline Offline
Brattain Member
*****
Karma: 638
Posts: 50304
Seattle, WA USA
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

Code:
float sensorPin = A2;
Why?
Logged

Pages: 1 [2]   Go Up
Jump to: