Odd problem in ISR

Hello arduino community!
Wired question about a problem in my interrupt service routine.
I’m with Pulse Sensor www.pulsesensor.com, and I have a question about the weird way my firmware is working. The code can be found here
http://code.google.com/p/pulse-sensor/downloads/list
Called A_PulseSensor_06.
In it, I am using Timer1 to throw an interrupt every millisecond. In the isr, I am running a digital filter to help find the heartbeat and determine the heart rate. The filter needs to be seeded the first time through, or else it barfs out crazy noise. I’m using a Boolean called ‘first’ to test if it’s the first time through the filter, so it can be seeded, then the the ‘first’ flag is cleared.
My problem is in the if statement that tests ‘first’. I’m saying
If(first = true){

}
Which is wrong. Only one = is not a comparison. BUT it doesn’t work if I write it correctly.
If(first == true) and if(first) both cause the program to barf out illegible noise. Actually the noise looks the same as when I comment out the seeding routine completely.
I’ve tried on windows and Mac machines, also on UNO, decimillia, duemillanove, ardweenie, and my own compatible builds.always same result. The code works great with the one = , and more than 98% of my customers do not report this problem. Any insight into this?

I don’t download code from sites I don’t trust, and code.google is one of them. Post your code here, at least the declaration section and the ISR that is barfing.

If(first = true){
<seed filter>
}

May as well remove the if test, since it doesn’t do anything.

or else it barfs out crazy noise

both cause the program to barf out illegible noise

Perhaps a more scientific description?

Do you ever set “first” to true? Is it declared volatile? Where and how is it declared? Where is it changed? Hopefully an ISR is not so long you can’t post it here (and any relevant other parts). Preferably post the whole thing.

Here is the in full. The boolean variable ‘first’ is not being declared as a volatile variable. I have tried it as both volatile and not volatile in all the tests that I have did, with no noticeable difference. It looks like the last test was using a non-volatile ‘first’. A simple resistive sensor on Analog 0 will work as a stand in for the Pulse Sensor. The digital filter is stable, in that, once it starts working the input signal can go crazy, by disconnecting the sensor circuit, for example, and it will still recover when the sensor is re-attached. If the seed code section is commented out, or if it is written in correct syntax in all of my experience, a high frequency high amplitude noise is the result from which there can be no recovery. The noise is visible on pin 13, and if you look at the variable called ‘FSignal’.
This arduino code mates with a processing sketch called ‘P_PulseSensor_06’. I am not including it here because it is rather large. The question I have is not in the Processing sketch, which is a simple visualizer.
testing the timing using a digitalWrite at both ends of the ISR show ~150uS run time.

/*
>> Pulse Sensor Digital Filter <<
This code is the library prototype for Pulse Sensor by Yury Gitman and Joel Murphy
    www.pulsesensor.com 
    >>> Pulse Sensor purple wire goes to Analog Pin 0 <<<
Pulse Sensor sample aquisition and processing happens in the background via Timer 1 interrupt. 1mS sample rate.

*/
long Hxv[4]; // these arrays are used in the digital filter
long Hyv[4]; // H for highpass, L for lowpass
long Lxv[4];
long Lyv[4];

unsigned long readings; // used to help normalize the signal
unsigned long peakTime; // used to time the start of the heart pulse
unsigned long lastPeakTime = 0;// used to find the time between beats
volatile int Peak;     // used to locate the highest point in positive phase of heart beat waveform
int rate;              // used to help determine pulse rate
volatile int BPM;      // used to hold the pulse rate
int offset = 0;        // used to normalize the raw data
int sampleCounter;     // used to determine pulse timing
int beatCounter = 1;   // used to keep track of pulses
volatile int Signal;   // holds the incoming raw data
int NSignal;           // holds the normalized signal 
volatile int FSignal;  // holds result of the bandpass filter
volatile int HRV;      // holds the time between beats
volatile int Scale = 13;  // used to scale the result of the digital filter. range 12<>20 : high<>low amplification
volatile int Fade = 0;

