package org.physionet.challenge2011;
import java.io.IOException;
import java.io.InputStream;
import java.io.ObjectInputStream;
/*
* Class: ChallengeEntry
* Algorithm to compute classification of ECG records.
* Returned value = 1 (acceptable) or 0 (unacceptable)
*
* @ authors: Philip Langley, Philip.Langley@ncl.ac.uk
* David Duncan, david.duncan@nuth.nhs.uk
* Note: This is a partial implementation of the algorithm developed
* for part 1. Flat line, Saturation, Baseline drift and Low Amplitude
* checks are implemented. High amplitude and Steep slope are not implemented.
*/
public class ChallengeEntry {
public static final String DEBUGTAG = ChallengeEntry.class.toString();
final static int FS = 500; // Sampling Frequency
final static int DW = FS * 2; // 2s (skip first DW samples in BL analysis)
final static int DW_Sat = Math.round(200 * FS / 1000);
// ======= THRESHOLDS ============
final static int SAT_THR = 400;
final static int D1_Thr = 50;
final static int dE_min_Thr = 25; // units [LSB] max elongation of ECG from BL
static int D1_Thr_Count = 0;
final static int BLW_max_Thr = 500; // units [LSB] Baseline amplitude swing
// (max)
final static int CH = 12;
final static int MAX_RT = 220; // Max expected beats in minutes
final static int WIN = FS * 10;
static int[][] Thr_Count_Array = new int[12][2];
static double[][] filtered_ecg_array = new double[WIN][CH];
short[] data = new short[WIN * CH];
synchronized public int get_result(InputStream iFile,
final ECG_MetaData m_MetaData) throws IOException {
ObjectInputStream in = new ObjectInputStream(iFile);
try {
data = (short[]) in.readObject();
} catch (ClassNotFoundException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
// preset result (1=acceptable)
int result = 1;
int bUnacc = 0;
// Load data into 2D array
int COLS = CH;
int ROWS = WIN;
short[][] ecg_array = loadArrayData(data, COLS, ROWS);
// verify flat line
FlatLine: for (int c = 0; c < ecg_array[1].length; c++) {
short[] col_array = new short[ROWS];
for (int r = 0; r < ecg_array.length; r++) {
col_array[r] = ecg_array[r][c];
}
// Detect Flat Line
if (detectFlat(col_array, c)) {
bUnacc = 1;
break FlatLine;
}
}
// test for saturation
if (bUnacc == 0) {
Saturation: for (int c = 0; c < ecg_array[1].length; c++) {
short[] col_array = new short[ROWS];
for (int r = 0; r < ecg_array.length; r++) {
col_array[r] = ecg_array[r][c];
}
// Detect Saturation
if (detectSaturation(col_array, c)) {
bUnacc = 1;
break Saturation;
}
}
}
if (bUnacc == 0) {
// determine BL for all leads
for (int c = 0; c < ecg_array[1].length; c++) {
double[] col_array = new double[ROWS];
for (int r = 0; r < ecg_array.length; r++) {
col_array[r] = ecg_array[r][c];
}
// filter
setFilteredECGData(c, filter(col_array), filtered_ecg_array);
}
// ECG elongation from BL (dE)
double[][] dE = subtractArray(ecg_array, filtered_ecg_array);
// LPF for PM rejection
//double[][] ecgLPF = LPF_ecg(dE);
// collect baseline wander (BLW) and elongation (dE) for all leads
double[][] v = new double[12][3];
for (int c = 0; c < filtered_ecg_array[1].length; c++) {
double[] blwander_array = new double[filtered_ecg_array.length];
double[] v2_array = new double[filtered_ecg_array.length];
for (int r = 0; r < filtered_ecg_array.length; r++) {
blwander_array[r] = filtered_ecg_array[r][c];
v2_array[r] = dE[r][c];
}
v[c][0] = getBLWander(blwander_array); // BLW units [LSB]
v[c][1] = getAbsMax(v2_array); // dE max, units [LSB]
// v(ld, 3) = max(abs(yLPF(:,ld))); % dE-LPF max, units [LSB]
}
double min_val = v[0][0];
for (int r = 0; r < 12; r++) {
if (v[r][0] < min_val) {
min_val = v[r][0];
}
}
if ((min_val < BLW_max_Thr)) {
// stable baseline
// rejection criterion:: 1) amplitude lower than dE min threshold in ANY lead
int dE_min_Thr2 = dE_min_Thr + 10; // units [uV]
double min_value = v[0][1];
int cnt = 0;
for (int r = 0; r < 12; r++) {
if (v[r][1] < min_value) {
min_value = v[r][1];
}
if (v[r][1] < dE_min_Thr2) {
cnt++;
}
}
if( min_value<dE_min_Thr ){
// mark unacceptable
bUnacc = 1;
}else if (min_value<dE_min_Thr2){
if( cnt>=3 ){
// mark unacceptable
bUnacc = 1;
}
}
}else{
//unstable baseline
bUnacc = 1;
}
}
if (bUnacc == 1){
// set result to unacceptable
result = 0;
}
// System.out.println("Result = " + result + " Flat = " + flatFlag
// + ", Saturation = " + saturationFlag + ", SBL = " + isStableBL+", MinThresh1 = "+minThresh1
// + ", MinThresh2 = "+minThresh2);
// clean-up
iFile.close();
return result;
}
/*
* Load 1D array of data into a 2D array
*/
private short[][] loadArrayData(short[] data, int cols, int rows)
throws NumberFormatException, IOException {
short[][] ecg_array = new short[rows][cols];
int linear_index = 0;
for (int r = 0; r < rows; r++) {
for (int c = 0; c < cols; c++) {
ecg_array[r][c] = data[linear_index];
linear_index++;
}
}
return ecg_array;
}
private double getAbsMax(double[] dE){
double max_val = dE[0];
for (int r = 0; r < dE.length; r++) {
if (Math.abs(dE[r]) > max_val) {
max_val = Math.abs(dE[r]);
}
}
return max_val;
}
private double getBLWander(double[] ecg) {
int Le = ecg.length;
double max_val = 0;
for (int i = DW; i < Le; i++) {
if (ecg[i] > max_val) {
max_val = ecg[i];
}
}
double min_val = max_val;
for (int i = DW; i < Le; i++) {
if (ecg[i] < min_val) {
min_val = ecg[i];
}
}
double BLW = max_val - min_val;
return BLW;
}
private double[][] setFilteredECGData(int col, double[] filtered1DArray,
double[][] filteredecg) {
for (int r = 0; r < filtered1DArray.length; r++) {
filteredecg[r][col] = filtered1DArray[r];
}
return null;
}
private double[][] subtractArray(short[][] a, double[][] b) {
double[][] dE = new double[a.length][a[1].length];
for (int c = 0; c < a[1].length; c++) {
for (int r = 0; r < a.length; r++) {
dE[r][c] = a[r][c] - b[r][c];
}
}
return dE;
}
private double[][] LPF_ecg(double[][] ecg) {
int DS_80MS = Math.round(80 * FS / 1000);
int DS_100MS = Math.round(100 * FS / 1000);
int DS_180MS = Math.round(180 * FS / 1000);
// preset
double[][] yLPF = ecg;
double Le = ecg.length;
double[] col_array = new double[ecg.length];
for (int c = 0; c < ecg[1].length; c++) {
for (int r = 0; r < ecg.length; r++) {
col_array[r] = ecg[r][c];
}
double[] d = diff_abs(col_array);
int[] vI = findD1Threshold(d);
int L = vI.length;
double T_refr = 1;
for (int j = 0; j < L; j++) {
if (vI[j] > T_refr) {
int t1 = vI[j] - DS_80MS;
int t2 = (int) Math.min(Le, t1 + DS_180MS);
if (t1 > 0) {
for (int r = t1; r < t2; r++) {
yLPF[r][c] = ecg[t1][c];
}
} else {
for (int r = 0; r <= t2; r++) {
yLPF[r][c] = ecg[t2][c];
}
}
T_refr = t2 + DS_100MS;
}
}
}
return yLPF;
}
private int[] diff(int[] ecg_column) {
// maintain count of non-zero differences
int[] diff_array = new int[ecg_column.length - 1];
for (int row = 1; row < ecg_column.length; row++) {
short df = (short) (ecg_column[row] - ecg_column[row - 1]);
diff_array[row - 1] = df;
}
return diff_array;
}
private double[] diff_abs(double[] ecg_column) {
D1_Thr_Count = 0;
// maintain count of non-zero differences
double[] diff_array = new double[ecg_column.length - 1];
for (int row = 1; row < ecg_column.length; row++) {
double df = (double) (ecg_column[row] - ecg_column[row - 1]);
diff_array[row - 1] = Math.abs(df);
if (diff_array[row - 1] > D1_Thr) {
D1_Thr_Count++;
}
}
return diff_array;
}
private int[] diff_find(short[] ecg_column, int data_number) {
int abs_thr_count = 0;
// maintain count of non-zero differences
int non_zero_count = 0;
short[] diff_array = new short[ecg_column.length - 1];
for (int row = 1; row < ecg_column.length; row++) {
short df = (short) (ecg_column[row] - ecg_column[row - 1]);
diff_array[row - 1] = df;
if (df != 0) {
non_zero_count++;
}
// also collect number of instances where abs(ecg) < 400
if (Math.abs(ecg_column[row]) < 400) {
abs_thr_count++;
}
}
// set number of abs_thr finds for future use
if (Math.abs(ecg_column[0]) < 400) {
abs_thr_count++;
}
Thr_Count_Array[data_number][1] = abs_thr_count;
int[] find_array = find(diff_array, non_zero_count);
return find_array;
}
// get all non-zero indices from diff array
private int[] find(short[] diff_array, int count) {
int[] find_array = new int[count];
int index_counter = 0;
int find_array_counter = 0;
for (int row = 0; row < diff_array.length; row++) {
if ((short) diff_array[row] != 0) {
find_array[find_array_counter] = index_counter;
find_array_counter++;
}
index_counter++;
}
index_counter++;
return find_array;
}
// get all non-zero indices from diff array
private int[] findThreshold(short[] diff_array, int count) {
int[] find_array = new int[count];
int index_counter = 0;
int find_array_counter = 0;
for (int row = 0; row < diff_array.length; row++) {
if (Math.abs(diff_array[row]) < SAT_THR) {
find_array[find_array_counter] = index_counter;
find_array_counter++;
}
index_counter++;
}
index_counter++;
return find_array;
}
private int[] findD1Threshold(double[] diff_array) {
int[] find_array = new int[D1_Thr_Count];
int index_counter = 0;
int find_array_counter = 0;
for (int row = 0; row < diff_array.length; row++) {
if (diff_array[row] > D1_Thr) {
find_array[find_array_counter] = index_counter;
find_array_counter++;
}
index_counter++;
}
index_counter++;
return find_array;
}
private int getMaxDiff(int[] find_array) {
int max_val = 0;
for (int i = 1; i < find_array.length; i++) {
int temp_val = find_array[i] - find_array[i - 1];
if (temp_val > max_val) {
max_val = temp_val;
}
}
return max_val;
}
private boolean detectFlat(short[] col_array, int c_number) {
int[] non_zero_indice_array = diff_find(col_array, c_number);
if (non_zero_indice_array.length > 0) {
int[] diff_array = diff(non_zero_indice_array);
int max = getMaxDiff(diff_array);
if (max > FS) {
return true;
} else {
return false;
}
} else {
return true;
}
}
private boolean detectSaturation(short[] col_array, int c_number) {
int[] below_threshold_array = findThreshold(col_array,
Thr_Count_Array[c_number][1]);
if (below_threshold_array.length > 0) {
int[] diff_array = diff(below_threshold_array);
int max = getMaxDiff(diff_array);
if (max > DW_Sat) {
return true;
} else {
return false;
}
} else {
return true;
}
}
private double[] filter(double[] ecg_unfilt) {
// order 6 filter
int order = 6;
double[] b = new double[order + 1];
double[] a = new double[order + 1];
double[] x = ecg_unfilt;
int nDataPoints = x.length;
b[0] = 0.0000000000000600231669922735;
b[1] = 0.000000000000360139001953641;
b[2] = 0.0000000000009003475048841025;
b[3] = 0.00000000000120046333984547;
b[4] = 0.0000000000009003475048841025;
b[5] = 0.000000000000360139001953641;
b[0] = 0.0000000000000600231669922735;
a[0] = 1.0000;
a[1] = -5.951447334984479;
a[2] = 14.758413952404254;
a[3] = -19.5191643872012;
a[4] = 14.521482981667374;
a[5] = -5.761891935232455;
a[6] = 0.952606723350349;
double[] y = new double[nDataPoints];
double[] filteredEcg_1 = filterECG(order, a, b, nDataPoints, x, y);
return filteredEcg_1;
}
private double[] filterECG(int order, double[] a, double[] b,
int nDataPoints, double[] x, double[] y) {
y[0] = b[0] * x[0];
for (int i = 0; i < order; i++) {
y[i] = 0.0;
for (int j = 0; j < i + 1; j++)
y[i] = y[i] + b[j] * x[i - j];
for (int j = 0; j < i; j++)
y[i] = y[i] - a[j + 1] * y[i - j - 1];
}
for (int i = order; i < nDataPoints; i++) {
y[i] = 0.0;
for (int j = 0; j < order + 1; j++)
y[i] = y[i] + b[j] * x[i - j];
for (int j = 0; j < order; j++)
y[i] = y[i] - a[j + 1] * y[i - j - 1];
}
return y;
}
}