SEIS742-FIR.c

/*	fir.c
*   SEIS742 - May2010
*/
#include <math.h>
//#include <filter.h>
#include "normal.h"
#include "flat.h"
#include "tachy.h"
#include "defib.h"
#include <stdio.h>
#include <string.h>

#define ARRAY_SIZE 1024	// length of the input vector
#define NUM_TAPS 15	// number of filter coefficients
//#define NORMAL_IN_SIZE 1024
typedef enum { false, true } bool;
const long double coefs_doub[NUM_TAPS] = {
	0.000247809018765355,
	0.002273715947318348,
	0.010468194559522206,
	0.031816282592572021,
	0.070964079390059054,
	0.122467424372272940,
	0.168279086570308960,
	0.186803245199647250,
	0.168279086570308960,
	0.122467424372272940,
	0.070964079390059054,
	0.031816282592572021,
	0.010468194559522206,
	0.002273715947318348,
	0.000247809018765355
};

void runFIR(char inputType[50], long double dataArray[], int dataArraySize, 
			int samplesPerSecond);
long double filteredDataMean(long double dataArray[], int dataArraySize);
long double filteredDataStdDeviation(long double dataArray[], 
									 int dataArraySize, long double dataMean);
long double fir_split(long double input, int ntaps, 
					  const long double FIRCoeffs[], 
					  long double delayBuffer[], 
					  int *p_state);

int numTaps = NUM_TAPS;
int main()
{
	char normalLabel[100] = "NORMAL_";
	char defibLabel[100] = "DEFIB_";
	char flatLabel[100] = "FLAT_";
	char tachyLabel[100] = "TACY_";
	printf("************************************************************\r\n");
	printf("***                 NORMAL                               ***\r\n"); 
	runFIR(normalLabel, IN_NORMAL, NORMAL_IN_SIZE, NORMAL_SAMPLES_PER_SECOND);
	printf("************************************************************\r\n");
	printf("\r\n\r\n");
	printf("************************************************************\r\n");
	printf("***                 DEFIB                                ***\r\n");
	runFIR(defibLabel, IN_DEFIB, DEFIB_IN_SIZE, DEFIB_SAMPLES_PER_SECOND);
	printf("************************************************************\r\n");
	printf("\r\n\r\n");
	printf("************************************************************\r\n");
	printf("***                 FLAT                                 ***\r\n");
	runFIR(flatLabel, IN_FLAT, FLAT_IN_SIZE, FLAT_SAMPLES_PER_SECOND);
	printf("************************************************************\r\n");
	printf("\r\n\r\n");
	printf("************************************************************\r\n");
	printf("***                 TACHY                                ***\r\n");
	runFIR(tachyLabel, IN_TACHY, TACHY_IN_SIZE, TACHY_SAMPLES_PER_SECOND);
	printf("************************************************************\r\n");
	printf("%s", defibLabel);

	scanf(" ");
}

