mrkd
1/6/2011 - 10:47 PM
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)));
}