Data Analysis questions

I've made good progress on measurement of my tremors using a LSM6DS33 and ATmega2560. I can log both gyro and accelerometer at 250hz. Offline I managed to determine the frequency of my tremors and can apply that code to the Arduino to measure frequency in real time.

Now come the problems. The signal is a bit noisy so I had to use moving averages to smooth the curve but now I want to continue to measure my tremor when I move. Up until now, all measurements have be in a static pose. This going to require a High-pass digital filter that will remove all signal and noise in the range of 0-3 Hz. Nice to have would be to filter off anything over 20hz. Nice to have would be a power calculation of the 3-20hz band.

I've read a bit about digital filters and even tried some online design applications but so far have failed to get my head around digital filter design and implementation. Pointers to a "dummies start here" design / software would be appreciated.

Don't know much about the LSM6DS33 but I did see in Table 64 of the spec a register to enable a low pass filter LPF2 Enable. Have you been using that?

Thank you for your reply. I did have a look at these filters but didn't pursue them as they appear to lack the range that I require for my particular application.

But their like chicken soup..... couldn't hurt.

I don't know if this will be of much use but look here a FFT

The FFT (Fast Fourier Transform) converts data in real time to the frequency domain.
To apply the "Fast" in Fast Fourier Transform, the number of samples must be a factor of 2 (i.e. 2,4,8....).

Although its a bit more complex, the theory is: once you convert the data to the frequency domain you can remove the frequencies you don't want. Then convert the data back to the time domain.

Its been a long time since I worked with the FFT and not in a micro so I can't be of much more help than this. However I'll guess there is a lot of info and likely libraries on the internet.

Question: How does knowing the frequency become a benefit?

Good luck :slight_smile:

Long story short: I have Essential Tremors of both hands. My blood pressure medicine, a beta blocker, partially controls the amplitude of my tremors. There is a very expensive device, Cala Trio, that helps with tremors of the hand. It does so by applying an electrical stimulation to the Medial and Radial nerve of the wrist, alternating between the two nerves at the frequency of the tremor.

sooooooooo knowing the frequency of my hand tremor allows be to build a similar device from an Arduino plus some other COTS h/w for about 2% of the cost of the original.

The DFT (Discrete Fourier Transform) does that analysis. First, study the LSM6DS33 data sheet to learn how to turn on the appropriate low- and high-pass filter options.

Here is a simple, slow example of the DFT that demonstrates how to calculate the band-limited power spectrum.

float a[8] = {0, 0, 0, 0, 10, 10, 10, 10}; //example square wave signal, period = 8 x sample time

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

  //sample time = 10 ms, signal frequency in this case = 1/period = 12.5 Hz.
  //DFT calculated for 0 to 40 Hz in steps of 0.5 Hz, assuming signal repeats 4 times, no window
  // expect peaks at 0Hz (DC offset), 12.5 Hz, 37.5 Hz, ...

  dft (a, 8, 10, 0, 40, 0.5, 4, 0);

}

void loop() {
}

// from https://www.instructables.com/Arduino-Frequency-Transform-DFT/
// some errors corrected SJR 10/2021
/*
  8 TERMS THAT NEED TO BE SPECIFIED
  1: an array of which dft need to be taken
  2:size of the array
  3:time interval between 2 reading in array in milliSECONDS
  4:lower value of frequency range in Hz
  5:upper value of frequency range in Hz
  6:size of steps for frequency range
  7:repetitions of signal array (minimum 1) higher number better accuracy but increased solution time
  8:  0 for no window, 1 for flat-top window, 2 for Hann window, 3 for Hamming window
      (if you do not have any idea about selecting window keep default 3)
  example:
      dft(a,110,0.5,0,30,0.5,10,3);
      here a is an array of size 110 element to be checked for 0 Hz to 30 Hz with 0.5 step (0,0.5,1,1.5, ... ,29,29.5,30)
      10 repetition and hamming window
      by- ABHILASH PATEL
*/

float dft(float a[], int arraysize, float interval, float f0, float fn, float stepsize, int reps, int window)
{
  float mag, sumi, sumr, ti, tr;
  int j, k;
  float twopi = 2.0 * PI;
  Serial.print("data array repetitions = ");
  Serial.println(reps);

  // apply windowing function
  if (window == 1) //flat-top window
  {
    for (int i = 0; i < arraysize; i++)
    {
      float b = PI * i / (arraysize);
      a[i] = a[i] * ( 1 - (1.93 * cos(2 * b)) + (1.29 * cos(4 * b)) - (0.388 * cos(6 * b)) + (0.028 * cos(8 * b)));
      // Serial.println(a[i]);
    }
  }

  if (window == 2) //hann window
  {
    for (int i = 0; i < arraysize; i++)
    {
      float b = twopi * i / (arraysize);
      a[i] = a[i] * 0.5 * (1 - cos(b));
      //Serial.println(a[i]);
    }
  }

  if (window == 3) //hamming window
  {
    for (int i = 0; i < arraysize; i++)
    {
      float b = twopi * i / (arraysize);
      a[i] = a[i] * (0.54 - 0.46 * cos(b));
    //  Serial.println(a[i]);
    }
  }

  // DFT calculation

  for (float f = f0; f <= fn; f = f + stepsize)
  {

    sumi = sumr = 0.0;
    k = 0;
    for (int i = 0; i < (arraysize * reps); i++)
    {
      j = i - k;
      if (j >= arraysize) {
        k = k + arraysize;  //signal repeat index offset
      }
      ti = a[j] * (sin(twopi * f * i * interval * 0.001)); //0.001 => time in seconds
      tr = a[j] * (cos(twopi * f * i * interval * 0.001));
      sumi = sumi + ti;
      sumr = sumr + tr;
    }
    mag = sqrt(sumi * sumi + sumr * sumr) / (arraysize * reps);
    Serial.print(f);
    Serial.print("\t");
    Serial.println(mag);
  }

}