/***********************************************************
/* Data analysis
/**********************************************************/
void runFIR(char inputType[50], long double dataArray[], int dataArraySize, 
			int samplesPerSecond){

	int i;
	int filterState=0;	//tracks filter coeficient location
	long double filteredDataArray[ARRAY_SIZE]={0.0};
	long double delayLineArray[NUM_TAPS] = {0};

	int rWaveCount=0; //count of peaks detected in array
	long double rWavePeakLocation[100]={0.0};
	long double deltaT[25]={0.0};
	bool atPeak=false;
	bool slopePositive=false;
	long double slopeValue=0;
	long double peakLevel;	//peak of R wave (mu+sigma)
	long double muValue=0;
	long double sigmaValue=0;

	long double muDeltaT;
	long double sigmaDeltaT;

	/*DEBUG FILE*/
	FILE *ifp, *ofp;
	char filteredOutput[] = "filteredOutput.txt";
	char inputSamples[] = "inputSamples.txt";
	char temp[100] = "";
	char outFileName[100] = "";
	char inFileName[100] = "";
	strcat(temp, inputType);
	strcat(inputType, filteredOutput);
	strcat(outFileName, inputType);
	strcat(temp, inputSamples);
	strcat(inFileName, temp);
	printf("%s", inputSamples);
	ofp = fopen(outFileName, "w");
	ifp = fopen(inFileName, "w");
	if (ifp == NULL) printf("Can't open input file\r\n");
	if (ofp == NULL) printf("Can't open output file\r\n");
	/*END DEBUG FILE*/

	for(i=0;i<dataArraySize;i++)
	{
		//RUN FILTER
		filteredDataArray[i] = fir_split(dataArray[i], numTaps, coefs_doub, 
			delayLineArray, &filterState);

		/* DEBUG - Output to file - Input Data and Filtered Data */ 
		//fprintf(ofp, "%g, %g\n",IN_NORMAL[i], myOutput[i]);
		fprintf(ofp, "%g\n",filteredDataArray[i]);
		fprintf(ifp, "%g\n", dataArray[i]);
		//printf("NormalIn: %f  - Out: %.17g\r\n", dataArray[i], 
		//filteredDataArray[i]);
	}
	fclose(ofp);
	fclose(ifp);

	//find mean - mu
	muValue = filteredDataMean(filteredDataArray, dataArraySize);
	printf("muValue: %.18g\r\n", muValue);
	//find std deviation - sigma
	sigmaValue = filteredDataStdDeviation(filteredDataArray, dataArraySize, 
		muValue);
	printf("sigmaValue: %.18g\r\n", sigmaValue);
	//peak level becomes (mu+sigma)
	peakLevel = muValue+sigmaValue;
	printf("peakLevel: %.18g\r\n", peakLevel);
	/* for loop */
	for(i=NUM_TAPS;i<dataArraySize;i++){
		if((filteredDataArray[i] > peakLevel) && atPeak==0){	
			//found R wave, but not at peak
			//printf("FOUND R WAVE at index: %i -", i);
			//calculate slopes until negative slope is found; 
			//tracing up the r wave
			if(i<dataArraySize) 
				slopeValue = (filteredDataArray[i+1]-filteredDataArray[i]);
			//printf("slopeValue %g\r\n", slopeValue);
			if(slopeValue>0) slopePositive=true;
			else slopePositive = false;

			if(slopePositive && atPeak){
				printf("Climbing R Wave\r\n");
			}
			else { //at peak

				if (!slopePositive && !atPeak){
					//printf("At Peak\r\n");
					if((filteredDataArray[i]-peakLevel)>muValue){
						rWavePeakLocation[rWaveCount] = 
							//holds time of rwave peak
							((long double)i/samplesPerSecond);	
						/*printf("Peak found at %i, time %f, amplitude %g\r\n", 
							i, 
							rWavePeakLocation[rWaveCount], 
							filteredDataArray[i]);*/
						atPeak=true;
						rWaveCount++;
					}
				}
			}
		}
		if(atPeak && (filteredDataArray[i] < peakLevel)){
			atPeak=false;
			slopeValue=0.0;
		}
	}
	/*for loop*/
	//
	printf("R Wave Count: %i\r\n", rWaveCount);

	for(i=0;i<(rWaveCount-1);i++){
		deltaT[i] = rWavePeakLocation[i+1]-rWavePeakLocation[i];
	}

	muDeltaT = filteredDataMean(deltaT, rWaveCount-1);
	sigmaDeltaT = filteredDataStdDeviation(deltaT, rWaveCount-1, muDeltaT);

	printf("Mean R wave peaks: %g\r\n", muDeltaT);
	printf("Std Deviation of dT R wave peaks: %g\r\n", sigmaDeltaT);


	if (sigmaDeltaT > 10) {
		printf("SHOCK!\r\n");
	}
	if (muDeltaT < 0.5 || muDeltaT > 6) {
		printf("SHOCK!\r\n");
	}else
	{
		printf("no shock\r\n");
	}
	return;
}

/****************************************************************************
* 
* Copyright 2000 by Grant R. Griffin
*
* Thanks go to contributors of comp.dsp for teaching me some of these
* techniques, and to Jim Thomas for his review and great suggestions.
*
*                          The Wide Open License (WOL)
*
* Permission to use, copy, modify, distribute and sell this software and its
* documentation for any purpose is hereby granted without fee, provided that
* the above copyright notice and this license appear in all source copies. 
* THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
* ANY KIND. See http://www.dspguru.com/wol.htm for more information.
*
* fir_split: This splits the calculation into two parts so the circular    
* buffer logic doesn't have to be done inside the calculation loop.
*****************************************************************************/
long double fir_split(long double input, int ntaps, 
					  const long double firCoeffs[], long double delayBuffer[],
					  int *p_state)
{
	int ii, end_ntaps, state = *p_state;
	long double accum;
	long double const *p_h;
	long double *p_z;

	/* setup the filter */
	accum = 0;
	p_h = firCoeffs;  //h = filterCoeffs

	/* calculate the end part */
	p_z = delayBuffer + state;
	*p_z = input;
	end_ntaps = ntaps - state;
	for (ii = 0; ii < end_ntaps; ii++) {
		accum += *p_h++ * *p_z++;
	}

	/* calculate the beginning part */
	p_z = delayBuffer;
	for (ii = 0; ii < state; ii++) {
		accum += *p_h++ * *p_z++;
	}

	/* decrement the state, wrapping if below zero */
	if (--state < 0) {
		state += ntaps;
	}
	*p_state = state;       /* pass new state back to caller */

	return accum;
}


//calculate mean of buffer data - mu
long double filteredDataMean(long double dataArray[], int dataArraySize){
	long double accum=0.0;
	int i=0;

	for(i=0;i<dataArraySize;i++){
		accum += dataArray[i];
	}
	return (accum/(long double)dataArraySize);
}

//calculate std deviation of buffer data - sigma
long double filteredDataStdDeviation(long double dataArray[], int dataArraySize,
									 long double dataMean){
	int i=0;
	long double accum=0.0, ii=0.0;

	for(i=0;i<dataArraySize;i++){
		ii = ((dataArray[i]-dataMean));
		accum += (ii*ii);
	}
	return (sqrt((accum/(long double)dataArraySize)));
}