QRS detector implementation in Java using the WFDB Swig library

I’ve decided that my first snippet will be a bit computer science based, well not a bit, a lot. Biomedical signal and image processing is a “realatively new breed” in computer science. New algorithms and researches are done on a daily basis, so that one day, computers could act as doctors themselves. Well we’re far from that, but the basics are here. The most basic thing we can do in biomedical signal processing is to write a QRS detector. In this post, I’ll explain some basics about signal processing and we’ll go through a quick ‘n’ dirty crash course of the LTST 1 database, reading files with the WFDB Java library and in the end, writing a QRS detector by HC Chen and SW Chen.

If anyone is more interested in the QRS detector that we’ll be implementing, here is the article (very light reading I might add lol - http://cinc.org/archives/2003/pdf/585.pdf). Let’s setup the environment. I recommend creating a 16GB virtual machine that has Ubuntu or some other equal distribution installed (VirtualBox of course). Why? Because we’ll be installing a lot of packages, and I really do mean A LOT! Download this file and run the following in the terminal:

sudo dpkg --set-selections < NamesceniPaketi-Pangolin.txt
sudo apt-get dselect-upgrad

Now that this is done (you’ll download about six gigs worth of packages); download WFDB from this site. Run the following in the terminal (yet again):

tar zxvf wfdb.tar.gz
cd wfdb.xx.yy.zz # xx.yy.zz is the package version
./configure
make
sudo make install

And we’re done, WFDB is installed. Now let’s install WFDB Swig (we’ll be using it for calling C functions from Java; you can also use it for Python or C#.) Download Swig from here. Now to the terminal again:

tar zxvf wfdb-swig.tar.gz
cd wfdb-swig.xx.yy.zz # where xx.yy.zz is the version of the current package
make
sudo make install

To test if the install was successfull, type “bxb” in the terminal and press enter. If the command was executed, everything is cool and ready. So now we need some data to test our QRS detector. The LTST database was made for this very purpose. Go to this page and download s20011.dat, s20011.hea and s20011.atr. The .dat file contains the ECG (electrocardiogram) data, hea contains a description of the header file and atr contains the reference annotations (done by doctors and 100% accurate); these are used for testing our detector.  Now open Netbeans and create a new project. Add the library WFDB.jar to the project and create a new main class. So now, let’s begin the fun stuff!

First we need to read the data file. (check the code bellow).

    public static int[] readData(int[] sig0, String record) {
        WFDB_SiginfoArray siarray = new WFDB_SiginfoArray(2);
        if(wfdb.isigopen(record, siarray.cast(), 2) < 2)
            System.exit(1);

        WFDB_SampleArray v = new WFDB_SampleArray(2);
        for(int i=0; i<siarray.cast().getNsamp(); i++) {
            if(wfdb.getvec(v.cast()) < 0)
                break;
            sig0[i] = v.getitem(0);
        }

        return sig0;
    }

I won’t go into details, but basically the parameter record is the name of the file (s20011 in our case), sig0 is the integer table that we’ll be writing to. Note that the program isn’t optimized, because we’ll be writing a lot of data in memory (so you’ll need at least 1.5 gigs of ram allocated). First we open the record and then we get the item and write it in the table. That’s it! Play around and output the data, so that you’ll see what it contains.

We need to send the data through some filters (just as the article says). First we need to put our data through a high pass filter, which is defined by the following formulas:

  y_1[n] = \frac{1}{M}  \sum^{M-1}_{m=0} x[n-m]

  y_2[n] = x[n - \frac{M+1}{2}]

  y[n] = y_2[n] - y_1[n]

Nothing so hard, but what is M? Well M is the window size and the article says that 5 or 7 is optimal. In our code, I selected 5. So the code isn’t very comlex, a for loop for calculating the elements and the sum. The only thing that we might miss is that you have to take the previous elements of the sum, if we start at position 0 (or 1-9), we have to take the previous ones (from the last indexes of the table). At the end we return a table.

public static float[] highPass(int[] sig0, int nsamp) {
        float[] highPass = new float[nsamp];
        float constant = (float) 1/M;

        for(int i=0; i<sig0.length; i++) {
            float y1 = 0;
            float y2 = 0;

            int y2_index = i-((M+1)/2);
            if(y2_index < 0) {
                y2_index = nsamp + y2_index;
            }
            y2 = sig0[y2_index];

            float y1_sum = 0;
            for(int j=i; j>i-M; j--) {
                int x_index = i - (i-j);
                if(x_index < 0) {
                    x_index = nsamp + x_index;
                }
                y1_sum += sig0[x_index];
            }

            y1 = constant * y1_sum;
            highPass[i] = y2 - y1;

        }        

        return highPass;
    }

We can output the first 2000 variables in a txt file (or as many as you like) and plot them using gnuplot.

echo 'plot "orig.txt"; pause 20' | gnuplot

Our output (if we output the first 2000 elements in a txt file and plot it ussing the upper code) now looks like this (awesome isn’t it?):

 

Now we have to pass our data through a filter yet again, this time a low pass one. It is implemented by the following formula (why the first 30 numbers, because the article says so – the size of the window is the most optimal one):

  y[0] = y_0^2 + y_1^2 + y_2^2 + \dots + y_{29}^2

y[n] = y_n^2 + y_{n+1}^2 + \dots + y_{n+29}^2

Well, what can be simpler? It’s only another for loop that iterates with a window size of 30 (we only have to watch that the index doesn’t go out of bounds).

    // Low pass filter; na n-to mesto zapiši kvadrat 30ih števil v oknu
    public static float[] lowPass(float[] sig0, int nsamp) {
        float[] lowPass = new float[nsamp];
        for(int i=0; i<sig0.length; i++) {
            float sum = 0;
            if(i+30 < sig0.length) {
                for(int j=i; j<i+30; j++) {
                    float current = sig0[j] * sig0[j];
                    sum += current;
                }
            }
            else if(i+30 >= sig0.length) {
                int over = i+30 - sig0.length;
                for(int j=i; j<sig0.length; j++) {
                    float current = sig0[j] * sig0[j];
                    sum += current;
                }
                for(int j=0; j<over; j++) {
                    float current = sig0[j] * sig0[j];
                    sum += current;
                }
            }

            lowPass[i] = sum;
        }

        return lowPass;

    }

Now we plot it again (we output the first 2000 elements in a txt file) and the data looks like this:

We can see, that the values are now normalized and above 0, and the beats can be detected easily.

All we need to do know is implement the “beat seeker”. The treshold formula is as follows:

Treshold = \alpha \times \gamma \times PEAK + (1 - \alpha) \times Treshold

0.01 \le \alpha \le 0.10

\gamma = 0.15, \gamma = 0.20

PEAK is the local maximum in the first 250 values we take. This is also our starting treshold. Alpha is “randomly” chosen between 0.01 and 0.10 and gamma can be 0.15 or 0.20 (also randomly chosen). We have a window of 250, we “walk” through all the values and we check if any of the values exceeds the treshold, if it does, we calculate a new treshold. The code:

    public static int[] QRS(float[] lowPass, int nsamp) {
        int[] QRS = new int[nsamp];

        double treshold = 0;

        for(int i=0; i<200; i++) {
            if(lowPass[i] > treshold) {
                treshold = lowPass[i];
            }
        }

        int frame = 250;

        for(int i=0; i<lowPass.length; i+=frame) {
            float max = 0;
            int index = 0;
            if(i + frame > lowPass.length) {
                index = lowPass.length;
            }
            else {
                index = i + frame;
            }
            for(int j=i; j<index; j++) {
                if(lowPass[j] > max) max = lowPass[j];
            }
            boolean added = false;
            for(int j=i; j<index; j++) {
                if(lowPass[j] > treshold && !added) {
                    QRS[j] = 1;
                    added = true;
                }
                else {
                    QRS[j] = 0;
                }
            }

            double gama = (Math.random() > 0.5) ? 0.15 : 0.20;
            double alpha = 0.01 + (Math.random() * ((0.1 - 0.01)));

            treshold = alpha * gama * max + (1 - alpha) * treshold;

        }

        return QRS;
    }

Now all we have to do, is write the values to a file. The following snippet does so (I won’t go into the details, because it’s all in the samples folder of the WFDB Swig library):

    public static void writeResults(int[] QRS, String record) {
        WFDB_Anninfo annInfo = new WFDB_Anninfo();
        WFDB_Annotation annot = new WFDB_Annotation();

        annInfo.setName("qrs");
        annInfo.setStat(wfdb.WFDB_WRITE);
        if(wfdb.annopen(record, annInfo, 1) < 0) {
            System.err.println("Problem pri pisanju v datoteko.");
            return;
        }

        annot.setSubtyp(0);
        annot.setChan(0);
        annot.setNum(0);
        annot.setAux(null);

        for(int i=0; i<QRS.length; i++) {
            if(QRS[i] != 0) {
                annot.setAnntyp(QRS[i]);
                annot.setTime(i);
                int error;
                if((error=wfdb.putann(0, annot)) < 0) {
                    System.err.println("Napaka pri pisanju v QRS datoteko.");
                    return;
                }
            }
        }

        wfdb.wfdbquit();

    }

And this is it! Now all we have to do is test our detector, we run it with the argument s20011, then it calculates for some time (about 60 seconds), and the output file s20011.qrs should be written in the same folder. Next we have to run the following command in the terminal to check our results:

bostjan@FaksLinux:~ bxb -r s20011 -a atr qrs
Beat-by-beat comparison results for record s20011
Reference annotator: atr
     Test annotator: qrs

               Algorithm
        n    v    f    q    o    x
   _______________________________
 N | 81889    0    0    0 17783    0
 V |    2    0    0    0    0    0
 F |    0    0    0    0    0    0
 Q |    0    0    0    0    0    0
 O |  181    0    0    0
 X |    0    0    0    0

           QRS sensitivity:  82.16% (81891/99674)
 QRS positive predictivity:  99.78% (81891/82072)
           VEB sensitivity:   0.00% (0/2)
 VEB positive predictivity:      -  (0/0)
   VEB false positive rate:  0.000% (0/82070)

  Beats missed in shutdown:   0.00% (0/99674)
      N missed in shutdown:   0.00% (0/99672)
      V missed in shutdown:   0.00% (0/2)
      F missed in shutdown:      -  (0/0)
       Total shutdown time:     0 seconds

So what does it say? Our QRS detector positively predicts 99.78% cases as heartbeats (QRS positive predictivity) and misses 0.23% of them. Pretty good. But we can optimize things even more. I’ve taken a default window of 250 in the QRS detection stage, and increment it by 250, we could increment it by 1 and make the detection even better.

Well, a zis is it! The whole source:

/* OBSS 1. seminarska naloga
 * QRS detector
 *
 * Author: Bostjan Cigan (http://zerocool.is-a-geek.net)
 *
 */

import wfdb.*;

public class QRS {

    public static final int M = 5;

    public static void main(String[] args) {

        int nsamp = getNrOfSamps(args[0]);
        int[] sig0 = new int[nsamp];
        sig0 = readData(sig0, args[0]);

        float[] highPass = highPass(sig0, nsamp);
        float[] lowPass = lowPass(highPass, nsamp);
        int[] QRS = QRS(lowPass, nsamp);

        writeResults(QRS, args[0]);

        System.out.println("Rezultati so zapisani v datoteki "+args[0]+".qrs ...");

    }

    public static float[] highPass(int[] sig0, int nsamp) {
        float[] highPass = new float[nsamp];
        float constant = (float) 1/M;

        for(int i=0; i<sig0.length; i++) {
            float y1 = 0;
            float y2 = 0;

            int y2_index = i-((M+1)/2);
            if(y2_index < 0) {
                y2_index = nsamp + y2_index;
            }
            y2 = sig0[y2_index];

            float y1_sum = 0;
            for(int j=i; j>i-M; j--) {
                int x_index = i - (i-j);
                if(x_index < 0) {
                    x_index = nsamp + x_index;
                }
                y1_sum += sig0[x_index];
            }

            y1 = constant * y1_sum;
            highPass[i] = y2 - y1;

        }        

        return highPass;
    }

    public static float[] lowPass(float[] sig0, int nsamp) {
        float[] lowPass = new float[nsamp];
        for(int i=0; i<sig0.length; i++) {
            float sum = 0;
            if(i+30 < sig0.length) {
                for(int j=i; j<i+30; j++) {
                    float current = sig0[j] * sig0[j];
                    sum += current;
                }
            }
            else if(i+30 >= sig0.length) {
                int over = i+30 - sig0.length;
                for(int j=i; j<sig0.length; j++) {
                    float current = sig0[j] * sig0[j];
                    sum += current;
                }
                for(int j=0; j<over; j++) {
                    float current = sig0[j] * sig0[j];
                    sum += current;
                }
            }

            lowPass[i] = sum;
        }

        return lowPass;

    }

    public static int[] QRS(float[] lowPass, int nsamp) {
        int[] QRS = new int[nsamp];

        double treshold = 0;

        for(int i=0; i<200; i++) {
            if(lowPass[i] > treshold) {
                treshold = lowPass[i];
            }
        }

        int frame = 250;

        for(int i=0; i<lowPass.length; i+=frame) {
            float max = 0;
            int index = 0;
            if(i + frame > lowPass.length) {
                index = lowPass.length;
            }
            else {
                index = i + frame;
            }
            for(int j=i; j<index; j++) {
                if(lowPass[j] > max) max = lowPass[j];
            }
            boolean added = false;
            for(int j=i; j<index; j++) {
                if(lowPass[j] > treshold && !added) {
                    QRS[j] = 1;
                    added = true;
                }
                else {
                    QRS[j] = 0;
                }
            }

            double gama = (Math.random() > 0.5) ? 0.15 : 0.20;
            double alpha = 0.01 + (Math.random() * ((0.1 - 0.01)));

            treshold = alpha * gama * max + (1 - alpha) * treshold;

        }

        return QRS;
    }

    public static int getNrOfSamps(String record) {
        WFDB_SiginfoArray siarray = new WFDB_SiginfoArray(2);
        if(wfdb.isigopen(record, siarray.cast(), 2) < 2)
            System.exit(1);

        return siarray.cast().getNsamp();
    }

    public static int[] readData(int[] sig0, String record) {
        WFDB_SiginfoArray siarray = new WFDB_SiginfoArray(2);
        if(wfdb.isigopen(record, siarray.cast(), 2) < 2)
            System.exit(1);

        WFDB_SampleArray v = new WFDB_SampleArray(2);
        for(int i=0; i<siarray.cast().getNsamp(); i++) {
            if(wfdb.getvec(v.cast()) < 0)
                break;
            sig0[i] = v.getitem(0);
        }

        return sig0;
    }

    public static void writeResults(int[] QRS, String record) {
        WFDB_Anninfo annInfo = new WFDB_Anninfo();
        WFDB_Annotation annot = new WFDB_Annotation();

        annInfo.setName("qrs");
        annInfo.setStat(wfdb.WFDB_WRITE);
        if(wfdb.annopen(record, annInfo, 1) < 0) {
            System.err.println("Problem pri pisanju v datoteko.");
            return;
        }

        annot.setSubtyp(0);
        annot.setChan(0);
        annot.setNum(0);
        annot.setAux(null);

        for(int i=0; i<QRS.length; i++) {
            if(QRS[i] != 0) {
                annot.setAnntyp(QRS[i]);
                annot.setTime(i);
                int error;
                if((error=wfdb.putann(0, annot)) < 0) {
                    System.err.println("Napaka pri pisanju v QRS datoteko.");
                    return;
                }
            }
        }

        wfdb.wfdbquit();

    }    

}
References:
  1. Franc Jager, Alessandro Taddei, George B. Moody, Michele Emdin, Gorazd Antolic, Roman Dorn, Ales Smrdel, Carlo Marchesi, and Roger G. Mark. Long-term ST database: a reference for the development and evaluation of automated ischaemia detectors and for the study of the dynamics of myocardial ischaemia. Medical & Biological Engineering & Computing 41(2):172-183 (2003) []

About this article