Neat ! Kind of like noise canceling headphones :slight_smile:

Thanks for the code. I think it be very slow on an ATmega2560. I was considering moving to one of the ARM based processor such as a Teensy 4.1 where its speed should be ok. I'll give it a try when I move processors.

Even the 'inventors' don't know exactly how the Trio helps to suppress the amplitude of ETs in the hand. Current speculation is that it somehow stimulates or distracts the ventral intermediate (VIM) nucleus of the thalamus of the brain. More research needs to be done before this can be established.

The Trio differs from noise canceling headphones in that you treat your tremors for 40 minutes and then you get upwards of 2-3 hours partial relief from your tremors. Somewhere between 40% and 50% of ET patients respond to this treatment and there is a large variance in the degree of attenuation of amplitude of tremor.

Depends on what you mean by "slow" -- try it. But from the application you describe, it is hard to see how speed could be an issue.

If needed, execution speed can be vastly increased by a few very simple tricks, like a table lookup for the sine and cosine functions. And unlike the FFT, you don't have to calculate the entire spectrum, only the segment of interest.

I think there is a problem with the above code section where you check if j is greater or equal to the arraysize and then a couple of lines later index the a[] array with j.

Doesn't j need to be reset if equal to or greater than the arraysize?

No, look how j is calculated. The code works fine, and gets the correct answer with test cases.

That is not my code. The download link is in the code.

Last time I checked 'C' arrays have a range of 0 to (arraysize -1). The size of a[] is 8. The valid indexes range from 0-7. a[8] is out of range. Please correct me if I'm wrong. thx

You are correct. The code is missing a line.

Corrected code:

float a[8] = {0, 0, 0, 0, 10, 10, 10, 10}; //example square wave signal, period = 8 x sample time

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

  //sample time = 10 ms, signal frequency in this case = 1/period = 12.5 Hz.
  //DFT calculated for 0 to 40 Hz in steps of 0.5 Hz, assuming signal repeats 4 times, no window
  // expect peaks at 0Hz (DC offset), 12.5 Hz, 37.5 Hz, ...

  dft (a, 8, 10, 0, 40, 0.5, 4, 0);

}

void loop() {
}

// from https://www.instructables.com/Arduino-Frequency-Transform-DFT/
// some errors corrected SJR 10/2021
/*
  8 TERMS THAT NEED TO BE SPECIFIED
  1: an array of which dft need to be taken
  2:size of the array
  3:time interval between 2 reading in array in milliSECONDS
  4:lower value of frequency range in Hz
  5:upper value of frequency range in Hz
  6:size of steps for frequency range
  7:repetitions of signal array (minimum 1) higher number better accuracy but increased solution time
  8:  0 for no window, 1 for flat-top window, 2 for Hann window, 3 for Hamming window
      (if you do not have any idea about selecting window keep default 3)
  example:
      dft(a,110,0.5,0,30,0.5,10,3);
      here a is an array of size 110 element to be checked for 0 Hz to 30 Hz with 0.5 step (0,0.5,1,1.5, ... ,29,29.5,30)
      10 repetition and hamming window
      by- ABHILASH PATEL
*/

void dft(float a[], int arraysize, float interval, float f0, float fn, float stepsize, int reps, int window)
{
  float mag, sumi, sumr, ti, tr;
  int j, k;
  float twopi = 2.0 * PI;
  Serial.print("data array repetitions = ");
  Serial.println(reps);

  // apply windowing function
  if (window == 1) //flat-top window
  {
    for (int i = 0; i < arraysize; i++)
    {
      float b = PI * i / (arraysize);
      a[i] = a[i] * ( 1 - (1.93 * cos(2 * b)) + (1.29 * cos(4 * b)) - (0.388 * cos(6 * b)) + (0.028 * cos(8 * b)));
      // Serial.println(a[i]);
    }
  }

  if (window == 2) //hann window
  {
    for (int i = 0; i < arraysize; i++)
    {
      float b = twopi * i / (arraysize);
      a[i] = a[i] * 0.5 * (1 - cos(b));
      //Serial.println(a[i]);
    }
  }

  if (window == 3) //hamming window
  {
    for (int i = 0; i < arraysize; i++)
    {
      float b = twopi * i / (arraysize);
      a[i] = a[i] * (0.54 - 0.46 * cos(b));
    //  Serial.println(a[i]);
    }
  }

  // DFT calculation

  for (float f = f0; f <= fn; f = f + stepsize)
  {

    sumi = sumr = 0.0;
    k = 0;
    for (int i = 0; i < (arraysize * reps); i++)
    {
      j = i - k;

      if (j >= arraysize) {
        k = k + arraysize;  //signal repeat index offset
        j = i - k;
      }
      ti = a[j] * (sin(twopi * f * i * interval * 0.001)); //0.001 => time in seconds
      tr = a[j] * (cos(twopi * f * i * interval * 0.001));
      sumi = sumi + ti;
      sumr = sumr + tr;
    }
    mag = sqrt(sumi * sumi + sumr * sumr) / (arraysize * reps);
    Serial.print(f);
    Serial.print("\t");
    Serial.println(mag);
  }

}
2 Likes

thx