/*
lib.cpp - ECG annotation library based on continuous (CWT) and fast (FWT) wavelet transforms.
Copyright (C) 2006 YURIY V. CHESNOKOV
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 Street, Fifth Floor, Boston, MA 02110-1301, USA.
You may contact the author by e-mail (chesnokov.yuriy@gmail.com or
chesnokov_yuriy@mail.ru) or postal mail (Unilever Centre for Molecular
Sciences Informatics, University Chemical Laboratory, Cambridge
University, Lensfield Road, Cambridge, CB2 1EW, UK)
*/
#include "stdafx.h"
#include "LIB.h"
//----------------------------------------------------------------------------
//////////////////////////////////////////////////////////////////////////////
////////////CLASS SIGNAL//////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
Signal::Signal(): data(0), fp(0),fpmap(0),lpMap(0),
SR(0.0),Bits(0),UmV(0),Lead(0),Len(0),
hh(0),mm(0),ss(0)
{
wcscpy(fname,L"");
}
Signal::~Signal()
{
if(vdata.size())
{
for(int i=0; i<(int)vdata.size(); i++)
{
data = vdata[i];
delete[] data;
}
}
}
//----------------------------------------------------------------------------
long double * Signal::ReadFile(char *name)
{
wchar_t ustr[256]=L"";
for(int i=0; i<(int)strlen(name); i++)
mbtowc(ustr+i,name+i, MB_CUR_MAX);
return ReadFile(ustr);
}
long double * Signal::ReadFile(wchar_t *name)
{
wcscpy(fname,name);
if(!GetInfo(name)) return 0;
if(dat)
{
if(!ReadDat()) return 0;
}
else
{
if(!ReadTxt()) //read text file
{
if(!ReadMit()) //read mit-bih file
return 0;
}
}
return data;
}
//----------------------------------------------------------------------------
bool Signal::GetInfo(wchar_t *name)
{
fp = CreateFileW(name,GENERIC_WRITE | GENERIC_READ,0,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
return false;
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,sizeof(DATAHDR),0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,sizeof(DATAHDR));
phdr = (PDATAHDR)lpMap;
if(!memcmp(phdr->hdr,"DATA",4))
{
Len = phdr->size;
SR = phdr->sr;
Bits = phdr->bits;
Lead = phdr->lead;
UmV = phdr->umv;
hh = phdr->hh;
mm = phdr->mm;
ss = phdr->ss;
dat = true;
}
else
dat = false;
CloseFile();
return true;
}
//----------------------------------------------------------------------------
bool Signal::ReadDat()
{
short tmp;
fp = CreateFileW(fname,GENERIC_WRITE|GENERIC_READ,0,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
return false;
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,0,0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,0);
phdr = (PDATAHDR)lpMap;
hdrs.push_back(*phdr);
lpf = (float *)lpMap;
lps = (short *)lpMap;
lpu = (unsigned short *)lpMap;
lpc = (char *)lpMap;
Len = phdr->size;
SR = phdr->sr;
Bits = phdr->bits;
Lead = phdr->lead;
UmV = phdr->umv;
lpf += sizeof(DATAHDR)/sizeof(float);
lps += sizeof(DATAHDR)/sizeof(short);
lpc += sizeof(DATAHDR);
data = new long double[Len];
vdata.push_back(data);
for(int i=0; i<Len; i++)
switch(Bits)
{
case 12: //212 format 12bit
if(i%2 == 0)
{
tmp = MAKEWORD(lpc[0],lpc[1]&0x0f);
if(tmp > 0x7ff)
tmp |= 0xf000;
data[i] = (long double)tmp/(long double)UmV;
}
else
{
tmp = MAKEWORD(lpc[2],(lpc[1]&0xf0)>>4);
if(tmp > 0x7ff)
tmp |= 0xf000;
data[i] = (long double)tmp/(long double)UmV;
lpc += 3;
}
break;
case 16: //16format
data[i] = (long double)lps[i]/(long double)UmV;
break;
case 0:
case 32: //32bit float
data[i] = (long double)lpf[i];
break;
}
CloseFile();
return true;
}
//----------------------------------------------------------------------------
bool Signal::ReadTxt()
{
vector<long double> vec1;
long double tmp;
int res;
if((in = _wfopen(fname,L"rt")) == 0)
return false;
for(;;)
{
res = fscanf(in,"%Lf",&tmp);
if(res == EOF || res == 0)
break;
else
vec1.push_back(tmp);
}
fclose(in);
Len = vec1.size();
if(Len < 2)
return false;
data = new long double[Len];
for(int i=0; i<Len; i++)
data[i] = vec1[i];
vdata.push_back(data);
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
hdr.size = Len;
hdr.umv = 1;
hdr.bits = 32;
hdrs.push_back(hdr);
return true;
}
//-----------------------------------------------------------------------------
bool Signal::ReadMit()
{
wchar_t hdrfile[256];
wcscpy(hdrfile,fname);
changeext(hdrfile,L".hea");
FILE *fh = _wfopen(hdrfile,L"rt");
if(!fh) return false;
if(parseHdr(fh))
{
fclose(fh);
phdr = &hdrs[0];
int lNum = (int)hdrs.size();
int size = phdr->size;
fp = CreateFileW(fname,GENERIC_WRITE|GENERIC_READ,0,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0);
if(fp==INVALID_HANDLE_VALUE || (GetFileSize(fp,0)!=(lNum*phdr->size*phdr->bits)/8))
return false;
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,0,0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,0);
lps = (short *)lpMap;
lpc = (char *)lpMap;
for(int i=0; i<lNum; i++)
{
data = new long double[phdr->size];
vdata.push_back(data);
}
short tmp;
for(int s=0; s<size; s++)
{
for(int n=0; n<lNum; n++)
{
data = vdata[n];
phdr = &hdrs[n];
switch(phdr->bits)
{
case 12: //212 format 12bit
if((s*lNum+n)%2 == 0)
{
tmp = MAKEWORD(lpc[0],lpc[1]&0x0f);
if(tmp > 0x7ff)
tmp |= 0xf000;
tmp -= phdr->bline;
data[s] = (long double)tmp/(long double)phdr->umv;
}
else
{
tmp = MAKEWORD(lpc[2],(lpc[1]&0xf0)>>4);
if(tmp > 0x7ff)
tmp |= 0xf000;
tmp -= phdr->bline;
data[s] = (long double)tmp/(long double)phdr->umv;
lpc += 3;
}
break;
case 16: //16format
data[s] = (long double)(*lps++ - phdr->bline)/(long double)phdr->umv;
break;
}
}
}
GetData(0);
CloseFile();
return true;
}
else
{
fclose(fh);
return false;
}
}
int Signal::GetData(int index)
{
if(!vdata.size()) return 0;
if(index > (int)vdata.size()-1) index=(int)vdata.size()-1;
else if(index < 0) index = 0;
phdr = &hdrs[index];
data = vdata[index];
SR = phdr->sr;
Lead = phdr->lead;
UmV = phdr->umv;
Bits = phdr->bits;
Len = phdr->size;
hh = phdr->hh;
mm = phdr->mm;
ss = phdr->ss;
return index;
}
int Signal::parseHdr(FILE *fh)
{
char leads[18][6] = {"I","II","III","aVR","aVL","aVF","v1","v2",
"v3","v4","v5","v6","MLI","MLII","MLIII","vX","vY","vZ"};
char str[10][256];
if(readline(fh, str[0]) <= 0)
return false;
int sNum,size;
float sr;
int res = sscanf(str[0],"%s %s %s %s",str[1],str[2],str[3],str[4]);
if(res != 4)
return 0;
sNum = atoi(str[2]);
sr = atof(str[3]);
size = atoi(str[4]);
int eNum=0;
for(int i=0; i<sNum; i++)
{
if(readline(fh, str[0]) <= 0)
return 0;
int umv,bits,bline;
memset(str[9],0,256);
res = sscanf(str[0],"%s %s %s %s %s %s %s %s %s",str[1],str[2],str[3],str[4],str[5],str[6],str[7],str[8],str[9]);
if(res < 5) return 0;
bits = atoi(str[2]);
umv = atoi(str[3]);
bline = atoi(str[5]);
int offs=strlen(str[1]),j=0;
for(j=0; j<offs; j++)
{
if(fname[(wcslen(fname)-offs)+j] != str[1][j]) //wctomb
break;
}
if(j == offs)
{
eNum++;
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
hdr.sr = sr;
hdr.bits = (bits == 212)? 12: bits;
hdr.umv = (umv == 0)? 200: umv;
hdr.bline = bline;
hdr.size = size;
for(int l=0; l<18; l++)
{
if(!stricmp(leads[l],str[9]))
{
hdr.lead = l+1;
break;
}
}
hdrs.push_back(hdr);
}
}
return eNum;
}
//-----------------------------------------------------------------------------
bool Signal::SaveFile(char *name, long double *buff, PDATAHDR hdr)
{
wchar_t ustr[256]=L"";
for(int i=0; i<(int)strlen(name); i++)
mbtowc(ustr+i,name+i, MB_CUR_MAX);
return SaveFile(ustr, buff, hdr);
}
bool Signal::SaveFile(wchar_t *name, long double *buff, PDATAHDR hdr)
{
int filesize;
int tmp;
switch(hdr->bits)
{
case 12:
if(hdr->size%2 != 0)
filesize = int((long double)(hdr->size+1)*1.5);
else
filesize = int((long double)(hdr->size)*1.5);
break;
case 16:
filesize = hdr->size*2;
break;
case 0:
case 32:
filesize = hdr->size*4;
break;
}
fp = CreateFileW(name,GENERIC_WRITE|GENERIC_READ,0,0,CREATE_ALWAYS,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
return false;
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,filesize+sizeof(DATAHDR),0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,filesize+sizeof(DATAHDR));
lpf = (float *)lpMap;
lps = (short *)lpMap;
lpu = (unsigned short *)lpMap;
lpc = (char *)lpMap;
memset(lpMap,0,filesize+sizeof(DATAHDR));
memcpy(lpMap,hdr,sizeof(DATAHDR));
lpf += sizeof(DATAHDR)/sizeof(float);
lps += sizeof(DATAHDR)/sizeof(short);
lpc += sizeof(DATAHDR);
for(unsigned int i=0; i<hdr->size; i++)
{
switch(hdr->bits)
{
case 12: //212 format 12bit
tmp = int(buff[i]*(long double)hdr->umv);
if(tmp > 2047) tmp = 2047;
if(tmp < -2048) tmp = -2048;
if(i%2 == 0)
{
lpc[0] = LOBYTE((short)tmp);
lpc[1] = 0;
lpc[1] = HIBYTE((short)tmp) & 0x0f;
}
else
{
lpc[2] = LOBYTE((short)tmp);
lpc[1] |= HIBYTE((short)tmp) << 4;
lpc += 3;
}
break;
case 16: //16format
tmp = int(buff[i]*(long double)hdr->umv);
if(tmp > 32767) tmp = 32767;
if(tmp < -32768) tmp = -32768;
*lps++ = (short)tmp;
break;
case 0:
case 32: //32bit float
*lpf++ = (float)buff[i];
break;
}
}
CloseFile();
return true;
}
//-----------------------------------------------------------------------------
bool Signal::ToTxt(wchar_t *name, long double *buff, int size)
{
in = _wfopen(name,L"wt");
if(in)
{
for(int i=0; i<size; i++)
fwprintf(in,L"%f\n",(float)buff[i]);
fclose(in);
return true;
}
else
return false;
}
void Signal::getfname(wchar_t *path, wchar_t *name)
{
int sl=0,dot=(int)wcslen(path);
int i;
for(i=0; i<(int)wcslen(path); i++)
{
if(path[i] == '.') break;
if(path[i] == '\\') break;
}
if(i >= (int)wcslen(path))
{
wcscpy(name,path);
return;
}
for(i=(int)wcslen(path)-1; i>=0; i--)
{
if(path[i] == '.')
dot = i;
if(path[i] == '\\')
{
sl = i+1;
break;
}
}
memcpy(name,&path[sl],(dot-sl)*2);
name[dot-sl]=0;
}
void Signal::changeext(wchar_t *path, wchar_t *ext)
{
for(int i=(int)wcslen(path)-1; i>0; i--)
{
if(path[i] == '.')
{
path[i] = 0;
wcscat(path, ext);
return;
}
}
wcscat(path, ext);
}
int Signal::readline(FILE *f, char *buff)
{
int res=0;
char *pbuff = buff;
while((short)res != EOF)
{
res = fgetc(f);
if(res == 0xD || res == 0xA)
{
if(pbuff == buff) continue;
*pbuff = 0;
return 1;
}
if((short)res != EOF)
{
*pbuff++ = (char)res;
}
}
return (short)res;
}
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
void Signal::mSecToTime(int msec, int &h, int &m, int &s, int &ms)
{
ms = msec%1000;
msec /= 1000;
if(msec < 60)
{
h = 0;
m = 0;
s = msec; //sec to final
}
else
{
long double tmp;
tmp = (long double)(msec%60) / 60;
tmp *= 60;
s = int(tmp);
msec /= 60;
if(msec < 60)
{
h = 0;
m = msec;
}
else
{
h = msec/60;
tmp = (long double)(msec%60)/60;
tmp *= 60;
m = int(tmp);
}
}
}
//-----------------------------------------------------------------------------
void Signal::MinMax(long double *buff, int size, long double &min, long double &max)
{
max=buff[0];
min=buff[0];
for(int i=1; i<size; i++)
{
if(buff[i] > max)max=buff[i];
if(buff[i] < min)min=buff[i];
}
}
//-----------------------------------------------------------------------------
void Signal::CloseFile()
{
if(lpMap)
{
UnmapViewOfFile(lpMap);
if(fpmap)
{
CloseHandle(fpmap);
if(fp && fp != INVALID_HANDLE_VALUE)
CloseHandle(fp);
}
}
}
///////////////////////////////////////////////////////////////////////////////
//-------MATH ROUTINES--------------------------------------------------------
long double Signal::Mean(long double *mass, int size)
{
long double mean = 0;
for(int i=0; i<size; i++)
mean += mass[i];
return mean /(long double)size;
}
long double Signal::Disp(long double *mass, int size)
{
long double mean,disp=0;
mean = Mean(mass,size);
for(int i=0; i<size; i++)
disp += (mass[i] - mean)*(mass[i] - mean);
return (sqrtl(disp/(long double)(size-1)));
}
//----------------------------------------------------------------------------
void Signal::nMinMax(long double *buff, int len, long double a, long double b)
{
long double min,max;
MinMax(buff,len,min,max);
for(int i=0; i<len; i++)
{
if(max-min)
buff[i] = (buff[i]-min)*((b-a)/(max-min)) + a;
else
buff[i] = a;
}
}
void Signal::nMean(long double *buff, int len)
{
long double mean = Mean(buff,len);
for(int i=0; i<len; i++)
buff[i] = buff[i]-mean;
}
void Signal::nZscore(long double *buff, int len)
{
long double mean = Mean(buff,len);
long double disp = Disp(buff,len);
if(disp==0.0) disp=1.0;
for(int i=0; i<len; i++)
buff[i] = (buff[i]-mean)/disp;
}
void Signal::nSoftmax(long double *buff, int len)
{
long double mean = Mean(buff,len);
long double disp = Disp(buff,len);
if(disp==0.0) disp=1.0;
for(int i=0; i<len; i++)
buff[i] = 1.0/( 1 + expl(-((buff[i]-mean)/disp)) );
}
void Signal::nEnergy(long double *buff, int len, int L)
{
long double enrg=0.0;
for(int i=0; i<len; i++)
enrg += powl(fabsl(buff[i]), (long double)L);
enrg = powl(enrg, 1.0/(long double)L);
if(enrg==0.0) enrg=1.0;
for(int i=0; i<len; i++)
buff[i] /= enrg;
}
//----------------------------------------------------------------------------
long double Signal::MINIMAX(long double *mass, int size)
{
return Disp(mass,size)*(0.3936 + 0.1829*logl((long double)size));
}
long double Signal::FIXTHRES(long double *mass, int size)
{
return Disp(mass,size)*sqrtl(2.0*logl((long double)size));
}
long double Signal::SURE(long double *mass, int size)
{
return Disp(mass,size)*sqrtl(2.0*logl((long double)size*logl((long double)size)));
}
void Signal::Denoise(long double *mass, int size, int window, int type, bool soft)
{
long double TH;
for(int i=0; i<size/window; i++)
{
switch(type)
{
case 0:
TH = MINIMAX(mass,window);
break;
case 1:
TH = FIXTHRES(mass,window);
break;
case 2:
TH = SURE(mass,window);
break;
}
if(soft)
SoftTH(mass, window, TH);
else
HardTH(mass, window, TH);
mass += window;
}
if(size%window > 5) //skip len=1
{
switch(type)
{
case 0:
TH = MINIMAX(mass,size%window);
break;
case 1:
TH = FIXTHRES(mass,size%window);
break;
case 2:
TH = SURE(mass,size%window);
break;
}
if(soft)
SoftTH(mass, size%window, TH);
else
HardTH(mass, size%window, TH);
}
}
void Signal::HardTH(long double *mass, int size, long double TH, long double l)
{
for(int i=0; i<size; i++)
if(fabsl(mass[i]) <= TH)
mass[i] *= l;
}
void Signal::SoftTH(long double *mass, int size, long double TH, long double l)
{
for(int i=0; i<size; i++)
{
if(fabsl(mass[i]) <= TH)
{
mass[i] *= l;
}
else
{
if(mass[i] > 0)
{
mass[i] -= TH*(1-l);
}
else
{
mass[i] += TH*(1-l);
}
}
}
}
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
bool Signal::AutoCov(long double *mass, int size)
{
long double *rk, mu;
int t;
rk = new long double[size];
mu = Mean(mass,size);
for(int k=0; k<size; k++)
{
rk[k] = 0;
t=0;
while(t+k != size)
{
rk[k] += (mass[t]-mu)*(mass[t+k]-mu);
t++;
}
rk[k] /= (long double)t; // rk[k] /= t ? autocovariance
}
for(int i=0; i<size; i++)
mass[i] = rk[i];
delete[] rk;
return true;
}
//-----------------------------------------------------------------------------
//////////////////////////////////////////////////////////////////////////////
////////////////////////CLASS CWT/////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
CWT::CWT(): type(0), cwtBuff(0),ReW(0),ImW(0)
{
}
CWT::~CWT()
{
if(ReW) free(ReW);
if(ImW) free(ImW);
if(cwtBuff) free(cwtBuff);
}
//----------------------------------------------------------------------------
float * CWT::CWTCreateFileHeader(wchar_t *name, PCWTHDR hdr, int index, long double w)
{
wchar_t tmp[256];
int filesize;
switch(index)
{
case 0:
wcscat(name,L"(mHat).w");
break;
case 1:
wcscat(name,L"(Inv).w");
break;
case 2:
wcscat(name,L"(Morl).w");
break;
case 3:
wcscat(name,L"(MPow).w");
break;
case 4:
wcscat(name,L"(MComp");
swprintf(tmp,L"%d",(int)w);
wcscat(name,tmp);
wcscat(name,L").w");
break;
case 5:
wcscat(name,L"(Gaussian).w");
break;
case 6:
wcscat(name,L"(1Gauss).w");
break;
case 7:
wcscat(name,L"(2Gauss).w");
break;
case 8:
wcscat(name,L"(3Gauss).w");
break;
case 9:
wcscat(name,L"(4Gauss).w");
break;
case 10:
wcscat(name,L"(5Gauss).w");
break;
case 11:
wcscat(name,L"(6Gauss).w");
break;
case 12:
wcscat(name,L"(7Gauss).w");
break;
}
if(hdr->type)
filesize = hdr->size * (int)ceill( (logl(hdr->fmax)+hdr->fstep-logl(hdr->fmin))/hdr->fstep );
else
filesize = hdr->size * (int)ceill( (hdr->fmax+hdr->fstep-hdr->fmin)/hdr->fstep );
filesize = sizeof(float)*filesize+sizeof(CWTHDR);
fp = CreateFileW(name,GENERIC_WRITE | GENERIC_READ,0,0,CREATE_ALWAYS,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
return 0;
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,filesize,0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,filesize);
lpf = (float *)lpMap;
memset(lpMap,0,filesize);
memcpy(lpMap,hdr,sizeof(CWTHDR));
return (lpf+sizeof(CWTHDR)/sizeof(float));
}
//----------------------------------------------------------------------------
float * CWT::CWTReadFile(wchar_t *name)
{
fp = CreateFileW(name,GENERIC_WRITE | GENERIC_READ,0,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
return 0;
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,0,0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,0);
phdr = (PCWTHDR)lpMap;
lpf = (float *)lpMap;
if(memcmp(phdr->hdr,"WLET",4))
return 0;
fmin = phdr->fmin;
fmax = phdr->fmax;
fstep = phdr->fstep;
Len = phdr->size;
SR = phdr->sr;
type = phdr->type;
return (lpf+sizeof(CWTHDR)/sizeof(float));
}
//----------------------------------------------------------------------------
long double CWT::HzToScale(long double f, long double sr, int index, long double w)
{
long double k;
switch(index)
{
case 0:
k = 0.22222 * sr;
break;
case 1:
k = 0.15833 * sr;
break;
case 2:
case 3:
k = sr;
break;
case 4:
k = sr * w * 0.1589;
break;
case 5: //Gauss
k = 0.2 * sr;
break;
case 6: //1Gauss
k = 0.16 * sr;
break;
case 7: //2Gauss
k = 0.224 * sr;
break;
case 8: //3Gauss
k = 0.272 * sr;
break;
case 9: //4Gauss
k = 0.316 * sr;
break;
case 10: //5Gauss
k = 0.354 * sr;
break;
case 11: //6Gauss
k = 0.388 * sr;
break;
case 12: //7Gauss
k = 0.42 * sr;
break;
default:
k = 0;
}
return (k/f);
}
int CWT::GetFreqRange()
{
if(type == 0)
return (int)((fmax+fstep-fmin)/fstep);
else
return (int)((logl(fmax)+fstep-logl(fmin))/fstep);
}
void CWT::ConvertName(wchar_t *name, int index, long double w)
{
wchar_t tmp[256];
switch(index)
{
case 0:
wcscat(name,L"(mHat).w");
break;
case 1:
wcscat(name,L"(Inv).w");
break;
case 2:
wcscat(name,L"(Morl).w");
break;
case 3:
wcscat(name,L"(MPow).w");
break;
case 4:
wcscat(name,L"(MComp");
swprintf(tmp,L"%d",(int)w);
wcscat(name,tmp);
wcscat(name,L").w");
break;
case 5:
wcscat(name,L"(Gaussian).w");
break;
case 6:
wcscat(name,L"(1Gauss).w");
break;
case 7:
wcscat(name,L"(2Gauss).w");
break;
case 8:
wcscat(name,L"(3Gauss).w");
break;
case 9:
wcscat(name,L"(4Gauss).w");
break;
case 10:
wcscat(name,L"(5Gauss).w");
break;
case 11:
wcscat(name,L"(6Gauss).w");
break;
case 12:
wcscat(name,L"(7Gauss).w");
break;
}
}
//----------------------------------------------------------------------------
void CWT::InitCWT(int sz, int index, long double w, long double sr)
{
cwtSize = sz;
if(sr)
SR = sr;
w0 = w;
ReW = (long double *)malloc(sizeof(long double)*(2*cwtSize-1));
ImW = (long double *)malloc(sizeof(long double)*(2*cwtSize-1));
cwtBuff = (long double *)malloc(sizeof(long double)*(cwtSize));
cwIndex = index;
for(int i=0; i<2*cwtSize-1; i++)
{
ReW[i]=0;
ImW[i]=0;
}
}
//----------------------------------------------------------------------------
void CWT::CloseCWT()
{
if(ReW) { free(ReW); ReW=0; };
if(ImW) { free(ImW); ImW=0; };
if(cwtBuff) {free(cwtBuff); cwtBuff=0; };
}
//----------------------------------------------------------------------------
long double *CWT::CWTrans(long double *buff, long double freq, bool mr,long double lv,long double rv)
{
mirror = mr;
lval=lv;
rval=rv;
prec = false;
precision=0; //0,0000001 float prsision
long double t,sn,cs,scale;
scale = HzToScale(freq,SR,cwIndex,w0);
///////////wavelet calculation//////////////////////////////////////////////////////////////////
///////// center = csigLen-1 in wavelet mass////////////////////////////////////////////////////
for(int i=0; i<cwtSize; i++) //positive side
{
t = ((long double)i)/scale;
if(cwIndex > 1 && cwIndex < 4)
{
sn = sinl(6.28*t);
cs = cosl(6.28*t);
}
if(cwIndex == 4)
{
sn = sinl(w0*t);
cs = cosl(w0*t);
}
switch(cwIndex)
{
case 0:
ReW[(cwtSize-1)+i] = expl(-t*t/2)*(-t*t+1);
break;
case 1:
ReW[(cwtSize-1)+i] = t*expl(-t*t/2);
break;
case 2:
ReW[(cwtSize-1)+i] = expl(-t*t/2)*(cs-sn);
break;
case 3:
ReW[(cwtSize-1)+i] = expl(-t*t/2)* cs;
ImW[(cwtSize-1)+i] = expl(-t*t/2)* sn;
break;
case 4:
ReW[(cwtSize-1)+i] = expl(-t*t/2)* (cs - expl(-w0*w0/2));
ImW[(cwtSize-1)+i] = expl(-t*t/2)* (sn - expl(-w0*w0/2));
break;
case 5:
ReW[(cwtSize-1)+i] = expl(-t*t/2);
break;
case 6:
ReW[(cwtSize-1)+i] = -t*expl(-t*t/2);
break;
case 7:
ReW[(cwtSize-1)+i] = (t*t-1)*expl(-t*t/2);
break;
case 8:
ReW[(cwtSize-1)+i] = (2*t+t-t*t*t)*expl(-t*t/2);
break;
case 9:
ReW[(cwtSize-1)+i] = (3-6*t*t+t*t*t*t)*expl(-t*t/2);
break;
case 10:
ReW[(cwtSize-1)+i] = (-15*t+10*t*t*t-t*t*t*t*t)*expl(-t*t/2);
break;
case 11:
ReW[(cwtSize-1)+i] = (-15+45*t*t-15*t*t*t*t+t*t*t*t*t*t)*expl(-t*t/2);
break;
case 12:
ReW[(cwtSize-1)+i] = (105*t-105*t*t*t+21*t*t*t*t*t-t*t*t*t*t*t*t)*expl(-t*t/2);
break;
}
if(fabsl(ReW[(cwtSize-1)+i]) < 0.0000001)
precision++;
if(precision > 15)
{
precision = i;
prec = true;
break;
}
}
if(prec == false)
precision = cwtSize;
for(int i=-(precision-1); i<0; i++) //negative side
{
t = ((long double)i)/scale;
if(cwIndex > 1 && cwIndex < 4)
{
sn = sinl(6.28*t);
cs = cosl(6.28*t);
}
if(cwIndex == 4)
{
sn = sinl(w0*t);
cs = cosl(w0*t);
}
switch(cwIndex)
{
case 0:
ReW[(cwtSize-1)+i] = expl(-t*t/2)*(-t*t+1);
break;
case 1:
ReW[(cwtSize-1)+i] = t*expl(-t*t/2);
break;
case 2:
ReW[(cwtSize-1)+i] = expl(-t*t/2)*(cs-sn);
break;
case 3:
ReW[(cwtSize-1)+i] = expl(-t*t/2)* cs;
ImW[(cwtSize-1)+i] = expl(-t*t/2)* sn;
break;
case 4:
ReW[(cwtSize-1)+i] = expl(-t*t/2)* (cs - expl(-w0*w0/2));
ImW[(cwtSize-1)+i] = expl(-t*t/2)* (sn - expl(-w0*w0/2));
break;
case 5:
ReW[(cwtSize-1)+i] = expl(-t*t/2); //gauss
break;
case 6:
ReW[(cwtSize-1)+i] = -t*expl(-t*t/2); //gauss1
break;
case 7:
ReW[(cwtSize-1)+i] = (t*t-1)*expl(-t*t/2); //gauss2
break;
case 8:
ReW[(cwtSize-1)+i] = (2*t+t-t*t*t)*expl(-t*t/2); //gauss3
break;
case 9:
ReW[(cwtSize-1)+i] = (3-6*t*t+t*t*t*t)*expl(-t*t/2); //gauss4
break;
case 10:
ReW[(cwtSize-1)+i] = (-15*t+10*t*t*t-t*t*t*t*t)*expl(-t*t/2); //gauss5
break;
case 11:
ReW[(cwtSize-1)+i] = (-15+45*t*t-15*t*t*t*t+t*t*t*t*t*t)*expl(-t*t/2); //gauss6
break;
case 12:
ReW[(cwtSize-1)+i] = (105*t-105*t*t*t+21*t*t*t*t*t-t*t*t*t*t*t*t)*expl(-t*t/2); //gauss7
break;
}
}
///////end wavelet calculations////////////////////////////////////////////
csig = buff;
for (int x = 0; x < cwtSize; x++)
cwtBuff[x] = cWave(x,scale);
return cwtBuff;
}
//----------------------------------------------------------------------------
long double CWT::cWave(int x, long double scale)
{
long double res=0, Re=0, Im=0;
for(int t = 0; t < cwtSize; t++) //main
{
if(prec==true)
{
if(t < x-precision)t=x-(precision-1);//continue;
if(t >= precision+x)break;
}
Re += ReW[((cwtSize-1)-x) + t] * csig[t];
if(cwIndex == 3 || cwIndex == 4)
Im += ImW[((cwtSize-1)-x) + t] * csig[t];
}
////////////////////boundaries///////////////////////////////////////////////
int p=0;
for(int i=(cwtSize-precision); i<(cwtSize-1)-x; i++) // Left edge calculations
{
if(mirror)
{
Re += ReW[i]*csig[(cwtSize-1)-i-x]; //mirror
}
else
{
if(lval)
Re += ReW[i]*lval;
else
Re += ReW[i]*csig[0];
}
if(cwIndex == 3 || cwIndex == 4) //Im part for complex wavelet
{
if(mirror)
{
Im += ImW[i]*csig[(cwtSize-1)-i-x];
}
else
{
if(lval)
Im += ImW[i]*lval;
else
Im += ImW[i]*csig[0];
}
}
}
int q=0;
for(int i=2*cwtSize-(x+1); i<cwtSize+precision-1; i++) // Right edge calculations
{
if(mirror)
Re += ReW[i]*csig[(cwtSize-2)-q]; //mirror
else
{
if(rval)
Re += ReW[i]*rval;
else
Re += ReW[i]*csig[cwtSize-1];
}
if(cwIndex == 3 || cwIndex == 4)
{
if(mirror)
{
Im += ImW[i]*csig[(cwtSize-2)-q];
}
else
{
if(rval)
Im += ImW[i]*rval;
else
Im += ImW[i]*csig[cwtSize-1];
}
}
q++;
}
////////////////////boundaries///////////////////////////////////////////////
switch(cwIndex)
{
case 2:
res = (1/sqrtl(6.28)) * Re;
break;
case 3:
res = sqrtl(Re*Re + Im*Im);
res *= (1/sqrtl(6.28));
break;
case 4:
res = sqrtl(Re*Re + Im*Im);
res *= (1/powl(3.14,0.25));
break;
default:
res = Re;
}
res = (1/sqrtl(scale)) * res;
return res;
}
//----------------------------------------------------------------------------*/
//////////////////////////////////////////////////////////////////////////////
////////////////////////CLASS FWT/////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
FWT::FWT(): tH(0),tG(0),H(0),G(0),
thL(0),tgL(0),hL(0),gL(0), thZ(0),tgZ(0),hZ(0),gZ(0),
fwtBuff(0),tempBuff(0),
J(0),Jnumbs(0)
{
}
FWT::~FWT()
{
if(tH) delete[] tH;
if(tG) delete[] tG;
if(H) delete[] H;
if(G) delete[] G;
if(fwtBuff) free(fwtBuff);
if(tempBuff) free(tempBuff);
if(Jnumbs) delete[] Jnumbs;
}
//----------------------------------------------------------------------------
//////////////////// filters init ///////////////////////////////////////////////////
bool FWT::InitFWT(wchar_t *fltname, long double *sig, int len)
{
filter = _wfopen(fltname,L"rt");
if(filter)
{
tH = loadfilter(thL, thZ);
tG = loadfilter(tgL, tgZ);
H = loadfilter(hL, hZ);
G = loadfilter(gL, gZ);
fclose(filter);
size = len;
sigLen = len;
fwtBuff = (long double *)malloc(sizeof(long double)*len);
tempBuff = (long double *)malloc(sizeof(long double)*len);
lo = tempBuff;
hi = tempBuff+len;
for(int i=0; i<len; i++)
fwtBuff[i] = sig[i];
memset(tempBuff,0,sizeof(long double)*len);
J = 0;
return true;
}
else
return false;
}
long double *FWT::loadfilter(int &L, int &Z)
{
fscanf(filter,"%d",&L);
fscanf(filter,"%d",&Z);
long double *flt = new long double[L];
for(int i=0; i<L; i++)
fscanf(filter,"%Lf",&flt[i]);
return flt;
}
void FWT::CloseFWT()
{
if(tH) {delete[] tH; tH=0;}
if(tG) {delete[] tG; tG=0;}
if(H) {delete[] H; H=0;}
if(G) {delete[] G; G=0;}
if(fwtBuff) {free(fwtBuff); fwtBuff=0;}
if(tempBuff) {free(tempBuff); tempBuff=0;}
if(Jnumbs) {delete[] Jnumbs; Jnumbs = 0;}
}
//////////////////// filters init ///////////////////////////////////////////////////
//----------------------------------------------------------------------------
void FWT::fwtHiLoTrans()
{
int n;
long double s,d;
for(int k=0; k<size/2; k++)
{
s = 0;
d = 0;
for(int m=-thZ; m<thL-thZ; m++)
{
n = 2*k + m;
if(n < 0) n = 0 - n;
if(n >= size) n -= (2 + n-size);
s += tH[m+thZ] * fwtBuff[n];
}
for(int m=-tgZ; m<tgL-tgZ; m++)
{
n = 2*k + m;
if(n < 0) n = 0 - n;
if(n >= size) n -= (2 + n-size);
d += tG[m+tgZ] * fwtBuff[n];
}
lo[k] = s;
hi[k] = d;
}
for(int i=0; i<sigLen; i++)
fwtBuff[i] = tempBuff[i];
}
void FWT::fwtTrans(int scales)
{
for(int j=0; j<scales; j++)
{
hi -= size/2;
fwtHiLoTrans();
size /= 2;
J++;
}
}
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
void FWT::fwtHiLoSynth()
{
int n;
long double s2k, s2k1;
for(int i=0; i<sigLen; i++)
tempBuff[i] = fwtBuff[i];
for(int k=0; k<size; k++)
{
s2k = 0;
s2k1 = 0;
for(int m=-hZ; m<hL-hZ; m++) //s2k and s2k1 for H[]
{
n = k-m;
if(n < 0) n = 0 - n;
if(n >= size) n -= (2 + n-size);
if(2*m >= -hZ && 2*m < hL-hZ)
s2k += H[(2*m)+hZ]*lo[n];
if((2*m+1) >= -hZ && (2*m+1) < hL-hZ)
s2k1 += H[(2*m+1)+hZ]*lo[n];
}
for(int m=-gZ; m<gL-gZ; m++) //s2k and s2k1 for G[]
{
n = k-m;
if(n < 0) n = 0 - n;
if(n >= size) n -= (2 + n-size);
if(2*m >= -gZ && 2*m < gL-gZ)
s2k += G[(2*m)+gZ]*hi[n];
if((2*m+1) >= -gZ && (2*m+1) < gL-gZ)
s2k1 += G[(2*m+1)+gZ]*hi[n];
}
fwtBuff[2*k] = 2.0*s2k;
fwtBuff[2*k+1] = 2.0*s2k1;
}
}
void FWT::fwtSynth(int scales)
{
for(int j=0; j<scales; j++)
{
fwtHiLoSynth();
hi += Jnumbs[j];
size *= 2;
J--;
}
}
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
int *FWT::GetJnumbs(int j, int len)
{
if(Jnumbs) delete[] Jnumbs;
Jnumbs = new int[j];
for(int i=0; i<j; i++)
Jnumbs[i] = len/(int)powl(2,(long double)(j-i));
return Jnumbs;
}
//-----------------------------------------------------------------------------
void FWT::HiLoNumbs(int j, int len, int &hinum, int &lonum)
{
lonum = 0;
hinum = 0;
for(int i=0; i<j; i++)
{
len /= 2;
hinum += len;
}
lonum = len;
}
//-----------------------------------------------------------------------------
bool FWT::fwtSaveFile(wchar_t *name, long double *hipass, long double *lopass, PFWTHDR hdr)
{
int filesize;
short tmp;
HiLoNumbs(hdr->J,hdr->size,hiNum,loNum);
switch(hdr->bits)
{
case 12:
if((hiNum+loNum)%2 != 0)
filesize = (long double)((hiNum+loNum)+1)*1.5;
else
filesize = (long double)(hiNum+loNum)*1.5;
break;
case 16:
filesize = (hiNum+loNum)*2;
break;
case 0:
case 32:
filesize = (hiNum+loNum)*4;
break;
}
fp = CreateFileW(name,GENERIC_WRITE | GENERIC_READ,0,0,CREATE_ALWAYS,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
{
fp = 0;
return false;
}
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,filesize+sizeof(FWTHDR),0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,filesize+sizeof(FWTHDR));
lpf = (float *)lpMap;
lps = (short *)lpMap;
lpc = (char *)lpMap;
memset(lpMap,0,filesize+sizeof(FWTHDR));
memcpy(lpMap,hdr,sizeof(FWTHDR));
lpf += sizeof(FWTHDR)/sizeof(float);
lps += sizeof(FWTHDR)/sizeof(short);
lpc += sizeof(FWTHDR);
for(int i=0; i<hiNum+loNum; i++)
{
if(i < hiNum)
tmp = short( hipass[i]*(long double)hdr->umv );
else
tmp = short( lopass[i-hiNum]*(long double)hdr->umv );
switch(hdr->bits)
{
case 12:
if(i%2 == 0)
{
lpc[0] = LOBYTE(tmp);
lpc[1] = 0;
lpc[1] = HIBYTE(tmp) & 0x0f;
}
else
{
lpc[2] = LOBYTE(tmp);
lpc[1] |= HIBYTE(tmp) << 4;
lpc += 3;
}
break;
case 16: //16format
*lps++ = tmp;
break;
case 0:
case 32:
if(i < hiNum) //32bit float
*lpf++ = (float)hipass[i];
else
*lpf++ = (float)lopass[i-hiNum];
break;
}
}
CloseFile();
return true;
}
//----------------------------------------------------------------------------
bool FWT::fwtReadFile(wchar_t *name, char *appdir)
{
short tmp;
fp = CreateFileW(name,GENERIC_WRITE | GENERIC_READ,0,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0);
if(fp == INVALID_HANDLE_VALUE)
{
fp = 0;
return false;
}
fpmap = CreateFileMapping(fp,0,PAGE_READWRITE,0,0,0);
lpMap = MapViewOfFile(fpmap,FILE_MAP_WRITE,0,0,0);
phdr = (PFWTHDR)lpMap;
lpf = (float *)lpMap;
lps = (short *)lpMap;
lpc = (char *)lpMap;
Len = phdr->size;
SR = phdr->sr;
Bits = phdr->bits;
Lead = phdr->lead;
UmV = phdr->umv;
strcpy(fname,phdr->wlet);
J = phdr->J;
lpf += sizeof(FWTHDR)/sizeof(float);
lps += sizeof(FWTHDR)/sizeof(short);
lpc += sizeof(FWTHDR);
HiLoNumbs(J,Len,hiNum,loNum);
sigLen = loNum+hiNum;
size = loNum;
fwtBuff = new long double[sigLen];
tempBuff = new long double[sigLen];
memset(tempBuff,0,sizeof(long double)*sigLen);
lo = fwtBuff;
hi = fwtBuff+loNum;
for(int i=0; i<sigLen; i++)
{
switch(Bits)
{
case 12: //212 format 12bit
if(i%2 == 0)
{
tmp = MAKEWORD(lpc[0],lpc[1]&0x0f);
if(tmp > 0x7ff)
tmp |= 0xf000;
}
else
{
tmp = MAKEWORD(lpc[2],(lpc[1]&0xf0)>>4);
if(tmp > 0x7ff)
tmp |= 0xf000;
lpc += 3;
}
if(i<hiNum)
{hi[i] = (long double)tmp/(long double)UmV;}
else
{lo[i-hiNum] = (long double)tmp/(long double)UmV;}
break;
case 16: //16format
if(i<hiNum)
{hi[i] = (long double)lps[i]/(long double)UmV;}
else
{lo[i-hiNum] = (long double)lps[i]/(long double)UmV;}
break;
case 0:
case 32: //32bit float
if(i<hiNum)
{hi[i] = (long double)lpf[i];}
else
{lo[i-hiNum] = (long double)lpf[i];}
break;
}
}
lo = tempBuff;
hi = tempBuff+loNum;
if(appdir) //load filter for synthesis
{
char flt[256];
strcpy(flt,appdir);
strcat(flt,"\\filters\\");
strcat(flt,fname);
strcat(flt,".flt");
filter = fopen(flt,"rt");
if(filter)
{
tH = loadfilter(thL, thZ);
tG = loadfilter(tgL, tgZ);
H = loadfilter(hL, hZ);
G = loadfilter(gL, gZ);
fclose(filter);
}
else
{
CloseFile();
return false;
}
}
CloseFile();
return true;
}
//-----------------------------------------------------------------------------
//////////////////////////////////////////////////////////////////////////////
////////////////////////CLASS ANN/////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
Annotation::Annotation(PANNHDR p): qrsNum(0),annNum(0),auxNum(0),
ANN(0),qrsANN(0),AUX(0)
{
if(p)
{
memcpy(&ahdr, p, sizeof(ANNHDR));
}
else //defaults
{
ahdr.minbpm = 40;
ahdr.maxbpm = 200;
ahdr.minQRS = 0.04; //min QRS duration
ahdr.maxQRS = 0.2; //max QRS duration
ahdr.minUmV = 0.2; //min UmV of R,S peaks
ahdr.minPQ = 0.07; //min PQ duration
ahdr.maxPQ = 0.20; //max PQ duration
ahdr.minQT = 0.21; //min QT duration
ahdr.maxQT = 0.48; //max QT duration
ahdr.pFreq = 9.0; //cwt Hz for P wave
ahdr.tFreq = 3.0; //cwt Hz for T wave
}
}
Annotation::~Annotation()
{
if(ANN)
{
for (int i = 0; i < annNum; i++)
delete[] ANN[i];
delete[] ANN;
}
if(qrsANN)
{
for (int i = 0; i < qrsNum; i++)
delete[] qrsANN[i];
delete[] qrsANN;
}
if(AUX)
{
for (int i = 0; i < auxNum; i++)
delete[] AUX[i];
delete[] AUX;
}
}
//-----------------------------------------------------------------------------
// 10Hz cwt trans of signal
// spectrum in cwt class
// create qrsANN array with qrsNum records num of heart beats = qrsNum/2
//
int ** Annotation::GetQRS(long double *mass, int len, long double sr, wchar_t *fltdir, int stype)
{
long double *pmass = (long double *)malloc(len*sizeof(long double));
for(int i=0; i<len; i++)
pmass[i] = mass[i];
if(fltdir)
{
if(flt30Hz(pmass,len,sr,fltdir,stype) == false) //pmass filed with filterd signal
{
delete[] pmass;
return 0;
}
// 2 continuous passes of CWT filtering
/*if(flt30Hz(pmass,len,sr,fltdir,qNORM) == false) //pmass filed with filterd signal
return 0;
if(stype == qNOISE)
if(flt30Hz(pmass,len,sr,fltdir,qNORM) == false)
return 0;*/
}
long double eCycle = (60.0/(long double)ahdr.maxbpm)-ahdr.maxQRS; //secs
if(int(eCycle*sr) <= 0)
{
eCycle = 0.1;
ahdr.maxbpm = int(60.0/(ahdr.maxQRS+eCycle));
}
//////////////////////////////////////////////////////////////////////////////
int lqNum=0;
vector <int> QRS; //clean QRS detected
int add=0;
while(pmass[add]) add+=int(0.1*sr); //skip QRS in begining
while(pmass[add]==0) add++; //get 1st QRS
QRS.push_back(add-1);
/////////////////////////////////////////////////////////////////////////////
for(int m=add; m<len; m++) //MAX 200bpm [0,3 sec min cario len]
{
m += int(ahdr.maxQRS*sr); //smpl + 0,20sec [0,20 max QRS length]
if(m>=len) m=len-1;
add=0;
//noise checking///////////////////////////////////////
if(m + int(eCycle*sr) >= len) //near end of signal
{
QRS.pop_back();
break;
}
if(chkNoise(&pmass[m],eCycle*sr)) //smpl(0.10sec)+0,20sec in noise
{
if(lqNum != (int)QRS.size()-1)
MA.push_back(QRS[QRS.size()-1]); //push MA noise location
QRS.pop_back();
lqNum = QRS.size();
//find for next possible QRS start
while(chkNoise(&pmass[m],eCycle*sr))
{
m += int(eCycle*sr);
if(m >= len-int(eCycle*sr)) break; //end of signal
}
if(m >= len-int(eCycle*sr)) break;
else
{
while(pmass[m]==0)
{
m++;
if(m >= len) break;
}
if(m >= len) break;
else
{
QRS.push_back(m-1);
continue;
}
}
}
////////////////////////////////////////////////////////
while(pmass[m-add]==0.0) add++; //find back for QRS end
if((m-add+1) - QRS[QRS.size()-1] > ahdr.minQRS*sr) //QRS size > 0.04 sec
QRS.push_back(m-add+2); //QRS end
else
QRS.pop_back();
m += int(eCycle*sr); //smpl + [0,20+0,10]sec 200bpm MAX
if(len-m < int(sr/2)) break;
while(pmass[m]==0 && (len-m)>=int(sr/2)) //find nearest QRS
m++;
if(len-m<int(sr/2)) break; //end of data
QRS.push_back(m-1); //QRS begin
}
/////////////////////////////////////////////////////////////////////////////
delete[] pmass;
qrsNum = QRS.size()/2;
if(qrsNum > 0) // 46: ?
{ // 1: N -1: nodata in **AUX
qrsANN = new int*[2*qrsNum]; // [samps] [type] [?aux data]
for(int i=0; i<2*qrsNum; i++)
qrsANN[i] = new int[3];
for(int i=0; i<2*qrsNum; i++)
{
qrsANN[i][0] = QRS[i]; //samp
if(i%2 == 0)
qrsANN[i][1] = 1; //type N
else
{
qrsANN[i][1] = 40; //type QRS)
//if( (qrsANN[i][0]-qrsANN[i-1][0]) >= int(sr*0.12) || (qrsANN[i][0]-qrsANN[i-1][0]) <= int(sr*0.03) )
// qrsANN[i-1][1] = 46; //N or ? beat (0.03?-0.12secs)
}
qrsANN[i][2] = -1; //no aux data
}
return qrsANN;
}
else
return 0;
}
bool Annotation::chkNoise(long double *mass, int window)
{
for(int i=0; i<window; i++)
if(mass[i]) return true;
return false;
}
bool Annotation::flt30Hz(long double *mass, int len, long double sr, wchar_t *fltdir, int stype)
{
CWT cwt;
///////////CWT 10Hz transform//////////////////////////////////////////
cwt.InitCWT(len,6,0,sr); //gauss1 wavelet 6-index
long double *pspec = cwt.CWTrans(mass,13); //10-13 Hz transform?
for(int i=0; i<len; i++)
mass[i] = pspec[i];
cwt.CloseCWT();
//debug
//ToTxt(L"10hz.txt",mass,len);
wchar_t flt[256]=L"";
wcscpy(flt,fltdir);
switch(stype)
{
case qNORM:
wcscat(flt,L"\\inter1.flt");
break;
case qNOISE:
for(int i=0; i<len; i++) //ridges
mass[i] *= (fabsl(mass[i])/2.0);
wcscat(flt,L"\\bior13.flt");
break;
}
FWT fwt;
////////////FWT 0-30Hz removal//////////////////////////////////////////
if(fwt.InitFWT(flt,mass,len) == false)
return false;
int J = ceill(log2(sr/23.0))-2;
//trans///////////////////////////////////////////////////
fwt.fwtTrans(J);
int *Jnumbs = fwt.GetJnumbs(J,len);
int hinum,lonum;
fwt.HiLoNumbs(J,len,hinum,lonum);
long double *lo = fwt.getSpec();
long double *hi = fwt.getSpec() + (len-hinum);
int window; //2.0sec window
for(int j=J; j>0; j--)
{
window = (2.0*sr)/powl(2.0,(long double)j); //2.0sec interval
Denoise(hi,Jnumbs[J-j],window,0,false); //hard,MINIMAX denoise [30-...Hz]
hi += Jnumbs[J-j];
}
for(int i=0; i<lonum; i++) //remove [0-30Hz]
lo[i] = 0.0;
//synth/////////////////////////////
fwt.fwtSynth(J);
for(int i=0; i<fwt.getSpecSize(); i++)
mass[i] = lo[i];
for(int i=len-(len-fwt.getSpecSize()); i<len; i++)
mass[i] = 0.0;
fwt.CloseFWT();
//debug
//ToTxt(L"10hz(intr1).txt",mass,len);
return true;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// find ectopic beats in HRV data
void Annotation::GetECT(int **ann, int qrsnum, long double sr)
{
if(qrsnum < 3)
return;
vector<long double> RRs;
for(int n=0; n<qrsnum-1; n++)
RRs.push_back((long double)(ann[n*2+2][0] - ann[n*2][0])/sr); //qrsNum-1 rr's
RRs.push_back(RRs[RRs.size()-1]);
// [RR1 RR2 RR3] RR2 beat classification
long double rr1,rr2,rr3;
for(int n=-2; n<(int)RRs.size()-2; n++)
{
if(n == -2)
{
rr1 = RRs[1]; //2
rr2 = RRs[0]; //1
rr3 = RRs[0]; //0
}
else if(n == -1)
{
rr1 = RRs[1];
rr2 = RRs[0];
rr3 = RRs[1];
}
else
{
rr1 = RRs[n];
rr2 = RRs[n+1];
rr3 = RRs[n+2];
}
if(60.0/rr1 < ahdr.minbpm || 60.0/rr1 > ahdr.maxbpm) //if RR's within 40-200bpm
continue;
if(60.0/rr2 < ahdr.minbpm || 60.0/rr2 > ahdr.maxbpm) //if RR's within 40-200bpm
continue;
if(60.0/rr3 < ahdr.minbpm || 60.0/rr3 > ahdr.maxbpm) //if RR's within 40-200bpm
continue;
if(1.15*rr2 < rr1 && 1.15*rr2 < rr3)
{
ann[n*2+4][1] = 46;
continue;
}
if( fabsl(rr1-rr2)<0.3 && rr1<0.8 && rr2<0.8 && rr3>2.4*(rr1+rr2) )
{
ann[n*2+4][1] = 46;
continue;
}
if( fabsl(rr1-rr2)<0.3 && rr1<0.8 && rr2<0.8 && rr3>2.4*(rr2+rr3) )
{
ann[n*2+4][1] = 46;
continue;
}
}
}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
//out united annotation with QRS PT////////////////////////////////////////////
// **ann [PQ,JP] pairs
int ** Annotation::GetPTU(long double *mass, int len, long double sr, wchar_t *fltdir, int **ann, int qrsnum)
{
int size, annPos;
int T1=-1,T=-1,T2=-1, Twaves=0;
int P1=-1,P=-1,P2=-1, Pwaves=0;
CWT cwt;
vector <int> Pwave;
vector <int> Twave; //Twave [ ( , T , ) ]
long double *buff; // [ 0 , 0 , 0 ] 0 no Twave
long double min,max; //min,max for gaussian1 wave, center is zero crossing
long double rr; //rr interval
bool sign;
bool twave,pwave;
int add = 0;//int(sr*0.04); //prevent imprecise QRS end detection
int maNum=0;
for(int n=0; n<qrsnum-1; n++)
{
annPos = ann[n*2+1][0]; //i
size = ann[n*2+2][0]-ann[n*2+1][0]; //i size of (QRS) <----> (QRS)
bool maNs=false;
for(int i=maNum; i<(int)MA.size(); i++)
{
if(MA[i]>ann[n*2+1][0] && MA[i]<ann[n*2+2][0])
{
Pwave.push_back(0);
Pwave.push_back(0);
Pwave.push_back(0);
Twave.push_back(0);
Twave.push_back(0);
Twave.push_back(0);
maNum++;
maNs=true;
break;
}
}
if(maNs) continue;
rr = (long double)(ann[n*2+2][0]-ann[n*2][0])/sr;
if(60.0/rr < ahdr.minbpm || 60.0/rr > ahdr.maxbpm-20) //check if normal RR interval (40bpm - 190bpm)
{
Pwave.push_back(0);
Pwave.push_back(0);
Pwave.push_back(0);
Twave.push_back(0);
Twave.push_back(0);
Twave.push_back(0);
continue;
}
///////////////search for TWAVE///////////////////////////////////////////////////////////
if(sr*ahdr.maxQT-(ann[n*2+1][0]-ann[n*2+0][0]) > size-add)
size = size - add;
else
size = sr*ahdr.maxQT-(ann[n*2+1][0]-ann[n*2+0][0]) - add;
//long double avg = Mean(mass+annPos+add,size); //avrg extension on boundaries
//long double lvl,rvl;
//lvl = mass[annPos+add];
//rvl = mass[annPos+add+size-1];
cwt.InitCWT(size,6,0,sr); //6-Gauss1 wlet
buff = cwt.CWTrans(mass+annPos+add,ahdr.tFreq);//,false,lvl,rvl); //3Hz transform buff = size-2*add
//cwt.ToTxt(L"debugS.txt",mass+annPos+add,size); //T wave
//cwt.ToTxt(L"debugC.txt",buff,size); //T wave spectrum
MinMax(buff,size,min,max);
for(int i=0; i<size; i++)
{
if(buff[i]==min) T1 = i+ annPos + add;
if(buff[i]==max) T2 = i+ annPos + add;
}
if(T1>T2)swap(T1,T2);
//additional constraints on [T1 T T2] duration, symmetry, QT interval
twave = false;
if((buff[T1-annPos-add]<0 && buff[T2-annPos-add]>0) || (buff[T1-annPos-add]>0 && buff[T2-annPos-add]<0))
twave = true;
if(twave)
{
if((long double)(T2-T1)>=0.09*sr)// && (long double)(T2-T1)<=0.24*sr) //check for T wave duration
{
twave = true; //QT interval = .4*sqrt(RR)
if((long double)(T2-ann[n*2+0][0])>=ahdr.minQT*sr && (long double)(T2-ann[n*2+0][0])<=ahdr.maxQT*sr)
twave = true;
else
twave = false;
}
else
twave = false;
}
if(twave)
{
if(buff[T1-annPos-add] > 0) sign = true;
else sign = false;
for(int i=T1-annPos-add; i<T2-annPos-add; i++)
{
if(sign)
{ if(buff[i] > 0) continue; }
else
{ if(buff[i] < 0) continue; }
T = i+annPos + add;
break;
}
//check for T wave symetry//////////////////////////
long double ratio;
if(T2-T < T-T1) ratio = (long double)(T2-T)/(long double)(T-T1);
else ratio = (long double)(T-T1)/(long double)(T2-T);
////////////////////////////////////////////////////
if(ratio < 0.4) //not a T wave
{
Twave.push_back(0);
Twave.push_back(0);
Twave.push_back(0);
twave = false;
}
else
{
//adjust center of T wave
//smooth it with gaussian, find max ?
//cwt.ToTxt(L"debugS.txt",mass+annPos+add,size);
int Tcntr = findTmax(mass+T1, T2-T1);
if(Tcntr != -1)
{
Tcntr += T1;
if(abs((Tcntr-T1)-((T2-T1)/2))<abs((T-T1)-((T2-T1)/2))) //which is close to center 0-cross or T max
T = Tcntr;
}
Twaves++;
Twave.push_back(T1);
Twave.push_back(T);
Twave.push_back(T2);
}
}
else //no T wave??? empty (QRS) <-------> (QRS)
{
Twave.push_back(0);
Twave.push_back(0);
Twave.push_back(0);
}
T=-1;
///////////////search for TWAVE///////////////////////////////////////////////////////////
cwt.CloseCWT();
///////////////search for PWAVE///////////////////////////////////////////////////////////
size = ann[n*2+2][0]-ann[n*2+1][0]; //n size of (QRS) <----> (QRS)
if(sr*ahdr.maxPQ < size)
size = sr*ahdr.maxPQ;
if(twave)
{
if(T2 > ann[n*2+2][0]-size - int(0.04*sr) ) // pwave wnd far from Twave at least on 0.02sec
size -= T2 - (ann[n*2+2][0]-size - int(0.04*sr));
}
int size23 = (ann[n*2+2][0]-ann[n*2+1][0])-size;
//size -= 0.02*sr; //impresize QRS begin detection
if(size <= 0.03*sr)
{
Pwave.push_back(0);
Pwave.push_back(0);
Pwave.push_back(0);
continue;
}
//avg = Mean(mass+annPos+size23,size); //avrg extension on boundaries
//lvl = mass[annPos+size23];
//rvl = mass[annPos+size23+size-1];
cwt.InitCWT(size,6,0,sr); //6-Gauss1 wlet
buff = cwt.CWTrans(mass+annPos+size23,ahdr.pFreq);//,false,lvl,rvl); //9Hz transform buff = size-2/3size
//cwt.ToTxt(L"debugS.txt",mass+annPos+size23,size);
//cwt.ToTxt(L"debugC.txt",buff,size);
MinMax(buff,size,min,max);
for(int i=0; i<size; i++)
{
if(buff[i]==min) P1 = i+ annPos + size23;
if(buff[i]==max) P2 = i+ annPos + size23;
}
if(P1>P2)swap(P1,P2);
//additional constraints on [P1 P P2] duration, symmetry, PQ interval
pwave = false;
if((buff[P1-annPos-size23]<0 && buff[P2-annPos-size23]>0) || (buff[P1-annPos-size23]>0 && buff[P2-annPos-size23]<0))
pwave = true;
if(pwave)
{
if((long double)(P2-P1)>=0.03*sr && (long double)(P2-P1)<=0.15*sr) //check for P wave duration 9Hz0.03 5Hz0.05
{
pwave = true; //PQ interval = [0.07 - 0.12,0.20]
if((long double)(ann[n*2+2][0]-P1)>=ahdr.minPQ*sr && (long double)(ann[n*2+2][0]-P1)<=ahdr.maxPQ*sr)
pwave = true;
else
pwave = false;
}
else
pwave = false;
}
if(pwave)
{
if(buff[P1-annPos-size23] > 0) sign = true;
else sign = false;
for(int i=P1-annPos-size23; i<P2-annPos-size23; i++)
{
if(sign)
{ if(buff[i] > 0) continue; }
else
{ if(buff[i] < 0) continue; }
P = i+ annPos + size23;
break;
}
//check for T wave symetry//////////////////////////
long double ratio;
if(P2-P < P-P1) ratio = (long double)(P2-P)/(long double)(P-P1);
else ratio = (long double)(P-P1)/(long double)(P2-P);
////////////////////////////////////////////////////
if(ratio < 0.4f) //not a P wave
{
Pwave.push_back(0);
Pwave.push_back(0);
Pwave.push_back(0);
}
else
{
Pwaves++;
Pwave.push_back(P1);
Pwave.push_back(P);
Pwave.push_back(P2);
}
}
else
{
Pwave.push_back(0);
Pwave.push_back(0);
Pwave.push_back(0);
}
P1=-1;P=-1;P2=-1;
///////////////search for PWAVE///////////////////////////////////////////////////////////
cwt.CloseCWT();
}
/////////////////get q,r,s peaks//////////////////////////////////////////////////////////
// on a denoised signal
int peaksnum=0;
int Q,R,S;
vector <int> qrsPeaks; //q,r,s peaks [ q , r , s ]
// [ 0, R , 0 ] zero if not defined
vector <char> qrsTypes; //[q,r,s] or [_,R,s], etc...
for(int n=0; n<qrsnum; n++) //fill with zeros
{
for(int i=0; i<3; i++)
{
qrsPeaks.push_back(0);
qrsTypes.push_back(' ');
}
}
buff = (long double *)malloc(len*sizeof(long double));
for(int i=0; i<len; i++)
buff[i] = mass[i];
ECGDenoise enoise;
enoise.InitDenoise(fltdir,buff,len,sr);
if(enoise.LFDenoise())
{
//ToTxt(L"lf.txt",buff,len);
long double *pbuff;
for(int n=0; n<qrsnum; n++)
{
annPos = ann[n*2][0]; //PQ
size = ann[n*2+1][0]-ann[n*2][0] + 1; //PQ-Jpnt, including Jpnt
pbuff = &buff[annPos];
Q = -1;
findRS(pbuff,size,R,S,ahdr.minUmV);
if(R != -1) R += annPos;
if(S != -1) S += annPos;
if(R!=-1 && S!=-1) // Rpeak > 0mV Speak < 0mV
{
if(S<R) //check for S
{
if(buff[R]>-buff[S])
{
Q = S;
S = -1;
size = ann[n*2+1][0]-R + 1; //including Jpnt
pbuff = &buff[R];
S = finds(pbuff,size,0.05);
if(S != -1) S += R;
}
}
else //check for Q
{
size = R-annPos + 1; //including R peak
Q = findq(pbuff,size,0.05);
if(Q != -1) Q += annPos;
}
}
//else if only S
else if(S!=-1) //find small r if only S detected in rS large T lead
{
size = S-annPos + 1; //including S peak
pbuff = &buff[annPos];
R = findr(pbuff,size,0.05);
if(R != -1) R += annPos;
}
//else if only R
else if(R!=-1) //only R find small q,s
{
size = R-annPos + 1; //including R peak
Q = findq(pbuff,size,0.05);
if(Q != -1) Q += annPos;
size = ann[n*2+1][0]-R + 1; //including Jpnt
pbuff = &buff[R];
S = finds(pbuff,size,0.05);
if(S != -1) S += R;
}
//put peaks to qrsPeaks vector
if(R==-1 && S==-1) //no peaks
{
ann[n*2][1] = 16; //ARTEFACT
//remove P,T
if(n != 0)
{
if(Pwave[3*(n-1)])
{
Pwaves--;
Pwave[3*(n-1)] = 0;
Pwave[3*(n-1)+1] = 0;
Pwave[3*(n-1)+2] = 0;
}
}
if(n != qrsnum-1)
{
if(Twave[3*n])
{
Twaves--;
Twave[3*n] = 0;
Twave[3*n+1] = 0;
Twave[3*n+2] = 0;
}
}
}
if(Q != -1)
{
peaksnum++;
qrsPeaks[n*3] = Q;
if(fabsl(buff[Q]) > 0.5)
qrsTypes[n*3] = 17; //'Q';
else
qrsTypes[n*3] = 15; //'q';
}
if(R != -1)
{
peaksnum++;
qrsPeaks[n*3+1] = R;
if(fabsl(buff[R]) > 0.5)
qrsTypes[n*3+1] = 48; //'R';
else
qrsTypes[n*3+1] = 47; //'r';
}
if(S != -1)
{
peaksnum++;
qrsPeaks[n*3+2] = S;
if(fabsl(buff[S]) > 0.5)
qrsTypes[n*3+2] = 50; //'S';
else
qrsTypes[n*3+2] = 49; //'s';
}
}
}
free(buff);
/////////////////get q,r,s peaks//////////////////////////////////////////////////////////
///////////////////////// complete annotation array///////////////////////////////////////
maNum = 0;
//Pwave vec size = Twave vec size
annNum = Pwaves*3 + qrsnum*2+peaksnum + Twaves*3 + (int)MA.size(); //P1 P P2 [QRS] T1 T T2 noise annotation
if(annNum > qrsnum) //42-(p 43-p) 24-Pwave
{ //44-(t 45-t) 27-Twave
ANN = new int*[annNum]; // [samps] [type] [?aux data]
for(int i=0; i<annNum; i++)
ANN[i] = new int[3];
int index=0; //index to ANN
int qindex=0; //index to qrsANN
for(int i=0; i<(int)Twave.size(); i+=3) //Twave=Pwaves=qrsPeaks size
{
//QRS complex
ANN[index][0] = ann[qindex][0]; //(QRS
ANN[index][1] = ann[qindex][1]; //
ANN[index++][2] = ann[qindex++][2]; //aux
if(qrsPeaks[i]) //q
{
ANN[index][0] = qrsPeaks[i];
ANN[index][1] = qrsTypes[i];
ANN[index++][2] = -1; //no aux
}
if(qrsPeaks[i+1]) //r
{
ANN[index][0] = qrsPeaks[i+1];
ANN[index][1] = qrsTypes[i+1];
ANN[index++][2] = -1; //no aux
}
if(qrsPeaks[i+2]) //s
{
ANN[index][0] = qrsPeaks[i+2];
ANN[index][1] = qrsTypes[i+2];
ANN[index++][2] = -1; //no aux
}
ANN[index][0] = ann[qindex][0]; //QRS)
ANN[index][1] = ann[qindex][1]; //
ANN[index++][2] = ann[qindex++][2]; //aux
//T wave
if(Twave[i])
{
ANN[index][0] = Twave[i]; //(t
ANN[index][1] = 44;
ANN[index++][2] = -1; //no aux
ANN[index][0] = Twave[i+1]; //T
ANN[index][1] = 27;
ANN[index++][2] = -1; //no aux
ANN[index][0] = Twave[i+2]; //t)
ANN[index][1] = 45;
ANN[index++][2] = -1; //no aux
}
//P wave
if(Pwave[i])
{
ANN[index][0] = Pwave[i]; //(t
ANN[index][1] = 42;
ANN[index++][2] = -1; //no aux
ANN[index][0] = Pwave[i+1]; //T
ANN[index][1] = 24;
ANN[index++][2] = -1; //no aux
ANN[index][0] = Pwave[i+2]; //t)
ANN[index][1] = 43;
ANN[index++][2] = -1; //no aux
}
if(!Twave[i] && !Pwave[i]) //check for MA noise
{
for(int m=maNum; m<(int)MA.size(); m++)
{
if(MA[m]>ann[qindex-1][0] && MA[m]<ann[qindex][0])
{
ANN[index][0] = MA[m]; //Noise
ANN[index][1] = 14;
ANN[index++][2] = -1; //no aux
maNum++;
break;
}
}
}
}
//last QRS complex
int ii=3*(qrsnum-1);
ANN[index][0] = ann[qindex][0]; //(QRS
ANN[index][1] = ann[qindex][1]; //
ANN[index++][2] = ann[qindex++][2]; //aux
if(qrsPeaks[ii]) //q
{
ANN[index][0] = qrsPeaks[ii];
ANN[index][1] = qrsTypes[ii];
ANN[index++][2] = -1; //no aux
}
if(qrsPeaks[ii+1]) //r
{
ANN[index][0] = qrsPeaks[ii+1];
ANN[index][1] = qrsTypes[ii+1];
ANN[index++][2] = -1; //no aux
}
if(qrsPeaks[ii+2]) //s
{
ANN[index][0] = qrsPeaks[ii+2];
ANN[index][1] = qrsTypes[ii+2];
ANN[index++][2] = -1; //no aux
}
ANN[index][0] = ann[qindex][0]; //QRS)
ANN[index][1] = ann[qindex][1]; //
ANN[index++][2] = ann[qindex++][2]; //aux
//check if noise after last qrs
if(maNum < (int)MA.size())
{
if(MA[maNum]>ann[qindex-1][0])
{
ANN[index][0] = MA[maNum]; //Noise
ANN[index][1] = 14;
ANN[index++][2] = -1;
}
}
return ANN;
}
else
return 0;
}
void Annotation::findRS(long double *mass, int len, int &R, int &S, long double err) //find RS or QR
{
long double tmp,min,max;
MinMax(mass,len, min,max);
R = -1; S = -1;
if(!( max<0.0 || max==mass[0] || max==mass[len-1] || max<err ))//(fabsl(max-mass[0])<err && fabsl(max-mass[len-1])<err) ))
{
for(int i=1; i<len-1; i++)
if(mass[i] == max)
{
R=i;
break;
}
}
if(!( min>0.0 || min==mass[0] || min==mass[len-1] || -min<err ))//(fabsl(min-mass[0])<err && fabsl(min-mass[len-1])<err) ))
{
for(int i=1; i<len-1; i++)
if(mass[i] == min)
{
S=i;
break;
}
}
}
int Annotation::findTmax(long double *mass, int len) //find T max/min peak position
{
long double min,max;
MinMax(mass,len, min,max);
int tmin=-1,tmax=-1;
for(int i=0; i<len; i++)
{
if(mass[i] == max)
{
tmax = i;
break;
}
}
for(int i=0; i<len; i++)
{
if(mass[i] == min)
{
tmin = i;
break;
}
}
//max closest to the center
if(tmin==-1 || tmax==-1) //??no max min found
return -1;
else
{
if(abs(tmax-(len/2))<abs(tmin-(len/2)))
return tmax;
else
return tmin;
}
}
int Annotation::findr(long double *mass, int len, long double err) //find small r in PQ-S
{
long double tmp,min,max;
MinMax(mass,len, min,max);
if(max<0.0 || max==mass[0] || max==mass[len-1] || fabsl(max-mass[0])<err) return -1;
else tmp = max;
for(int i=1; i<len-1; i++)
{
if(mass[i] == tmp)
return i;
}
return -1;
}
int Annotation::findq(long double *mass, int len, long double err) //find small q in PQ-R
{
long double tmp,min,max;
MinMax(mass,len, min,max);
if(min>0.0 || min==mass[0] || min==mass[len-1] || fabsl(min-mass[0])<err) return -1;
else tmp = min;
for(int i=1; i<len-1; i++)
{
if(mass[i] == tmp)
return i;
}
return -1;
}
int Annotation::finds(long double *mass, int len, long double err) //find small s in R-Jpnt
{
long double tmp,min,max;
MinMax(mass,len, min,max);
if(min>0.0 || min==mass[0] || min==mass[len-1] || fabsl(min-mass[len-1])<err) return -1;
else tmp = min;
for(int i=1; i<len-1; i++)
{
if(mass[i] == tmp)
return i;
}
return -1;
}
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
void Annotation::AddANN(int add)
{
for(int i=0; i<qrsNum; i++)
qrsANN[i][0] += add;
for(int i=0; i<annNum; i++)
ANN[i][0] += add;
}
//-----------------------------------------------------------------------------
// SaveANN (**aux) aux data
//-----------------------------------------------------------------------------
bool Annotation::SaveANN(wchar_t *name, int **ann, int nums)
{
int samps;
unsigned short anncode=0;
unsigned short type;
DWORD bytes;
char buff[6];
fp = CreateFileW(name,GENERIC_WRITE | GENERIC_READ,0,0,CREATE_ALWAYS,FILE_ATTRIBUTE_NORMAL,0);
samps = ann[0][0];
sprintf(buff,"%c%c%c%c%c%c",0,0xEC,HIWORD(LOBYTE(samps)),HIWORD(HIBYTE(samps)),LOWORD(LOBYTE(samps)),LOWORD(HIBYTE(samps)));
WriteFile(fp,buff,6,&bytes,0);
type = ann[0][1];
anncode |= (type<<10);
sprintf(buff,"%c%c",(char)LOBYTE(anncode),(char)HIBYTE(anncode));
WriteFile(fp,buff,2,&bytes,0);
for(int i=1; i<nums; i++)
{
samps = ann[i][0] - ann[i-1][0];
anncode = 0;
type = ann[i][1];
if(samps > 1023)
{
sprintf(buff,"%c%c%c%c%c%c",0,0xEC,HIWORD(LOBYTE(samps)),HIWORD(HIBYTE(samps)),LOWORD(LOBYTE(samps)),LOWORD(HIBYTE(samps)));
WriteFile(fp,buff,6,&bytes,0);
}
else
anncode = samps;
anncode |= (type<<10);
sprintf(buff,"%c%c",(char)LOBYTE(anncode),(char)HIBYTE(anncode));
WriteFile(fp,buff,2,&bytes,0);
}
sprintf(buff,"%c%c",0,0);
WriteFile(fp,buff,2,&bytes,0);
CloseHandle(fp);
return true;
}
bool Annotation::SaveQTseq(wchar_t *name, int **ann, int annsize, long double sr, int len)
{
vector <long double> QT;
int q=0,t=0;
for(int i=0; i<annsize; i++)
{
switch(ann[i][1])
{
case 14: //noise
case 15: //q
case 16: //artifact
case 17: //Q
case 18: //ST change
case 19: //T change
case 20: //systole
case 21: //diastole
case 22:
case 23:
case 24: //P
case 25:
case 26:
case 27:
case 28:
case 29:
case 30:
case 31:
case 32:
case 33:
case 36:
case 37:
case 39:
case 40:
case 42: //(p
case 43: //p)
case 44: //(t
case 47: //r
case 48: //R
case 49: //s
case 50: //S
continue;
}
if(ann[i][1] == 45) //45 - t)
{
t = ann[i][0];
if(q < t)
QT.push_back((long double)(t-q)/sr);
}
else
{
/*if(i+1<annsize && (ann[i+1][1]==47 || ann[i+1][1]==48)) //r only
q = ann[i+1][0];
else if(i+2<annsize && (ann[i+2][1]==47 || ann[i+2][1]==48)) //q,r
q = ann[i+2][0];
else*/
q = ann[i][0];
}
}
if(QT.size())
{
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
memcpy(hdr.hdr,"DATA",4);
hdr.size = QT.size();
hdr.sr = float( (long double)QT.size() / ((long double)len/sr) );
hdr.bits = 32;
hdr.umv = 1;
SaveFile(name, &QT[0], &hdr);
return true;
}
else
return false;
}
bool Annotation::SavePQseq(wchar_t *name, int **ann, int annsize, long double sr, int len)
{
vector <long double> PQ;
int p=len,q=0;
for(int i=0; i<annsize; i++)
{
switch(ann[i][1])
{
case 14: //noise
case 15: //q
case 16: //artifact
case 17: //Q
case 18: //ST change
case 19: //T change
case 20: //systole
case 21: //diastole
case 22:
case 23:
case 24: //P
case 25:
case 26:
case 27:
case 28:
case 29:
case 30:
case 31:
case 32:
case 33:
case 36:
case 37:
case 39:
case 40:
case 43: //p)
case 44: //(t
case 45: //t)
case 47: //r
case 48: //R
case 49: //s
case 50: //S
continue;
}
if(ann[i][1] == 42) //42 - (p
p = ann[i][0];
else
{
q = ann[i][0];
if(p < q)
{
PQ.push_back((long double)(q-p)/sr);
p=len;
}
}
}
if(PQ.size())
{
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
memcpy(hdr.hdr,"DATA",4);
hdr.size = PQ.size();
hdr.sr = float( (long double)PQ.size() / ((long double)len/sr) );
hdr.bits = 32;
hdr.umv = 1;
SaveFile(name,&PQ[0],&hdr);
return true;
}
else
return false;
}
bool Annotation::SavePPseq(wchar_t *name, int **ann, int annsize, long double sr, int len)
{
vector <long double> PP;
int p1,p2;
for(int i=0; i<annsize; i++)
{
if(ann[i][1] == 42) //42 - (p
p1 = ann[i][0];
else if(ann[i][1] == 43) //43 - p)
{
p2 = ann[i][0];
PP.push_back((long double)(p2-p1)/sr);
}
}
if(PP.size())
{
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
memcpy(hdr.hdr,"DATA",4);
hdr.size = PP.size();
hdr.sr = float( (long double)PP.size() / ((long double)len/sr) );
hdr.bits = 32;
hdr.umv = 1;
SaveFile(name,&PP[0],&hdr);
return true;
}
else
return false;
}
bool Annotation::GetRRseq(int **ann, int nums, long double sr, vector<long double> *RR, vector<int> *RRpos)
{
int add=-1;
long double rr, r1,r2;
RR->clear();
RRpos->clear();
//R peak or S peak annotation
bool rrs=true;
int rNum=0,sNum=0;
for(int i=0; i<nums; i++)
{
if(ann[i][1]==47 || ann[i][1]==48) rNum++;
else if(ann[i][1]==49 || ann[i][1]==50) sNum++;
}
if(int(1.2f*(float)rNum) < sNum)
rrs=false; //R peaks less than S ones
for(int i=0; i<nums; i++)
{
switch(ann[i][1])
{
case 0: //non beats
case 15: //q
case 17: //Q
case 18: //ST change
case 19: //T change
case 20: //systole
case 21: //diastole
case 22:
case 23:
case 24:
case 25:
case 26:
case 27:
case 28:
case 29:
case 30:
case 31:
case 32:
case 33:
case 36:
case 37:
case 40: //')' J point
case 42: //(p
case 43: //p)
case 44: //(t
case 45: //t)
case 47: //r
case 48: //R
case 49: //s
case 50: //S
continue;
case 14: //noise
case 16: //artifact
add = -1;
continue;
}
if(add != -1)
{
//annotation on RRs peaks
if(rrs)
{
if(i+1<nums && (ann[i+1][1]==47 || ann[i+1][1]==48)) //r only
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==47 || ann[i+2][1]==48)) //q,r
r2 = ann[i+2][0];
else //(ann[i][1]==N,ECT,...) //no detected R only S
r2 = ann[i][0];
if(add+1<nums && (ann[add+1][1]==47 || ann[add+1][1]==48))
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==47 || ann[add+2][1]==48))
r1 = ann[add+2][0];
else //(ann[add][1]==N,ECT,...) //no detected R only S
r1 = ann[add][0];
}
//annotation on S peaks
else
{
if(i+1<nums && (ann[i+1][1]==40)) //N)
r2 = ann[i][0];
else if(i+1<nums && (ann[i+1][1]==49 || ann[i+1][1]==50)) //Sr
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==49 || ann[i+2][1]==50)) //rS
r2 = ann[i+2][0];
else if(i+3<nums && (ann[i+3][1]==49 || ann[i+3][1]==50)) //errQ rS
r2 = ann[i+3][0];
else if(i+1<nums && (ann[i+1][1]==47 || ann[i+1][1]==48)) //no S
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==47 || ann[i+2][1]==48)) //no S
r2 = ann[i+2][0];
if(add+1<nums && (ann[add+1][1]==40)) //N)
r1 = ann[add][0];
else if(add+1<nums && (ann[add+1][1]==49 || ann[add+1][1]==50))
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==49 || ann[add+2][1]==50))
r1 = ann[add+2][0];
else if(add+3<nums && (ann[add+3][1]==49 || ann[add+3][1]==50))
r1 = ann[add+3][0];
else if(add+1<nums && (ann[add+1][1]==47 || ann[add+1][1]==48)) //no S
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==47 || ann[add+2][1]==48)) //no S
r1 = ann[add+2][0];
}
rr = 60.0 / ((r2-r1)/sr);
if(rr >= ahdr.minbpm && rr <= ahdr.maxbpm)
{
RR->push_back( rr ); //in bpm
RRpos->push_back( r1 );
}
}
add = i;
}
if(RR->size())
return true;
else
return false;
}
bool Annotation::SaveRRseq(wchar_t *name, int **ann, int nums, long double sr, int len)
{
vector <long double> RR;
int add=-1;
long double rr, r1,r2;
//R peak or S peak annotation////////////////////////////
bool rrs=true;
int rNum=0,sNum=0;
for(int i=0; i<nums; i++)
{
if(ann[i][1]==47 || ann[i][1]==48) rNum++;
else if(ann[i][1]==49 || ann[i][1]==50) sNum++;
}
if(int(1.1f*(float)rNum) < sNum)
{
rrs=false; //R peaks less than S ones
wcscat(name, L"_SS.dat");
}
else
wcscat(name, L"_RR.dat");
////////////////////////////////////////////////////////
for(int i=0; i<nums; i++)
{
switch(ann[i][1])
{
case 0: //non beats
case 15: //q
case 17: //Q
case 18: //ST change
case 19: //T change
case 20: //systole
case 21: //diastole
case 22:
case 23:
case 24:
case 25:
case 26:
case 27:
case 28:
case 29:
case 30:
case 31:
case 32:
case 33:
case 36:
case 37:
case 40: //')' J point
case 42: //(p
case 43: //p)
case 44: //(t
case 45: //t)
case 47: //r
case 48: //R
case 49: //s
case 50: //S
continue;
case 14: //noise
case 16: //artifact
add = -1;
continue;
}
if(add != -1)
{
//annotation on RRs peaks
if(rrs)
{
if(i+1<nums && (ann[i+1][1]==47 || ann[i+1][1]==48)) //r only
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==47 || ann[i+2][1]==48)) //q,r
r2 = ann[i+2][0];
else //(ann[i][1]==N,ECT,...) //no detected R only S
r2 = ann[i][0];
if(add+1<nums && (ann[add+1][1]==47 || ann[add+1][1]==48))
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==47 || ann[add+2][1]==48))
r1 = ann[add+2][0];
else //(ann[add][1]==N,ECT,...) //no detected R only S
r1 = ann[add][0];
}
//annotation on S peaks
else
{
if(i+1<nums && (ann[i+1][1]==40)) //N)
r2 = ann[i][0];
else if(i+1<nums && (ann[i+1][1]==49 || ann[i+1][1]==50)) //Sr
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==49 || ann[i+2][1]==50)) //rS
r2 = ann[i+2][0];
else if(i+3<nums && (ann[i+3][1]==49 || ann[i+3][1]==50)) //errQ rS
r2 = ann[i+3][0];
else if(i+1<nums && (ann[i+1][1]==47 || ann[i+1][1]==48)) //no S
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==47 || ann[i+2][1]==48)) //no S
r2 = ann[i+2][0];
if(add+1<nums && (ann[add+1][1]==40)) //N)
r1 = ann[add][0];
else if(add+1<nums && (ann[add+1][1]==49 || ann[add+1][1]==50))
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==49 || ann[add+2][1]==50))
r1 = ann[add+2][0];
else if(add+3<nums && (ann[add+3][1]==49 || ann[add+3][1]==50))
r1 = ann[add+3][0];
else if(add+1<nums && (ann[add+1][1]==47 || ann[add+1][1]==48)) //no S
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==47 || ann[add+2][1]==48)) //no S
r1 = ann[add+2][0];
}
rr = 60.0 / ((r2-r1)/sr);
if(rr >= ahdr.minbpm && rr <= ahdr.maxbpm)
RR.push_back( rr ); //in bpm
}
add = i;
}
if(RR.size())
{
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
memcpy(hdr.hdr,"DATA",4);
hdr.size = RR.size();
hdr.sr = float( (long double)RR.size() / ((long double)len/sr) );
hdr.bits = 32;
hdr.umv = 1;
SaveFile(name,&RR[0],&hdr);
return true;
}
else
return false;
}
bool Annotation::SaveRRnseq(wchar_t *name, int **ann, int nums, long double sr, int len)
{
vector <long double> RR;
int add=-1;
long double rr,r1,r2;
//R peak or S peak annotation////////////////////////////
bool rrs=true;
int rNum=0,sNum=0;
for(int i=0; i<nums; i++)
{
if(ann[i][1]==47 || ann[i][1]==48) rNum++;
else if(ann[i][1]==49 || ann[i][1]==50) sNum++;
}
if(int(1.1f*(float)rNum) < sNum)
{
rrs=false; //R peaks less than S ones
wcscat(name, L"_SSn.dat");
}
else
wcscat(name, L"_RRn.dat");
////////////////////////////////////////////////////////
for(int i=0; i<nums; i++)
{
switch(ann[i][1])
{
case 1: //N
if(add == -1)
{
add = i;
continue;
}
break;
case 0: //ectopic beats
case 4:
case 5:
case 6:
case 7:
case 8:
case 9:
case 10:
case 11:
case 12:
case 13:
case 14: //noise
case 16: //artefact
case 34:
case 35:
case 38:
case 46:
add = -1; //reset counter
continue;
default: //other types
continue;
}
//annotation on RRs peaks
if(rrs)
{
if(i+1<nums && (ann[i+1][1]==47 || ann[i+1][1]==48)) //r only
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==47 || ann[i+2][1]==48)) //q,r
r2 = ann[i+2][0];
else //(ann[i][1]==N,ECT,...) //no detected R only S
r2 = ann[i][0];
if(add+1<nums && (ann[add+1][1]==47 || ann[add+1][1]==48))
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==47 || ann[add+2][1]==48))
r1 = ann[add+2][0];
else //(ann[add][1]==N,ECT,...) //no detected R only S
r1 = ann[add][0];
}
//annotation on S peaks
else
{
if(i+1<nums && (ann[i+1][1]==40)) //N)
r2 = ann[i][0];
else if(i+1<nums && (ann[i+1][1]==49 || ann[i+1][1]==50)) //Sr
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==49 || ann[i+2][1]==50)) //rS
r2 = ann[i+2][0];
else if(i+3<nums && (ann[i+3][1]==49 || ann[i+3][1]==50)) //errQ rS
r2 = ann[i+3][0];
else if(i+1<nums && (ann[i+1][1]==47 || ann[i+1][1]==48)) //no S
r2 = ann[i+1][0];
else if(i+2<nums && (ann[i+2][1]==47 || ann[i+2][1]==48)) //no S
r2 = ann[i+2][0];
if(add+1<nums && (ann[add+1][1]==40)) //N)
r1 = ann[add][0];
else if(add+1<nums && (ann[add+1][1]==49 || ann[add+1][1]==50))
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==49 || ann[add+2][1]==50))
r1 = ann[add+2][0];
else if(add+3<nums && (ann[add+3][1]==49 || ann[add+3][1]==50))
r1 = ann[add+3][0];
else if(add+1<nums && (ann[add+1][1]==47 || ann[add+1][1]==48)) //no S
r1 = ann[add+1][0];
else if(add+2<nums && (ann[add+2][1]==47 || ann[add+2][1]==48)) //no S
r1 = ann[add+2][0];
}
rr = 60.0 / ((r2-r1)/sr);
if(rr >= ahdr.minbpm && rr <= ahdr.maxbpm)
RR.push_back( rr ); //in bpm
add = i;
}
if(RR.size())
{
DATAHDR hdr;
memset(&hdr,0,sizeof(DATAHDR));
memcpy(hdr.hdr,"DATA",4);
hdr.size = RR.size();
hdr.sr = float( (long double)RR.size() / ((long double)len/sr) );
hdr.bits = 32;
hdr.umv = 1;
SaveFile(name,&RR[0],&hdr);
return true;
}
else
return false;
}
//////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////CLASS DENOISE///////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
ECGDenoise::ECGDenoise(): tempBuff(0)
{
}
ECGDenoise::~ECGDenoise()
{
if(tempBuff) free(tempBuff);
}
//-----------------------------------------------------------------------------
void ECGDenoise::InitDenoise(wchar_t *fltdir, long double *data, int size, long double sr, bool mirror)
{
wcscpy(this->fltdir, fltdir);
ecgBuff = data;
SR = sr;
Len = size;
if(tempBuff) free(tempBuff);
tempBuff = (long double *)malloc(sizeof(long double)*(Len+2*SR)); // [SR add] [sig] [SR add]
for(int i=0; i<Len; i++) //signal
tempBuff[i+(int)SR] = ecgBuff[i];
if(Len < SR) mirror = false;
if(mirror) //mirror extension of signal
{
for(int i=0; i<(int)SR; i++) //left side
tempBuff[i] = ecgBuff[(int)SR-i];
for(int i=Len+SR; i<Len+2*SR; i++) //right side
tempBuff[i] = ecgBuff[(Len-2) - (i-(Len+(int)SR))];
}
else
{
for(int i=0; i<(int)SR; i++) //left side
tempBuff[i] = ecgBuff[0];
for(int i=Len+SR; i<Len+2*SR; i++) //right side
tempBuff[i] = ecgBuff[Len-1];
}
}
void ECGDenoise::CloseDenoise()
{
if(tempBuff)
{
free(tempBuff);
tempBuff = 0;
}
}
//-----------------------------------------------------------------------------
bool ECGDenoise::LFDenoise()
{
wchar_t filter[256];
//get base line J///////
int J = ceill(log2(SR/0.8)) - 1;
wcscpy(filter,fltdir);
wcscat(filter,L"\\daub2.flt");
if(InitFWT(filter,tempBuff,Len+2*SR) == false)
return false;
//transform////////////////////////
fwtTrans(J);
///////////////////////////////////
int *Jnumbs = GetJnumbs(J,Len+2*SR);
long double *lo = getSpec();
for(int i=0; i<Jnumbs[0]; i++)
lo[i] = 0.0;
//synth/////////////////////////////
fwtSynth(J);
////////////////////////////////////
for(int i=0; i<Len; i++)
ecgBuff[i] = lo[i+(int)SR];
CloseFWT();
return true;
}
bool ECGDenoise::HFDenoise()
{
wchar_t filter[256];
//get HF scale J///////
int J = ceill(log2(SR/23.0)) - 2; //[30Hz - ...] hf denoising
wcscpy(filter,fltdir);
wcscat(filter,L"\\bior97.flt");
if(InitFWT(filter,tempBuff,Len+2*SR) == false)
return false;
//transform////////////////////////
fwtTrans(J);
///////////////////////////////////
int *Jnumbs = GetJnumbs(J,Len+2*SR);
int hinum,lonum;
HiLoNumbs(J,Len+2*SR,hinum,lonum);
long double *lo = getSpec();
long double *hi = getSpec() + (int(Len+2*SR)-hinum);
int window; //3sec window
for(int j=J; j>0; j--)
{
window = 3.0*SR/powl(2.0,(long double)j);
Denoise(hi,Jnumbs[J-j],window);
hi += Jnumbs[J-j];
}
//synth/////////////////////////////
fwtSynth(J);
////////////////////////////////////
for(int i=0; i<Len; i++)
ecgBuff[i] = lo[i+(int)SR];
CloseFWT();
return true;
}
bool ECGDenoise::LFHFDenoise()
{
wchar_t filter[256];
int J,*Jnumbs;
long double *lo,*hi;
//get base line J///////
J = ceill(log2(SR/0.8)) - 1;
wcscpy(filter,fltdir);
wcscat(filter,L"\\daub2.flt");
if(InitFWT(filter,tempBuff,Len+2*SR) == false)
return false;
//transform////////////////////////
fwtTrans(J);
///////////////////////////////////
Jnumbs = GetJnumbs(J,Len+2*SR);
lo = getSpec();
for(int i=0; i<Jnumbs[0]; i++)
lo[i] = 0.0;
//synth/////////////////////////////
fwtSynth(J);
////////////////////////////////////
for(int i=0; i<Len+2*SR; i++)
tempBuff[i] = lo[i];
////////get min max///////////////////
long double min,max;
MinMax(&lo[(int)SR],Len,min,max);
CloseFWT();
//get HF scale J///////
J = ceill(log2(SR/23.0)) - 2; //[30Hz - ...] hf denoising
wcscpy(filter,fltdir);
wcscat(filter,L"\\bior97.flt");
if(InitFWT(filter,tempBuff,Len+2*SR) == false)
return false;
//transform////////////////////////
fwtTrans(J);
///////////////////////////////////
Jnumbs = GetJnumbs(J,Len+2*SR);
int hinum,lonum;
HiLoNumbs(J,Len+2*SR,hinum,lonum);
lo = getSpec();
hi = getSpec() + (int(Len+2*SR)-hinum);
int window; //3sec window
for(int j=J; j>0; j--)
{
window = 3.0*SR/powl(2.0,(long double)j);
Denoise(hi,Jnumbs[J-j],window);
hi += Jnumbs[J-j];
}
//synth/////////////////////////////
fwtSynth(J);
////////////////////////////////////
for(int i=0; i<Len; i++)
ecgBuff[i] = lo[i+(int)SR];
//renormalize
nMinMax(ecgBuff,Len,min,max);
CloseFWT();
return true;
}