Click here to Skip to main content
Click here to Skip to main content
 
Add your own
alternative version

Tagged as

ECG recording, storing, filtering and recognition

, 14 Apr 2004
Full open code project for making driver and application software for ECG medical measurements.
ecg_dsp_src.zip
drvECG_1
drvECG_1.clw
drvECG_1.def
drvECG_1.dsp
drvECG_1.dsw
drvECG_1.exp
drvECG_1.ilk
drvECG_1.lib
res
Ecg_1
ECG_1.clw
ECG_1.dsp
ECG_1.dsw
ECG_DRAW.obj
ECG_Statistic_View.lib
ECG_VIEW.exp
ECG_VIEW.ilk
ECG_VIEW.lib
TestOCX.ocx
drvECG_1.exp
drvECG_1.ilk
drvECG_1.lib
resource.hm
res
ECG_1.ico
bitmap1.bmp
bmp00001.bmp
setup.bmp
ECG_Statistic_View
ECG_Statistic_View.clw
ECG_Statistic_View.def
ECG_Statistic_View.dsp
ECG_Statistic_View.dsw
res
ECG_VIEW
ECG_DRAW.obj
ECG_VIEW.APS
ECG_VIEW.clw
ECG_VIEW.def
ECG_VIEW.dll
ECG_VIEW.dsp
ECG_VIEW.dsw
ECG_VIEW.exp
ECG_VIEW.lib
res
testrecord2.zip
georgi.ecg
test_file1.zip
test_file1.ecg
// DSP_Filter.cpp: implementation of the DSP_Filter class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ECG_1.h"
#include "DSP_Filter.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

DSP_Filter::DSP_Filter()
{
	PI = 3.141592;

}

DSP_Filter::~DSP_Filter()
{

}

void DSP_Filter::RFilter_Low(double *dest, int lenght, int fsmpl, double fk)
{
	//Low filter 1 stage
	int i=0;
	int lenght_1 = lenght + 100;
	//resize X and Y arrays
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	//Move array from dest in to X array
	MoveArr(X,dest,lenght);


	//Calculate filter coefficients recursive form
	double x = fk;
	double A0 = pow((1-x),4);
	double B1 = 4*x;
	double B2 = -6*x*x;
	double B3 = 4*x*x*x;
	double B4 = -x*x*x*x;
	//i rins from 2!!!
	for(i=4;i<lenght_1;i++)
	{
		*(Y+i) = (*(X+i)*A0 + *(Y+i-1)*B1 + *(Y+i-2)*B2
			+ *(Y+i-3)*B3 + *(Y+i-4)*B4)*1.005;
	}
	MoveArr(dest,Y,lenght);
	//Free the used memory
	Delete_Arrs();
}

void DSP_Filter::ZeroArrs(int lenght)
{
	for(int i = 0;i<lenght;i++)
	{
		*(Y+i) = 0;//Zero the filter array
		*(X+i) = 0;//Zero the remaine X array
	}

}

void DSP_Filter::MoveArr(double *dest, double *source, int lenght)
{
	//Move the array in to he *dest position
	for(int i=0;i<lenght;i++)
	{	
		*(dest+i) = *(source+i);
	}

}

void DSP_Filter::Message(double var)
{
	CString str;
	str.Format("%lf",var);
	MessageBox(NULL,str,"Message",MB_OK);

}

void DSP_Filter::RFilter_Noht(double *dest, int lenght, int fsmpl, double fcut,double bw)
{
	//Do not work good at cut of frequency =<0.1 fsmpl!!!
	int lenght_1 = lenght+100;
	Greate_Arrs(lenght_1);
	//Zerro the new maked arrays
	ZeroArrs(lenght_1);
	MoveArr(X,dest,lenght);
	//Calculate noht koefizients
	double K,R;
	fcut = fcut/double(fsmpl);
	R = 1-3*bw;
	K = (1-2*R*cos(2*PI*fcut)+R*R)/(2-2*cos(2*PI*fcut));
	
	double A0 = K;
	double A1 = -2*K*cos(2*PI*fcut);
	double A2 = K;
	double B1 = 2*R*cos(2*PI*fcut);
	double B2 = -R*R;
	//i rins from 2!!!
	for(int i=2;i<lenght_1;i++)
	{
		*(Y+i) = (*(X+i)*A0 + *(X+i-1)*A1 + *(X+i-2)*A2
			+ *(Y+i-1)*B1 + *(Y+i-2)*B2);
	}
	MoveArr(dest,Y,lenght);
	//Free memory
	Delete_Arrs();
}