boolean first = true; // reminds us to seed the filter on the first go
volatile boolean Pulse = false;  // becomes true when there is a heart pulse
volatile boolean B = false;     // becomes true when there is a heart pulse
volatile boolean QS = false;      // becomes true when pulse rate is determined. every 20 pulses

int pulsePin = 0;  // pulse sensor purple wire connected to analog pin 0


void setup(){
pinMode(13,OUTPUT);    // pin 13 will blink to your heartbeat!
Serial.begin(115200); // we agree to talk fast!
// this next bit will wind up in the library. it initializes Timer1 to throw an interrupt every 1mS.
TCCR1A = 0x00; // DISABLE OUTPUTS AND BREAK PWM ON DIGITAL PINS 9 & 10
TCCR1B = 0x11; // GO INTO 'PHASE AND FREQUENCY CORRECT' MODE, NO PRESCALER
TCCR1C = 0x00; // DON'T FORCE COMPARE
TIMSK1 = 0x01; // ENABLE OVERFLOW INTERRUPT (TOIE1)
ICR1 = 8000;   // TRIGGER TIMER INTERRUPT EVERY 1mS  
sei();         // MAKE SURE GLOBAL INTERRUPTS ARE ENABLED

}



void loop(){
  Serial.print("S");          // S tells processing that the following string is sensor data
  Serial.println(Signal);     // Signal holds the latest raw Pulse Sensor signal
  if (B == true){             //  B is true when arduino finds the heart beat
    Serial.print("B");        // 'B' tells Processing the following string is HRV data (time between beats in mS)
    Serial.println(HRV);      //  HRV holds the time between this pulse and the last pulse in mS
    B = false;                // reseting the B for next time
  }
  if (QS == true){            //  QS is true when arduino derives the heart rate by averaging HRV over 10 beats
    Serial.print("Q");        //  'Q' tells Processing that the following string is heart rate data
    Serial.println(BPM);      //  BPM holds the heart rate in beats per minute
    QS = false;               //  reset the QS for next time
  }
  Fade -= 15;                    // an LED on pin 11 will fade out this much if it has been turned on
  Fade = constrain(Fade,0,255);  // don't roll over the fade
  analogWrite(11,Fade);          
  
delay(20);                    //  take a break

}

