ekg.cs

心电信号处理库文件

using System;
using System.Collections.Generic;
using System.Text;

namespace DSProcessing
{
    /// <summary>
    /// ECG processing class, a part of DSProcessing library.
    /// author: Jan Sova
    /// mailto: [e-mail]
    ///
    /// -------------------------------------------------------------------------
    ///
    /// DSProcessing - C#/C++ library of signal processing, speech processing,
    ///              and communications classes and functions
    ///
    /// Copyright (C) 2007-2008
    ///
    /// This program is free software; you can redistribute it and/or modify
    /// it under the terms of the GNU General Public License as published by
    /// the Free Software Foundation; either version 2 of the License, or
    /// (at your option) any later version.
    ///
    /// This program is distributed in the hope that it will be useful,
    /// but WITHOUT ANY WARRANTY; without even the implied warranty of
    /// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    /// GNU General Public License for more details.
    ///
    /// You should have received a copy of the GNU General Public License
    /// along with this program; if not, write to the Free Software
    /// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
    ///
    ///-------------------------------------------------------------------------
    /// </summary>
    public class EKG
    {
        /// <summary>
        /// Implementation of Pan-Tompkinsonova algorithm for QRS detection
        /// </summary>
        /// <param name="data">EKG data</param>
        /// <returns>vektero s detekovanymi R-spickami</returns>
        /// <remarks>Jeste neumim filtfilt, tak delam konvoluci</remarks>
        public static double[] Rdetection(double[] data)
        {
            double[] celkMax = Tools.max(data);
            double[] returnData = new double[data.Length];
            double[] signal = new double[data.Length];
            Array.Copy(data, signal, data.Length);

            //odstraneni stejnosmerne slozky
            double mean = Tools.mean(signal);
            signal = Tools.minus(signal, mean);
            
            //pasmova filtrace mezi 4 a 30 Hz s nulovym fazovym posunem
            //koeficienty butter filteru
            //format long!!!! - vsimneme si, jak malo se meni koeficienty filtru a jak moc vysledek
            double[] a = { 1.00000000000000, -3.51001067675349, 4.65744210601898, -2.77721196640847, 0.63006969633312 };
            double[] b = { 0.02155836247819, 0, -0.04311672495637, 0, 0.02155836247819 };

            //double[] signalf = Tools.convolution(signal, butter);
            double[] signalf = Filter.filt(b, a, signal);
          
            //diference EKG signalu
            double[] sig_dif = Tools.diff(signalf);
            //umocneni signalu
            double[] sig_2  = Tools.pow(sig_dif, 2);

            //32-bodovy klouzavy prumer
            double[] b2 = new double[160];  
            double[] a2 = new double[160];
            for (int i = 0; i < b2.Length; i++)     
                b2[i] = 0.1563;
       
            a2[0] = 1;

            double[] sig_ma = Filter.filtfilt(b2, a2, sig_2);
            double prah = Tools.mean(sig_ma);

            //kpp=find(diff(sig_ma>prah)==1); 
            //kpn=find(diff(sig_ma>prah)==-1);
            int[] sig_ma_vetsi_prah = Tools.larger(sig_ma, prah);
            double[] sigVetsiDiff = Tools.diff(sig_ma_vetsi_prah);
            int[] kpp = Tools.find(sigVetsiDiff, 1);
            int[] kpn = Tools.find(sigVetsiDiff, -1);

            List<int> polohaMaxima = new List<int>();
            double[] maximum; //pro ulozeni jednoho kontretniho maxima
            
            //hledame pro oknech
            //vypocet maxima pro kazdy komplex
            for (int i = 0; i < kpn.Length; i++)
            {
                //priprava signalu, kde budu hledat maximum
                double[] okenko = new double[kpn[i] - kpp[i]];
                int n=0;
                for(int j = kpp[i]; j < kpn[i]; j++)
                {
                    okenko[n] = signal[j];
                    n++;
                }

                maximum = Tools.max(okenko);
                
                polohaMaxima.Add((int)maximum[1]+kpp[i]);
            }

            for (int i = 0; i < returnData.Length; i++)
                returnData[i] = 0;

            for (int i = 0; i < polohaMaxima.Count; i++)
                returnData[polohaMaxima[i]] = celkMax[0];
             
            return returnData;
        }
        /// <summary>
        /// Ze zadaneho vektrou vytahne casy
        /// </summary>
        /// <param name="data"></param>
        /// <returns></returns>
        public static double[] vratCasRSpicek(double[] data)
        {
            int counter = 0;
            int num = 0;

            for (int i = 0; i < data.Length; i++)
                if (data[i] > num)
                    counter++;

            double[] returnValue = new double[counter];
            counter = 0;

            for (int i = 0; i < data.Length; i++)
            {
                if (data[i] > num)
                {
                    returnValue[counter] = i;
                    counter++;
                }
            }

            return returnValue;
        }
        /// <summary>
        /// Hledani konce QRS komplexu
        /// </summary>
        /// <param name="x">vstupni signal</param>
        /// <param name="pma">pozice maxima</param>
        /// <param name="l">delka okna v ms</param>
        /// <param name="fs">vzorkovaci frekvence</param>
        /// <param name="pen">1 pro vzestupnou vlnu, -1 pro sestupnou</param>
        /// <param name="fac">prahovaci faktor</param>
        /// <returns>posledni bod vlny</returns>
        public static double[] konceQRS(double[] x, int[] pma, int l, int fs, int pen, int fac)
        {
            List<double> cum = new List<double>(); //sem se uklada vysledek 
            double dl=l*fs/1000;                                            // prevod ms _ vzorky

            for (int i = 0; i < pma.Length; i++)                                 // od prvni do posledni detekovane vlny
            {
                double[] xmu = Tools.fromToVector(x, pma[i], (int) (pma[i] + dl) );  //vyber aktualniho okna vlny
                double[] xmud = Tools.diff(xmu);                           //diferenciace
                double[] xmudXpen = Tools.mult(xmud, pen);                               //xmud*pen
                double[] maxima = Tools.max(xmudXpen);                                // nalezeni maxima derivace -- [m,n]
                xmud = Tools.fromToVector(xmud, (int) maxima[1], xmud.Length);    // omezeni okna
                double[] penXxmund = Tools.mult(xmud, pen);                  //pen*xmud
                int[] k = Tools.lesser(penXxmund, maxima[0] / fac);         // nalezeni vsech vzorku pod prahem
                int kpos = k[1];                                   // vzorek protinajici prah
                cum.Add(pma[i] + (int)maxima[1] + kpos -1);                                      // vypocet pozice
            }

            double[] returnValue = new double[cum.Count];

            for (int i = 0; i < returnValue.Length; i++)
            {
                returnValue[i] = cum[i];
            }

            return returnValue;
        }
    }
}