void DSP_Filter::WFilter_High(double *dest, double *source, int lenght,double amp, int fsmpl, double fk)
{
	//Win=-sync filter hight > 0.67hz
	int lenght_1 = lenght+400;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,dest,lenght);
	int M = 200;
	int corect=0;
	if(fsmpl==250)
		corect=0;
	if(fsmpl==500)
		corect=10;

	SlideArrRight(Y,X,lenght_1,(M/2));
	MoveArr(X,Y,lenght_1);
	for(int i=0;i<lenght_1;i++)
	{
		*(Y+i)=0;
	}
	for(int j=M;j<lenght_1-M;j++)
	{
		for(i=0;i<M;i++)
		{
			*(Y+j)=*(Y+j)+*(X+j-i)*H[i];
		//	*(Y+j) = *(Y+j)+*(X+j+i)*H[i];
		}
	}
	SlideArrLeft(Y,Y,lenght_1,M-corect);//
//	FilterLowFromTo(Y,Y,0,M,lenght);
	MoveArr(dest,Y,lenght);
	//Free memory
	Delete_Arrs();

}

void DSP_Filter::InvertArr(double *dest, double *source, int lenght)
{
	//Invert the array
	//Array source first nust be normalized to 0+-
	int lenght_1 = lenght+100;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,dest,lenght);
	for(int i=0;i<lenght_1;i++)
	{
		*(Y+i) = *(X+i)*-1; 

	}
	MoveArr(dest,Y,lenght);
	Delete_Arrs();
}

void DSP_Filter::Delta1(double *dest, double *source, int lenght)
{
	int lenght_1 = lenght+100;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,source,lenght);

	//finding the destination array 1delta
	for(int i=0;i<lenght;i++)
	{
		*(Y+i) = *(X+i+1) - *(X+i);
	}
	MoveArr(dest,Y,lenght);
	Delete_Arrs();
}

void DSP_Filter::Greate_Arrs(int new_lenght)
{
	//Creates new arr X&Y with lenght = mew_lenght
	b_arr_flag = 1;
	X = new double [new_lenght];
	Y = new double [new_lenght];
}

void DSP_Filter::Delete_Arrs()
{
	//deletes the created arrs
	if(b_arr_flag == 1)
	{
		delete X;
		delete Y;
		b_arr_flag = 0;
	}
}

double DSP_Filter::Mean(double *source, int lenght)
{
	//finding the mean
	double MEAN=0;
	for(int i=0;i<lenght;i++)
	{
		MEAN = MEAN + *(source+i);
	}
	MEAN = MEAN/lenght;
//	CString str;
//	str.Format("%lf",MEAN);
//	MessageBox(str,"Mean",MB_OK);
	return MEAN;
}

void DSP_Filter::NormalizeArr(double *ar, double mean, int lenght)
{
	double temp = 0;
	for(int i=0;i<lenght;i++)
	{
		if(*(ar+i)>mean)
			temp = *(ar+i)-mean;
		if(*(ar+i)==mean)
			temp = 0;
		if(*(ar+i)<mean)
			temp = *(ar+i)-mean;
		*(ar+i) = temp;
	}
}

void DSP_Filter::SetArrToNormalize(double *dest,int lenght)
{
	NormalizeArr(dest,Mean(dest,lenght),lenght);

}

void DSP_Filter::DeltaReverse1(double *dest, double *source, int lenght)
{
	int lenght_1 = lenght+100;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,source,lenght);

	//finding the destination array reverse delta
	for(int i=0;i<lenght;i++)
	{
		*(Y+i+1) = *(Y+i) + *(X+i+1);
	}
	MoveArr(dest,Y,lenght);
	Delete_Arrs();

}

void DSP_Filter::Delta2(double *dest, double *source, int lenght)
{
	//Finding delta2
	int lenght_1 = lenght+100;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,source,lenght);
	
	//finding the destination array 1delta
	for(int i=0;i<lenght;i++)
	{
		*(Y+i) = *(X+i+2) - 2*(*(X+i+1)) + *(X+i);
	}
	MoveArr(dest,Y,lenght);
	Delete_Arrs();

}

void DSP_Filter::DeltaReverse2(double *dest, double *source, int lenght)
{
	//Finding reverce delta2
	DeltaReverse1(dest,source,lenght);
	DeltaReverse1(dest,source,lenght);
}

double DSP_Filter::StandardDeviation(double *source, int lenght)
{
	double dev = 0;
	double mean = Mean(source,lenght);
	for(int i=0;i<lenght;i++)
	{
		dev = dev + (*(source+i)-mean)*(*(source+i)-mean);
	}
	dev = dev/(lenght-1);
	dev = sqrt(dev);
	return dev;

}

void DSP_Filter::SlideArrLeft(double *dest, double *source, int lenght, int slide)
{
	for(int i=0;i<lenght-slide;i++)
	{
		*(dest+i) = *(source+i+slide);
	}
}

