How to debug Error in Pan Tompkins?

I am getting the following error when the AF signal is passed in Pan Tompkins algorithm.

Error using filtfilt>getCoeffsAndInitialConditions (line 181)
Data length must be larger than 18 samples.

Error in filtfilt (line 96)
[b,a,zi,nfact,L] = getCoeffsAndInitialConditions(b,a,Npts);

Error in pan_tompkin2 (line 194)
ecg_h = filtfilt(a,b,ecg);

I have downloaded AF Signals (MIT-BIH Arrhythmia Database as .mat) from physiobank-atm from physionet.org. Could you help in debugging it?

[Walter Roberson]

Reference: MIT-BIH Arrhythmia Database at https://physionet.org/physiobank/database/mitdb/

Instructions on downloading: http://www.mathworks.com/matlabcentral/answers/10331-how-to-open-mit-bih-arrhythmia-database-file-in-matlab

[Walter Roberson]

For the Arrhythmia database, which "record" are you downloading? Are you downloading MLII or V5 ? On the MATLAB side, how are you loading the data into MATLAB, and how are you passing the data into pan_tomkin2?

My suspicion is that you might be trying to pass the name of the .mat instead of the content of the .mat

[Shraddha Joshi]

I have attached a signal sample which I have downloaded along with Pan tompkins code I am using. Signal is downloaded as

In the Pan Tompkins code attached below, in the lines 118 and 153, the command i have used to load signals is,

ecg=load('100.mat')

instead of

ecg=load('ECG_sample_noisy.mat')

Answer 1

[Walter Roberson]

matstruct = load('100m.mat');
ecg = matstruct.val;

That is, you need to use the .val variable within the .mat. When you load() a .mat like that, returning a value, the result is a struct in which the fields are the names of the loaded variables.

[Shraddha Joshi]

Sir, when I tried including the above instructions you said, in the lines 118 and 119 of the program attached below, I got following error:

Error using pan_tompkin2 (line 122) ecg must be a row or column vector

I am not getting where exactly to provide the above instructions in the original Pan Tompkin Algorithm code. Can u please help me?

[Walter Roberson]

matstruct = load('100m.mat');
ecg1 = matstruct.val(1,:);
ecg2 = matstruct.val(2,:);

Now ecg1 is the MLII signal and ecg2 is the V5 signal; process whichever one is appropriate for your purpose.

[Shraddha Joshi]

Thank You Sir. It did help in getting the output. One more help i needed, there are a few MITBIH Afib sample signals that I have downloaded which I have enclosed in the attached folder below. How can I pass all these sample signals one after the other through Pan Tompkins Algorithm and check and store their rr interval values?

[Shraddha Joshi]

Sir could u please help regarding above query at the earliset? How can the code correction be done to read all the mitbih afib signals one after the other through Pan tompkins algorithm and for each signal the rr interval calculated should be stored in a register or array?

[Walter Roberson]

My health is not good; assistance from me is irregular.

[Walter Roberson]

http://matlab.wikia.com/wiki/FAQ#How_can_I_process_a_sequence_of_files.3F

fs = 360;   %according to the .info files

dinfo = dir('*.mat');
filenames = {dinfo.name};
nfiles = length(filenames);
results = cell(nfiles, 1);
for K = 1 : nfiles
   thisfile = filenames{K};
   matstruct = load(thisfile);
   ecgs = matstruct.val;
   for rownum = 1 : size(ecgs,1)
      ecg = ecgs(rownum, :);
      results{K,rownum} = pan_tompkin2(ecg, fs, ... whatever);
    end
  end

Even better would be to go through the header files and read out the frequency, but I did not bother; all the ones I looked at were 360 Hz.

[Shraddha Joshi]

Sorry for the inconvenience Sir. And yes the frequency is 360Hz for all signals. I tried with the above code but the error I am getting is

Undefined function or variable "ecgs".

Error in pan_tompkin2 (line 130) for rownum=1:size(ecgs,1);

What is this ecgs and how can I define it?

[Walter Roberson]

It is defined on the line above where it is used.

[Shraddha Joshi]

Thank you so much sir!!

[Shraddha Joshi]

Is it possible to calculate sensitivity and specificity for the above signals? If so can u please tell me how to do so?

[Walter Roberson]

Sensitivity and specificity of what? You get back the R-wave indices and not much else. Do you have information about where the R-waves "really" are that you can compare against?

[Shraddha Joshi]

According to the code u have sent before which is as below, we are storing the results in results{K,rownum}.

fs = 360;   %according to the .info files
dinfo = dir('*.mat');
filenames = {dinfo.name};
nfiles = length(filenames);
results = cell(nfiles, 1);
for K = 1 : nfiles
   thisfile = filenames{K};
   matstruct = load(thisfile);
   ecgs = matstruct.val;
   for rownum = 1 : size(ecgs,1)
      ecg = ecgs(rownum, :);
      results{K,rownum} = pan_tompkin2(ecg, fs, ... whatever);
    end
  end

So as this results{K,rownum} gives information about RR intervals, can't we calculate Sensitivity and Specificity of the Afib signals we are passing? Or is there any other way from which we can calculate sensitivity and specificity of the signals?

[Walter Roberson]

You cannot calculate those without some information about what the accurate classification is. You could, for example, synthesize some data and put it through and then compare the results to the known values.

[Shraddha Joshi]

Suppose if I check these rr intervals stored in results{k,rownum} with some standard rr interval range for different kind of arrhythmia based on decision rules and classify them as normal or abnormal ecg signals and again in abnormal as Afib or other arrhythmia; will i be able to calculate sensitivity and specificity then? If so, how exactly can i do it through code?

[Walter Roberson]

No, you need some information known to be true, not information that just had a different approximate classification method used. That information an be known to be true because you synthesized the information, or it can be known to be true by expert analysis of the signals.

You have the additional difficulty that your signal bounds are probably slightly fuzzy rather than sample-for-sample exact. You are sampling at 360 Hz, so if you have an event that takes less than 1/720 of a second to start, it could might legitimately get associated with either the sample #K or sample #(K-1). An expert looking at the plot might point to one particular sample one time as being the boundary, and the same expert might point to an adjacent sample the next time as being the same boundary. Especially in the presence of noise, it is often impossible to say that a particular sample is exactly the right boundary. Right at the calculated boundary it probably does not matter much whether the sample is inside or outside the boundary -- but (for example) 10 samples inside the boundary if you haven't yet agreed the boundary is there, then You Have A Problem that needs to be measured.

Some classification methods account for this by generating a confidence in their calculation, and then the statistical measures can be weighted by the confidence. But you do not appear to have that situation here.

[Shraddha Joshi]

So there is absolutely no way of determining sensitivity and specificity for the signals above program? At least by considering no.of heart beats? Or else could you please suggest any other way which would help me to calculate them? I really needed this calculation in order for further analysis with my project studies. Hoping for the help.

[Shraddha Joshi]

Or is there a way to check if the above signals passed are atrial fibrillation or no, based on heart rate calculation from rr intervals obtained ?

[Walter Roberson]

Without some objective truth or expert analysis, your system might as well just be playing NumberWang

[Shraddha Joshi]

So are you trying to say that there is no way to calculate heart rate from rr intervals found?

[Walter Roberson]

I do not know anything about calculating heart rate from ecg signals. I do not know what rr intervals are. I was responding about sensitivity and specificity: to calculate those, you need information known by construction or expert opinion to be true.

[Shraddha Joshi]

Ok thank you!!