// THIS IS THE TIMER 1 INTERRUPT SERVICE ROUTINE. IT WILL BE PUT INTO THE LIBRARY
ISR(TIMER1_OVF_vect){ // triggered every time Timer 1 overflows
// Timer 1 makes sure that we take a reading every milisecond
Signal = analogRead(pulsePin);

// First normailize the waveform around 0
readings += Signal; // take a running total
sampleCounter++;     // we do this every milisecond. this timer is used as a clock
if ((sampleCounter %300) == 0){   // adjust as needed
  offset = readings / 300;        // average the running total
  readings = 0;                   // reset running total
}
NSignal = Signal - offset;        // normalizing here

// IF IT'S THE FIRST TIME THROUGH THE SKETCH, SEED THE FILTER WITH CURRENT DATA
if (first = true){
  for (int i=0; i<4; i++){
    Lxv[i] = Lyv[i] = NSignal <<10;  // seed the lowpass filter
    Hxv[i] = Hyv[i] = NSignal <<10;  // seed the highpass filter
  }
first = false;      // only seed once please
}
// THIS IS THE BANDPAS FILTER. GENERATED AT www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html
//  BUTTERWORTH LOWPASS ORDER = 3; SAMPLERATE = 1mS; CORNER = 5Hz
    Lxv[0] = Lxv[1]; Lxv[1] = Lxv[2]; Lxv[2] = Lxv[3];
    Lxv[3] = NSignal<<10;    // insert the normalized data into the lowpass filter
    Lyv[0] = Lyv[1]; Lyv[1] = Lyv[2]; Lyv[2] = Lyv[3];
    Lyv[3] = (Lxv[0] + Lxv[3]) + 3 * (Lxv[1] + Lxv[2])
          + (3846 * Lyv[0]) + (-11781 * Lyv[1]) + (12031 * Lyv[2]);
//  Butterworth; Highpass; Order = 3; Sample Rate = 1mS; Corner = .8Hz
    Hxv[0] = Hxv[1]; Hxv[1] = Hxv[2]; Hxv[2] = Hxv[3];
    Hxv[3] = Lyv[3] / 4116; // insert lowpass result into highpass filter
    Hyv[0] = Hyv[1]; Hyv[1] = Hyv[2]; Hyv[2] = Hyv[3];
    Hyv[3] = (Hxv[3]-Hxv[0]) + 3 * (Hxv[1] - Hxv[2])
          + (8110 * Hyv[0]) + (-12206 * Hyv[1]) + (12031 * Hyv[2]);
FSignal = Hyv[3] >> Scale;  // result of highpass shift-scaled

//PLAY AROUND WITH THE SHIFT VALUE TO SCALE THE OUTPUT ~12 <> ~20 = High <> Low Amplification.

if (FSignal >= Peak && Pulse == false){  // heart beat causes ADC readings to surge down in value.  
  Peak = FSignal;                        // finding the moment when the downward pulse starts
  peakTime = sampleCounter;              // recodrd the time to derive HRV. 
}
//  NOW IT'S TIME TO LOOK FOR THE HEART BEAT
if ((sampleCounter %20) == 0){// only look for the beat every 20mS. This clears out alot of high frequency noise.
  if (FSignal < 0 && Pulse == false){  // signal surges down in value every time there is a pulse
     Pulse = true;                     // Pulse will stay true as long as pulse signal < 0
     digitalWrite(13,HIGH);            // pin 13 will stay high as long as pulse signal < 0  
     Fade = 255;                       // set the fade value to highest for fading LED on pin 11 (optional)   
     HRV = peakTime - lastPeakTime;    // measure time between beats
     lastPeakTime = peakTime;          // keep track of time for next pulse
     B = true;                         // set the Quantified Self flag when HRV gets updated. NOT cleared inside this ISR     
     rate += HRV;                      // add to the running total of HRV used to determine heart rate
     beatCounter++;                     // beatCounter times when to calculate bpm by averaging the beat time values
     if (beatCounter == 10){            // derive heart rate every 10 beats. adjust as needed
       rate /= beatCounter;             // averaging time between beats
       BPM = 60000/rate;                // how many beats can fit into a minute?
       beatCounter = 0;                 // reset counter
       rate = 0;                        // reset running total
       QS = true;                       // set Beat flag when BPM gets updated. NOT cleared inside this ISR
     }
  }
  if (FSignal > 0 && Pulse == true){    // when the values are going up, it's the time between beats
    digitalWrite(13,LOW);               // so turn off the pin 13 LED
    Pulse = false;                      // reset these variables so we can do it again!
    Peak = 0;                           // 
  }
}

}// end isr

Before we even look at "first" let's look at the timers. How does this give you an interrupt ever 1 mS?

  TCCR1A = 0x00; // DISABLE OUTPUTS AND BREAK PWM ON DIGITAL PINS 9 & 10
  TCCR1B = 0x11; // GO INTO 'PHASE AND FREQUENCY CORRECT' MODE, NO PRESCALER
  TCCR1C = 0x00; // DON'T FORCE COMPARE
  TIMSK1 = 0x01; // ENABLE OVERFLOW INTERRUPT (TOIE1)
  ICR1 = 8000;   // TRIGGER TIMER INTERRUPT EVERY 1mS

For a start, at 62.5 nS per clock cycle, 8000 cycles is 500 uS, not 1000 uS. Second you have the interrupt on overflow, so the counter is rather irrelevant. It overflows every 65536 * 62.5 nS which is every 4.096 mS.

Having said that, I suspect there is something wrong with your filter, which I don't totally understand.

ok, let's look at the timer. notice that i am putting Timer1 in 'phase and frequency correct' mode. the timer counts up for 500uS, then counts down for 500uS, then trips the interrupt, which happens every 1000uS. for more information RTDS. the filter is complicated, there is a readme.TXT file located on the google code page in original post that tells the story about how i built it. i am hoping that someone would compile and upload the code with a simple analog circuit hanging off of pin 0 and help me trouble-shoot why it only works with one '=' sign in the if statement that measures 'first'. Please try it by initializing 'first' as a volatile boolean, if you like, I have found that the issue i'm asking about doesn't care about that. Someone on our forum had a problem when she loaded the code on her duemilanove. high frequency noise coming out of pin 13. I told her about the issue, and she corrected the if statement in question, and then it worked. http://pulsesensor.proboards.com/index.cgi?board=allthebitsarehere&action=display&thread=37 As I mentioned, there is a small percentage of people who have this problem with the code. for most people it works with only one '='.

OK, fair enough. Testing seems to confirm that.

testing the timing using a digitalWrite at both ends of the ISR show ~150uS run time.

250 uS more like.

To save me a lot of guessing can you specify what sort of signal on A0 would produce the results you are testing for?

… help me trouble-shoot why it only works with one ‘=’ sign in the if statement that measures ‘first’.

Let’s get rid of this “one ‘=’ sign” talk. By putting “first = true” you effectively have not this:

 if (first = true){
    for (int i=0; i<4; i++){
      Lxv[i] = Lyv[i] = NSignal <<10;  // seed the lowpass filter
      Hxv[i] = Hyv[i] = NSignal <<10;  // seed the highpass filter
    }
    first = false;      // only seed once please
  }

But rather:

   for (int i=0; i<4; i++){
      Lxv[i] = Lyv[i] = NSignal <<10;  // seed the lowpass filter
      Hxv[i] = Hyv[i] = NSignal <<10;  // seed the highpass filter
    }

So your question really is: “why does it only work if I seed the filter every time”?

Can you confirm what you are expecting to see here? "High frequency noise" could mean anything. Testing pin 13 with waving my hand over A0 results in it outputting a series of pulses around 20 to 80 mS long with a frequency of roughly 6 to 20 Hz. Hardly "high frequency noise".

I don't really know what your gadget is trying to achieve. That is, what output you are expecting for what input. If you clarify all that we might be able to move on.

So your question really is: "why does it only work if I seed the filter every time"?

well, that's weird enough. the program does indeed work well when the filter is seeded every time. the question becomes one more of DSP and digital filter design.

Can you confirm what you are expecting to see here?

The code is designed to interface with a Pulse Sensor, an optical heart-rate monitor. Should see a saw-tooth waveform approx 50mV to 250mV amplitude (sometimes larger) at ~1 to ~3 Hz. It also outputs a blink on pin 13 LED that is times with your heartbeat. When the seeding of the filter is taken out all together, and, i now see, when it is run only once, there is a high frequency flickering of pin 13 LED.

JoelSensor: the question becomes one more of DSP and digital filter design.

Yes, exactly. It's nothing to do with first-time flags and ISRs. And I personally don't know enough about filter design to help on this point.

I’m with Nick. It looks like your digital filter simply isn’t working. By having the first-time code execute all the time, you essentially end up passing an unfiltered value on to the next stage of code…

I don’t understand digital filters either, but a couple of things stood out:

  1. A good filter should stabilize relatively quickly even when it starts out with 0 values. Seeding could make that happen faster, but it should be NECESSARY.
  2. Your scale factors are inconsistent. You scale input<<10, but output>>Scale (and Scale==13)
  3. It’s be nice to see the original coeffecients in their FP form, at least in the comments. It’d be even nicer to use a macro to scale those coefficents to the scaled integers you’re using.
  4. You have
    Lyv[3] = (Lxv[0] + Lxv[3]) + 3 * (Lxv[1] + Lxv[2]) + (3846 * Lyv[0]) + (-11781 * Lyv[1]) + (12031 * Lyv[2]);
    It’s bothering me to see apparently unscaled constants (3) and scaled constants (3846, -11781) in the same equation. But I’m not quite seeing how the math should be, since the samples were scaled before they were put into the array as well…

Have you debugged the filter on a desktop, using created waveform data? In both the floating point and scaled integer forms? I created a 3rd order highpass filter with 0.5Hz corner freq using the website you reference, and called the code in its original form (floating point) with a constant value, and got approximately zero out. Which is about what you’d expect, since my input freq was 0. To get non-zero output I should send in a waveform whose frequency is significantly above 0.5Hz, right? (and indeed if I enter a 500Hz square wave, I get about the same thing out. (shifted from 0/.7 to +.35/-.35, which is sort of interesting…)

I didn’t see anything on the expected range of the input signal. I guess it shouldn’t matter unless I exceed the number system?

I think it bothers me that the sample history covers a much shorter timeframe than one cycle of the expected waveform, but I don’t know enough about DSP to say that that’s not OK. Samples seems to equal poles for whatever algorithm is in use, and a 100+ pole filter doesn’t seem practical…

studying... Why did you choose a low-pass filter followed by a high-pass filter, instead of just doing a bandpass filter to start with? (It doesn't look like it would be less calculation, just ... simpler?)

    Hxv[3] = Lyv[3] / 4116; // insert lowpass result into highpass filter

Isn't Lyv already scaled? I'm not understanding your fixed point scaled arithmetic at all (4116 !!??) (and I did read the .txt file, which mostly only says "I made up some numbers and it works well", except apparently it doesn't!)

Jeez westfw you must be a computer nerd, starting a bulleted list at 0 :)


Rob

you must be a computer nerd

Someone had a doubt?

starting a bulleted list at 0

I didn't, actually :-) I just decided that the item labeled 0 should come before the item labeled 1 :-) (now, to be a true old-fashioned computer nerd, I should have had 5, 10, 20, 30... ?)

The curves from http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html certainly look nicer when the sample frequency is less than 1000 times the input waveform frequency. Say, 50Hz...

I didn't, actually :-) I just decided that the item labeled 0 should come before the item labeled 1 :-)

Yep, always use a nasty patch rather then a rewrite, that's a true hacker tradition. ;)

westfw: (now, to be a true old-fashioned computer nerd, I should have had 5, 10, 20, 30... ?)

You know, I'd forgotten about doing that.

I think the problem is in your scaled math. If you scale a multiplicative constant, you can only use it ONCE, after which you have to start re-normalizing your results.

For instance, using a scale factor of 1000 (base 10 to make things easy), multiplying by 0.234 would mean multiplying by 234 instead. So far so good. But if you take that result and multiply by the scaled constant again, you've multiplied by 54756 (234*234), which is not at all the same as multiplying by a scaled (.234*.234*1000 != 55)

I didn't see the renormalization that you'd have to do. It might have been buried in the mysterious shifts and divisions, but 1) You should make such things obvious by letting the compiler do constant math

#define SCALEDFACTOR(a) ((long)(a*2.0^12))
#define SCALEDMULT(a,x) ((SCALEDFACTOR(a)*x)>>12

 .... + (SCALEDMULT(0.904347,  Lyv[0])

Compilers these days are very good at doing math on constants; don't make your code more obscure by doing it for them! 2) I think this would explain the results you're seeing. Multiply things by a few extra thousands a couple time in a row and you'll overflow your longs, as well as having intermediate incorrect results. 3) It'd be nice to see the actual equations, scaling factors, and etc CLEARLY DOCUMENTED. You're trying to do:

 yv[3] =   (xv[0] + xv[3]) + 3 * (xv[1] + xv[2])
                     + (  0.9390989403 * yv[0]) + ( -2.8762997235 * yv[1])
                     + (  2.9371707284 * yv[2]);

Using integers scaled by a factor of 4096 (12 bit shift), right?