void DSP_Filter::SlideArrRight(double *dest, double *source, int lenght, int slide)
{
	for(int i=0;i<lenght-slide;i++)
	{
		*(dest+i+slide) = *(source+i);
	}

}

int DSP_Filter::CreateWF50Noht(int fsmpl,int f_kernel)
{
	int M = f_kernel;
	double k=500/fsmpl;
	double sum = 0;

	//calculate band noht 50Hz first <45Hz
	double FC = 0.09*k;
	for(int i=0;i<M;i++)
	{
		if((i-M/2)==0)
			B1[i] = 2*PI*FC;
		if(0>(i-M/2)|0<(i-M/2))
			B1[i] = sin(2*PI*FC*(i-M/2))/(i-M/2);
		B1[i] = B1[i]*(0.54-0.46*cos(2*PI*i/M));
	}
	//Calculate unity gain at dc
	sum = 0;
	for(i=0;i<M;i++)
	{
		sum = sum+B1[i];
	}
	//Normalize for unity gain
	for(i=0;i<M;i++)
	{
		B1[i] = B1[i]/sum;
	}
	//Calculate >55Hz
	FC = 0.11*k;
	for(i=0;i<M;i++)
	{
		if((i-M/2)==0)
			B2[i] = 2*PI*FC;
		if(0>(i-M/2)|0<(i-M/2))
			B2[i] = sin(2*PI*FC*(i-M/2))/(i-M/2);
		B2[i] = B2[i]*(0.54-0.46*cos(2*PI*i/M));
	}
	sum = 0;
	for(i=0;i<M;i++)
	{
		sum = sum+B2[i];
	}
	for(i=0;i<M;i++)
	{
		B2[i] = B2[i]/sum;
	}
	for(i=0;i<M;i++)
	{
		B2[i] = -B2[i];
	}
	B2[M/2]=B2[M/2]+1;

	//Calculate summ kernel
	for(i=0;i<M;i++)
	{
		B1[i] = B1[i]+B2[i];
	}

	return M;
}

int DSP_Filter::CreateWFHigh(int fsmpl, int f_kernel)
{
	int M = f_kernel;
	double k=500/fsmpl;
	double FC = 0.00134*k;
	//Filter > 0.67Hz
	for(int i=0;i<M;i++)
	{
		if((i-M/2)==0)
			H[i] = 2*PI*FC;
		if(0>(i-M/2)|0<(i-M/2))
			H[i] = sin(2*PI*FC*(i-M/2))/(i-M/2);
		H[i] = H[i]*(0.54-0.46*cos(2*PI*i/M));
	}
	double sum=0;
	for(i=0;i<M;i++)
	{
		sum = sum+H[i];
	}
	for(i=0;i<M;i++)
	{
		H[i] = H[i]/sum;
	}
	for(i=0;i<M;i++)
	{
		H[i] = -H[i];
	}
	H[M/2]=H[M/2]+1;
	return M;
}

int DSP_Filter::CreateWFLow(int fsmpl, int f_kernel)
{
	int M = f_kernel;
	double k=500/fsmpl;
	//Filter < 100Hz
	double FC = 0.2*k;
	for(int i=0;i<M;i++)
	{
		if((i-M/2)==0)
			L[i] = 2*PI*FC;
		if(0>(i-M/2)|0<(i-M/2))
			L[i] = sin(2*PI*FC*(i-M/2))/(i-M/2);
		L[i] = L[i]*(0.54-0.46*cos(2*PI*i/M));
	}
	double sum=0;
	for(i=0;i<M;i++)
	{
		sum = sum+L[i];
	}
	for(i=0;i<M;i++)
	{
		L[i] = L[i]/sum;
	}
	return M;
}

void DSP_Filter::LoadFilterKernels(int fsmpl)
{
	if(fsmpl>500)
		MessageBox(NULL,"Sampling rate<500Hz","Stop",MB_OK);
	CreateWFLow(fsmpl,100);
	CreateWFHigh(fsmpl,200);
	CreateWF50Noht(fsmpl,100);
	CreateWF40Low(fsmpl,100);
}

void DSP_Filter::WFilter_Low100Hz(double *dest, double *source, int lenght, int fsmpl, double cf)
{
	//Win=-sync filter hight < 100hz
	int lenght_1 = lenght+200;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,dest,lenght);
	int M = 100;

	SlideArrRight(Y,X,lenght_1,M/2);
	MoveArr(X,Y,lenght_1);
	for(int i=0;i<lenght_1;i++)
	{
		*(Y+i)=0;
	}
	for(int j=M;j<lenght_1-M;j++)
	{
		for(i=0;i<M;i++)
		{
			*(Y+j)=*(Y+j)+*(X+j-i)*L[i];
		//	*(Y+j) = *(Y+j)+*(X+j+i)*H[i];
		}
	}
	SlideArrLeft(Y,Y,lenght_1,M);//
