qrsfilt-inl.h

00001 #ifndef ECG_THIRD_EPLIMITED_QRSFILT_INL_H_
00002 #define ECG_THIRD_EPLIMITED_QRSFILT_INL_H_
00003 // Modified to OOP class sturcture
00004 // Author: Kasturi Rangan Raghavan
00005 /*****************************************************************************
00006 FILE:  qrsfilt.cpp
00007 AUTHOR:    Patrick S. Hamilton
00008 REVISED:    5/13/2002
00009   ___________________________________________________________________________
00010 
00011 qrsfilt.cpp filter functions to aid beat detecton in electrocardiograms.
00012 Copywrite (C) 2000 Patrick S. Hamilton
00013 
00014 This file is free software; you can redistribute it and/or modify it under
00015 the terms of the GNU Library General Public License as published by the Free
00016 Software Foundation; either version 2 of the License, or (at your option) any
00017 later version.
00018 
00019 This software is distributed in the hope that it will be useful, but WITHOUT ANY
00020 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
00021 PARTICULAR PURPOSE.  See the GNU Library General Public License for more
00022 details.
00023 
00024 You should have received a copy of the GNU Library General Public License along
00025 with this library; if not, write to the Free Software Foundation, Inc., 59
00026 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00027 
00028 You may contact the author by e-mail (pat@eplimited.edu) or postal mail
00029 (Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville,
00030 MA 02143 USA).  For updates to this software, please visit our website
00031 (http://www.eplimited.com).
00032   __________________________________________________________________________
00033 
00034     This file includes QRSFilt() and associated filtering files used for QRS
00035     detection.  Only QRSFilt() and deriv1() are called by the QRS detector
00036     other functions can be hidden.
00037     Revisions:
00038         5/13: Filter implementations have been modified to allow simplified
00039             modification for different sample rates.
00040 *******************************************************************************/
00041 #include <cstdlib>
00042 
00043 namespace qrsdet {
00044 
00045 /******************************************************************************
00046 * Syntax:
00047 *    int QRSFilter(int datum, int init) ;
00048 * Description:
00049 *    QRSFilter() takes samples of an ECG signal as input and returns a sample of
00050 *    a signal that is an estimate of the local energy in the QRS bandwidth.  In
00051 *    other words, the signal has a lump in it whenever a QRS complex, or QRS
00052 *    complex like artifact occurs.  The filters were originally designed for data
00053 *  sampled at 200 samples per second, but they work nearly as well at sample
00054 *    frequencies from 150 to 250 samples per second.
00055 *
00056 *    The filter buffers and static variables are reset if a value other than
00057 *    0 is passed to QRSFilter through init.
00058 *******************************************************************************/
00059 
00060 class RingBuffer {
00061  public:
00062   RingBuffer(size_t buf_size) {
00063     this->buf_size = buf_size;
00064     data = new int[buf_size];
00065       for (ptr = 0; ptr < buf_size; ptr++)
00066           data[ptr] = 0;
00067     ptr = 0;
00068   }
00069 
00070   ~RingBuffer() {
00071     delete[] data;
00072   }
00073 
00074   int* GetRange(int begin, int end, int* n) {
00075     *n = end - begin;
00076     int* mem = new int[*n];
00077     for (int p = begin, i = 0; p < end; p++, i++) {
00078       mem[i] = GetValue(p);
00079     }
00080     return mem;
00081   }
00082 
00083   void SetValue(int value, int offset = 0) {
00084     int set_ptr = ptr + offset;
00085     if (set_ptr < 0) 
00086       set_ptr += static_cast<int>(buf_size);
00087     set_ptr %= buf_size;
00088     data[set_ptr] = value;
00089   }
00090 
00091   int GetValue(int offset = 0) {
00092     int get_ptr = ptr + offset;
00093     if (get_ptr < 0)
00094       get_ptr += static_cast<int>(buf_size);
00095     get_ptr %= buf_size;
00096     return data[get_ptr];
00097   }
00098     
00099   void Increment() {
00100     if (++ptr == buf_size)
00101       ptr = 0;
00102   }
00103 
00104  protected:
00105   size_t buf_size;
00106   int* data;
00107   int ptr;
00108 };
00109 
00110 class DelayFilter : public RingBuffer {
00111  public:
00112   explicit DelayFilter(size_t buf_size) : RingBuffer(buf_size) {}
00113   int Process(int datum) {
00114     int output = GetValue();
00115     SetValue(datum);
00116     Increment();
00117     return output;
00118   }
00119 };
00120 
00121 /*************************************************************************
00122 *  lpfilt() implements the digital filter represented by the difference
00123 *  equation:
00124 *
00125 *     y[n] = 2*y[n-1] - y[n-2] + x[n] - 2*x[t-24 ms] + x[t-48 ms]
00126 *
00127 *    Note that the filter delay is (LPBUFFER_LGTH/2)-1
00128 *
00129 **************************************************************************/
00130 class LowPassFilter : public RingBuffer {
00131  public:
00132   LowPassFilter(size_t buf_size) : RingBuffer(buf_size) {
00133     y1 = 0;
00134     y2 = 0;
00135   }
00136 
00137   int Process(int datum) {
00138       int length = static_cast<int>(buf_size);
00139       long y0;
00140       int output, halfPtr;
00141       halfPtr = ptr - length / 2;
00142       if (halfPtr < 0)  // to x[n-6].
00143       halfPtr += length;
00144       y0 = (y1 << 1) - y2 + datum - (data[halfPtr] << 1) + data[ptr];
00145       y2 = y1;
00146       y1 = y0;
00147       output = y0 / ((length * length) / 4);
00148       data[ptr] = datum;
00149     Increment();
00150       return output;
00151   }
00152 
00153  private:
00154   long y1, y2;
00155 };
00156 
00157 
00158 /******************************************************************************
00159 *  hpfilt() implements the high pass filter represented by the following
00160 *  difference equation:
00161 *
00162 *    y[n] = y[n-1] + x[n] - x[n-128 ms]
00163 *    z[n] = x[n-64 ms] - y[n] ;
00164 *
00165 *  Filter delay is (buf_size-1)/2
00166 ******************************************************************************/
00167 class HighPassFilter : public RingBuffer {
00168  public:
00169   HighPassFilter(size_t buf_size) : RingBuffer(buf_size) {
00170         y = 0;
00171   }
00172 
00173   int Process(int datum) {
00174       int output, halfPtr;
00175       int length = static_cast<int>(buf_size);
00176     y += datum - data[ptr];
00177       halfPtr = ptr - (length / 2);
00178       if (halfPtr < 0)
00179           halfPtr += length;
00180       output = data[halfPtr] - (y / length);
00181       data[ptr] = datum;
00182       Increment();
00183       return output;
00184   }
00185 
00186  private:
00187   long y;
00188 };
00189 
00190 
00191 /*****************************************************************************
00192 *  deriv1 and deriv2 implement derivative approximations represented by
00193 *  the difference equation:
00194 *
00195 *    y[n] = x[n] - x[n - 10ms]
00196 *
00197 *  Filter delay is DERIV_LENGTH/2
00198 *****************************************************************************/
00199 class DerivFilter : public RingBuffer {
00200  public:
00201   DerivFilter(size_t buf_size) : RingBuffer(buf_size) {}
00202 
00203   int Process(int datum) {
00204     int y;
00205     y = datum - data[ptr];
00206       data[ptr] = datum;
00207     Increment();
00208       return y;
00209   }
00210 };
00211 
00212 
00213 /*****************************************************************************
00214 * mvwint() implements a moving window integrator.  Actually, mvwint() averages
00215 * the signal values over the last WINDOW_WIDTH samples.
00216 *****************************************************************************/
00217 class MovingWindowAverager : public RingBuffer {
00218  public:
00219   MovingWindowAverager(size_t buf_size) : RingBuffer(buf_size) {  
00220     sum = 0;
00221   }
00222 
00223   int Process(int datum) {
00224     int output;
00225     int length = static_cast<int>(buf_size);
00226       sum += datum;
00227       sum -= data[ptr];
00228       data[ptr] = datum;
00229     Increment();
00230       if((sum / length) > 32000)
00231           output = 32000 ;
00232       else
00233           output = sum / length;
00234       return output;
00235   }
00236  private:
00237   long sum;
00238 };
00239 
00240 struct QRSFilterParams {
00241   QRSFilterParams(double ts) {
00242     this->ts = ts;
00243   }
00244 
00245   size_t get_ms(size_t t) const { return static_cast<double>(t) / ts + 0.5; }
00246 
00247   double get_ts() const { return ts; }
00248 
00249   size_t get_lp_length() const { return get_ms(50); }
00250 
00251   size_t get_hp_length() const { return get_ms(125); }
00252 
00253   size_t get_d_length() const { return get_ms(10); }
00254 
00255   size_t get_window_length() const { return get_ms(80); }
00256 
00257   size_t get_preblank() const { return get_ms(200); }
00258 
00259   size_t get_filter_delay() const {
00260     return get_d_length() / 2.0 +
00261            get_lp_length() / 2.0 - 1.0 +
00262            (get_hp_length() - 1.0) / 2.0 +
00263            get_preblank();
00264   }
00265 
00266   size_t get_der_delay() const {
00267     return get_filter_delay() + get_window_length() + get_ms(100);
00268   }
00269 
00270  private:
00271   double ts;
00272 };
00273 
00274 class QRSFilter {
00275  public:
00276   QRSFilter(QRSFilterParams params) :
00277   lp_filter(params.get_lp_length()),
00278   hp_filter(params.get_hp_length()),
00279   d_filter(params.get_d_length()),
00280   mwa_filter(params.get_window_length()) {}
00281 
00282   int Process(int datum) {
00283     datum = lp_filter.Process(datum);
00284     datum = hp_filter.Process(datum);
00285     datum = d_filter.Process(datum);
00286     datum = abs(datum);
00287     datum = mwa_filter.Process(datum);
00288     return datum;
00289   }
00290 
00291  private:
00292   LowPassFilter lp_filter;
00293   HighPassFilter hp_filter;
00294   DerivFilter d_filter;
00295   MovingWindowAverager mwa_filter;
00296 };
00297 
00298 /**************************************************************
00299 * peak() takes a datum as input and returns a peak height
00300 * when the signal returns to half its peak height, or 
00301 **************************************************************/
00302 class PeakDetector {
00303  public:
00304   PeakDetector(const QRSFilterParams& params) {
00305     max = 0;
00306     timeSinceMax = 0;
00307     lastDatum = 0;
00308     ms_95_ = params.get_ms(95);
00309   }
00310 
00311   // what is Dly used for?
00312   int Process(int datum) {
00313       int pk = 0;
00314       if(timeSinceMax > 0)
00315           ++timeSinceMax ;
00316 
00317       if((datum > lastDatum) && (datum > max)) {
00318           max = datum;
00319           if(max > 2)
00320               timeSinceMax = 1;
00321         } else if(datum < (max >> 1)) {
00322       pk = max ;
00323       max = 0 ;
00324       timeSinceMax = 0 ;
00325           Dly = 0 ;
00326         } else if(timeSinceMax > ms_95_) {
00327           pk = max;
00328           max = 0;
00329           timeSinceMax = 0;
00330           Dly = 3;
00331         }
00332       lastDatum = datum ;
00333       return pk;
00334   }
00335  private:
00336   size_t ms_95_;
00337   int Dly;  // unused?
00338   int max;
00339   int timeSinceMax;
00340   int lastDatum;
00341 };
00342 
00343 /********************************************************************
00344 mean returns the mean of an array of integers.  It uses a slow
00345 sort algorithm, but these arrays are small, so it hardly matters.
00346 ********************************************************************/
00347 int mean(int *array, int n)
00348     {
00349     long sum = 0;
00350     for(int i = 0; i < n; ++i)
00351         sum += array[i];
00352     sum /= n;
00353     return sum;
00354     }
00355 
00356 /****************************************************************************
00357  thresh() calculates the detection threshold from the qrs mean and noise
00358  mean estimates.
00359 ****************************************************************************/
00360 int detection_thresh(int qmean, int nmean) {
00361   static const float TH = .3125;
00362   int thrsh, dmed;
00363   float temp;
00364   dmed = qmean - nmean;
00365   temp = dmed;
00366   temp *= TH;
00367   dmed = temp;
00368   thrsh = nmean + dmed;
00369   return thrsh;
00370 }
00371 
00372 /***********************************************************************
00373     BLSCheck() reviews data to see if a baseline shift has occurred.
00374     This is done by looking for both positive and negative slopes of
00375     roughly the same magnitude in a 220 ms window.
00376 ***********************************************************************/
00377 class BLSCheckBuffer : RingBuffer {
00378  public:
00379   BLSCheckBuffer(const QRSFilterParams& params)
00380     : RingBuffer(params.get_der_delay()) {
00381     ms150_ = params.get_ms(150);
00382     ms220_ = params.get_ms(220);
00383   }
00384 
00385   int Process(int datum) {
00386      data[ptr] = datum;
00387      Increment();
00388      return datum;
00389   }
00390 
00391   int BLSCheck(int *maxder) {
00392     int max = 0, min = 0;
00393     int maxt, mint;
00394     int saved_ptr = ptr;
00395       for(int t = 0; t < ms220_; t++) {
00396           
00397           int x = data[ptr];
00398           if(x > max) {
00399               maxt = t;
00400               max = x;
00401           } else if(x < min) {
00402               mint = t;
00403               min = x;
00404           }
00405         Increment();
00406       }
00407       ptr = saved_ptr;
00408       *maxder = max;
00409       min = -min;
00410     
00411       /* Possible beat if a maximum and minimum pair are found
00412           where the interval between them is less than 150 ms. */
00413          
00414       if((max > (min >> 3)) && (min > (max >> 3)) &&
00415           (abs(maxt - mint) < ms150_))
00416           return 0;
00417       else
00418           return 1;
00419   }
00420  private:
00421   size_t ms150_, ms220_;
00422 };
00423 
00424 
00425 }  // namespace qrsdet
00426 
00427 #endif  // ECG_THIRD_EPLIMITED_QRSFILT_INL_H_