//	FilterLowFromTo(Y,Y,0,M,lenght);
	MoveArr(dest,Y,lenght);
	//Free memory
	Delete_Arrs();

}

void DSP_Filter::WFilter_Noht(double *dest, double *source, int lenght, int fsmpl, double cf)
{
	//Win=-sync filter hight < 47hz > 53hz
	int lenght_1 = lenght+200;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,dest,lenght);
	int M = 100;

	SlideArrRight(Y,X,lenght_1,M/2);
	MoveArr(X,Y,lenght_1);
	for(int i=0;i<lenght_1;i++)
	{
		*(Y+i)=0;
	}
	for(int j=M;j<lenght_1-M;j++)
	{
		for(i=0;i<M;i++)
		{
			*(Y+j)=*(Y+j)+*(X+j-i)*B1[i];
		//	*(Y+j) = *(Y+j)+*(X+j+i)*H[i];
		}
	}
	SlideArrLeft(Y,Y,lenght_1,M);//
//	FilterLowFromTo(Y,Y,0,M,lenght);
	MoveArr(dest,Y,lenght);
	//Free memory
	Delete_Arrs();
}

int DSP_Filter::CreateWF40Low(int fsmpl, int f_kernel)
{
	int M = f_kernel;
	double k=500/fsmpl;
	double sum = 0;

	//calculate band <41Hz
	double FC = 0.082*k;
	for(int i=0;i<M;i++)
	{
		if((i-M/2)==0)
			L40HZ[i] = 2*PI*FC;
		if(0>(i-M/2)|0<(i-M/2))
			L40HZ[i] = sin(2*PI*FC*(i-M/2))/(i-M/2);
		L40HZ[i] = L40HZ[i]*(0.54-0.46*cos(2*PI*i/M));
	}
	//Calculate unity gain at dc
	sum = 0;
	for(i=0;i<M;i++)
	{
		sum = sum+L40HZ[i];
	}
	//Normalize for unity gain
	for(i=0;i<M;i++)
	{
		L40HZ[i] = L40HZ[i]/sum;
	}

	return M;
}

void DSP_Filter::WFilter_Low40Hz(double *dest, double *source, int lenght, int fsmpl, double cf)
{
	int lenght_1 = lenght+200;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,dest,lenght);
	int M = 100;

	SlideArrRight(Y,X,lenght_1,M/2);
	MoveArr(X,Y,lenght_1);
	for(int i=0;i<lenght_1;i++)
	{
		*(Y+i)=0;
	}
	for(int j=M;j<lenght_1-M;j++)
	{
		for(i=0;i<M;i++)
		{
			*(Y+j)=*(Y+j)+*(X+j-i)*L40HZ[i];
		}
	}
	SlideArrLeft(Y,Y,lenght_1,M);//
	MoveArr(dest,Y,lenght);
	//Free memory
	Delete_Arrs();
}

double DSP_Filter::FindMax(double *source, int lenght)
{
	//Finding the max value from source array
	double max=0;
	max = *(source);
	for(int i=0;i<lenght;i++)
	{
		if(max<*(source+i))
			max = *(source+i);
	}
	return max;
}

double DSP_Filter::FindMin(double *source, int lenght)
{
	//Finding min value from source array
	double min=0;
	min = *(source);
	for(int i=0;i<lenght;i++)
	{
		if(min>*(source+i))
			min = *(source+i);
	}
	return min;

}

void DSP_Filter::Amplifier(double *dest, double *source, int lenght,double ampl)
{
	int lenght_1 = lenght+200;
	Greate_Arrs(lenght_1);
	ZeroArrs(lenght_1);
	MoveArr(X,source,lenght);
	for(int i=0;i<lenght;i++)
	{
		*(Y+i) = *(X+i)*ampl;
	}
	
	MoveArr(dest,Y,lenght);
	//Free memory
	Delete_Arrs();
}

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

License

This article has no explicit license attached to it but may contain usage terms in the article text or the download files themselves. If in doubt please contact the author via the discussion board below.

A list of licenses authors might use can be found here

Share

About the Author

Georgi Petrov
Instructor / Trainer
Bulgaria Bulgaria
PhD, Cum Laude in digital automation systems
M.S. in Telemommunication management
B.S. in Telecommunication systems engineering
Programming: CUDA, C/C++, VHDL
Software and Hardware development and consulting:
data acquisition, image processing, medical instrumentation

| Advertise | Privacy | Mobile
Web03 | 2.8.141015.1 | Last Updated 15 Apr 2004
Article Copyright 2003 by Georgi Petrov
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid