/* file: rrgen.c Entry number : 184 Author(s) : Phil Langley, Emma Bowers, Michael Drinnan, John Allan, Alan Murray Organization : Medical Physics, Freeman Hospital, Newcastle upon Tyne, UK email address : philip.langley at ncl.ac.uk, emma.bowers at nuth.northy.nhs.uk This program should compile without errors or warnings using: gcc -Wall rrgen.c -lm See http://www.physionet.org/challenge/2002/ for further information on the CinC Challenge 2002. This program was used to generate series rr17 and rr25 of the challenge dataset. */ #include #include #include #define SPS 128 /* sampling frequency (determines quantization of RR intervals; see main()) */ #define cSampleRate 5 #define cDT (1.0 / (double) cSampleRate) #define cSamplesPerDay (24 * 60 * 60 * cSampleRate) #define cMaxDays 2 #define cSamplesTotal (cSamplesPerDay * cMaxDays) #define cMaxAlloc (cSamplesTotal + 1000) #define cMaxCycles 1000 #define cPI (3.141592654) #define cDayActiveAmpMaxPos 0.15 #define cDayActiveAmpMinPos 0.05 #define cDayActiveAmpVarPos 0.25 #define cDayActiveFrqMaxPos 0.3 #define cDayActiveFrqMinPos 0.1666 #define cDayActiveFrqVarPos 0.15 #define cDayActiveAmpMaxNeg 0.15 #define cDayActiveAmpMinNeg 0.05 #define cDayActiveAmpVarNeg 0.25 #define cDayActiveFrqMaxNeg 0.3 #define cDayActiveFrqMinNeg 0.1666 #define cDayActiveFrqVarNeg 0.15 #define cDayPassiveAmpMaxPos 0.2 #define cDayPassiveAmpMinPos 0.05 #define cDayPassiveAmpVarPos 0.05 #define cDayPassiveFrqMaxPos 0.2666 #define cDayPassiveFrqMinPos 0.0666 #define cDayPassiveFrqVarPos 1.0 #define cDayPassiveAmpMaxNeg 0.3 #define cDayPassiveAmpMinNeg 0.05 #define cDayPassiveAmpVarNeg 0.15 #define cDayPassiveFrqMaxNeg 0.2666 #define cDayPassiveFrqMinNeg 0.0666 #define cDayPassiveFrqVarNeg 1.0 #define cNightActiveAmpMaxPos 0.3 #define cNightActiveAmpMinPos 0.05 #define cNightActiveAmpVarPos 0.05 #define cNightActiveFrqMaxPos 0.2666 #define cNightActiveFrqMinPos 0.0666 #define cNightActiveFrqVarPos 1.0 #define cNightActiveAmpMaxNeg 0.3 #define cNightActiveAmpMinNeg 0.05 #define cNightActiveAmpVarNeg 0.05 #define cNightActiveFrqMaxNeg 0.2666 #define cNightActiveFrqMinNeg 0.0666 #define cNightActiveFrqVarNeg 1.0 #define cNightPassiveAmpMaxPos 0.2 #define cNightPassiveAmpMinPos 0.05 #define cNightPassiveAmpVarPos 0.1 #define cNightPassiveFrqMaxPos 0.2 #define cNightPassiveFrqMinPos 0.0666 #define cNightPassiveFrqVarPos 0.15 #define cNightPassiveAmpMaxNeg 0.2 #define cNightPassiveAmpMinNeg 0.05 #define cNightPassiveAmpVarNeg 0.1 #define cNightPassiveFrqMaxNeg 0.2 #define cNightPassiveFrqMinNeg 0.0666 #define cNightPassiveFrqVarNeg 0.15 double *nT = NULL ; double *nRR = NULL ; double *nLFeffect = NULL ; double *nWhitenoise = NULL ; double *nPinknoise = NULL ; double *nResp = NULL ; double Min( double A, double B ) { return A < B ? A : B ; } double Max( double A, double B ) { return A > B ? A : B ; } /* Set part of the array (IndexFrom to IndexTo) to equal Value */ void SetArray( double *Array, int IndexFrom, int IndexTo, double Value ) { int Index ; if( ( IndexFrom < 0 ) || ( IndexTo >= cSamplesPerDay ) || ( IndexFrom > IndexTo ) ) fprintf( stderr, "Can't write to array elements %d to %d\n", IndexFrom, IndexTo ) ; else for( Index = IndexFrom ; Index <= IndexTo ; Index ++ ) Array[ Index ] = Value ; } /* Pick a random number between 0 and 1 */ /* rand( ) picks a number between 1 and RAND_MAX */ /* This is scaled down appropriately */ /* Numbers are from a flat distribution */ double Rand0To1( void ) { double Val ; Val = rand( ) ; return Val / (double) RAND_MAX ; } /* Pick a random number from a normal distribution */ /* Approximation to normal distribution works by adding 100 random values from a flat distribution */ /* The value is then scaled for mean of 0 and SD of 1 */ /* The factor 3.464 is square root of 12 */ double RandGauss( void ) { int Count ; double Total ; Total = 0.0 ; for( Count = 0 ; Count < 100 ; Count ++ ) Total += (double) rand( ) ; Total -= 50.0 * (double) RAND_MAX ; Total /= 10.0 * (double) RAND_MAX ; Total *= 3.464 ; return Total ; } /* Convert from Time to the corresponding Index into the array */ int TimeToIndex( double Time ) { double Index ; Index = (double ) Time * (double) cSampleRate ; return (int) Index ; } /* Convert from Index to the corresponding Time */ double IndexToTime( int Index ) { double Time ; Time = (double ) Index / (double) cSampleRate ; return Time ; } /* % a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) % - a(2)*y(n-1) - ... - a(na+1)*y(n-na) */ void Filter( double *X, int N, double *a, double *b, int Order, double *Y ) { int SampleNo ; int SourceSampleNo ; int TapNo ; double Output ; for( SampleNo = 0 ; SampleNo < N ; SampleNo ++ ) { Output = b[ 0 ] * X[ SampleNo ] ; for( TapNo = 1 ; TapNo <= Order ; TapNo ++ ) { SourceSampleNo = SampleNo - TapNo ; if( SourceSampleNo < 0 ) break ; Output += b[ TapNo ] * X[ SourceSampleNo ] - a[ TapNo ] * Y[ SourceSampleNo ] ; } Y[ SampleNo ] = Output ; } } /* Implement a first-order low-pass filter for 1/f rolloff */ /* Idea is first to calculate the corresponding time constant */ /* This leads to multiplication factor for exponential smoothing */ /* This is equivalent to a 1st order RC filter */ /* General idea is: NewOutput = ( Z * OldOutput ) + ( ( 1 - Z ) * CurrentInput ) */ /* If Z is nearly 1, the time constant is long */ void FirstOrderLowPassFilter( double *X, int N, double CutoffFreq, double *Y ) { double TimeConstant ; double DecayPerSec ; double DecayPerSample ; double OneMinusDecayPerSample ; double Output ; int SampleNo ; TimeConstant = 1.0 / ( 2.0 * cPI * CutoffFreq ) ; DecayPerSec = 1.0 / TimeConstant ; DecayPerSample = DecayPerSec / (double) cSampleRate ; OneMinusDecayPerSample = 1.0 - DecayPerSample ; Output = 0.0 ; for( SampleNo = 1 ; SampleNo <= N ; SampleNo ++ ) { Output = ( OneMinusDecayPerSample * Output ) + ( DecayPerSample * X[ SampleNo ] ) ; Y[ SampleNo ] = Output ; } } /* Sort the values in the Raw array */ /* Number of values is given by N */ /* Results are put in Cooked array */ /* Uses simple insertion sort algorithm */ /* Go thru Raw array N times */ /* Each time thru, find the smallest element */ /* Add this element to the end of the Cooked array */ /* Then set the element to a very big number so it doesn't get found again */ void Sort( double *Raw, int N, double *Cooked ) { int Iteration ; int Smallest ; int Test ; for( Iteration = 0 ; Iteration <= N - 1 ; Iteration ++ ) { Smallest = 1 ; for( Test = 1 ; Test <= N ; Test ++ ) if( Raw[ Test ] < Raw[ Smallest ] ) Smallest = Test ; Cooked[ Iteration ] = Raw[ Smallest ] ; Raw[ Smallest ] = 1e250 ; } } /* This section adapted from MATLAB code by Phil & Emma */ /* %development code for Computers in Cardiology Challenge 2002 %24 hr RR simulator for normal subjects %CVPER group %16th April 2002 */ /*Calculating the frequency and amplitude of the respiratory component*/ void choosefreqamp( double maxamp, double minamp, double percentagechangeamp, double previousamp, double maxfreq, double minfreq, double variance, double tlength, double i, double *numberoftimeelement, double *freq, double *amp ) { double meanfreq ; /*choosing the amplitude randomly within a percentage range of the previous value*/ *amp = ( Rand0To1( ) * 2.0 * percentagechangeamp * previousamp ) + ( ( 1.0 - percentagechangeamp ) * previousamp ) ; /*if the above value does not fall in the given range guess again */ *amp = Min( *amp, maxamp + Rand0To1( ) * 0.02 - 0.01 ) ; *amp = Max( *amp, minamp + Rand0To1( ) * 0.02 - 0.01 ) ; /*choosing the frequency randomly*/ meanfreq = ( maxfreq + minfreq ) / 2.0 ; *freq = RandGauss( ) * variance + meanfreq ; /*if the above value does not fall in the given range continue to guess until it does*/ while( ( *freq < minfreq ) || ( *freq > maxfreq ) ) *freq = RandGauss( ) * variance + meanfreq ; /*working out how many whole array elements are in a half period*/ *numberoftimeelement = floor( 1.0 / ( 2.0 * *freq * cDT ) ) ; if( i + *numberoftimeelement > tlength ) *numberoftimeelement = tlength - i ; } /* Note this only generates 1 day's worth of RR data */ void RRSim1Day( double RRday, double RRnight, double *T, double *RR, double *Resp ) { int Index ; int SleepIndex ; int WakeIndex ; /* int NumRemCycles ;*/ int NumFalls1 ; int NumFalls2 ; int NumArouse ; int NumArtifacts ; /* int Cycle ;*/ int Fall ; int Arouse ; int Artifact ; /* int RemIndex ;*/ int FallIndex1 ; int FallIndex2 ; int ArouseIndex ; int i ; int IsSleep ; int IsFall1 ; int IsFall2 ; int IsArousal ; /* int IsREM ;*/ double SleepTime ; double WakeTime ; /* double RemTime ;*/ double ArouseTime ; double RampDuration ; double Duration ; double FallTime ; double FallRR ; /* double RemRR ;*/ double ArouseRR ; double FreqWander ; double RandInd ; double amp ; double numberoftimeelement ; /* double CoeffA[ ] = { 1.0000, -4.0469, 6.3705, -4.8224, 1.7212, -0.2223 } ; double CoeffB[ ] = { 1.0000, -3.4673, 4.4317, -2.4428, 0.4608, 0.0177 } ;*/ /* double RemDuration[ cMaxCycles ] ; double StableDuration[ cMaxCycles ] ; double RemFrom[ cMaxCycles ] ; double RemTo[ cMaxCycles ] ;*/ double FallFrom1[ cMaxCycles ] ; double FallTo1[ cMaxCycles ] ; double FallFrom2[ cMaxCycles ] ; double FallTo2[ cMaxCycles ] ; double ArouseFrom[ cMaxCycles ] ; double ArouseTo[ cMaxCycles ] ; double OrderedFallFrom1[ cMaxCycles ] ; double OrderedFallTo1[ cMaxCycles ] ; double OrderedArouseFrom[ cMaxCycles ] ; double OrderedArouseTo[ cMaxCycles ] ; double OrderedFallFrom2[ cMaxCycles ] ; double OrderedFallTo2[ cMaxCycles ] ; double freqpos ; double freqneg ; double amppos; double ampneg ; static double previousamppos = 0.1 ; static double previousampneg = 0.1 ; /* %define time axis t = [0:dt:24*60*60]; */ for( Index = 0 ; Index < cSamplesPerDay ; Index ++ ) T[ Index ] = (double) Index * cDT ; /* %t = 0 is start of recording %Define sleep and awake times %sleep 12 to 14 hrs after start time sleepTime = (rand*2+12)*60*60; */ SleepTime = ( Rand0To1( ) * 2.0 + 12.0 ) * 60.0 * 60.0 ; /* %wake up between 2 and 3 hrs before end of recording next day wakeTime = (24-(rand*1+2))*60*60; */ WakeTime = ( 24.0 - ( Rand0To1( ) * 1.0 + 2.0 ) ) * 60.0 * 60.0 ; /* RR = ones(size(t))*RRday; sleepIndex = find(t > (sleepTime-dt/2) & t < (sleepTime+dt/2)); wakeIndex = find(t > (wakeTime-dt/2) & t < (wakeTime+dt/2)); RR(sleepIndex(1):wakeIndex(1)) = RRnight; */ SleepIndex = TimeToIndex( SleepTime ) ; WakeIndex = TimeToIndex( WakeTime ) ; for( Index = 0 ; Index < cSamplesPerDay ; Index ++ ) RR[ Index ] = RRday ; for( Index = SleepIndex ; Index <= WakeIndex ; Index ++ ) RR[ Index ] = RRnight ; /* %ramp up from day to night rr over 30 min period %y = mx+c if(1) rampDuration = 30*60; RR(sleepIndex(1):sleepIndex(1)+rampDuration/dt) = ((RRnight-RRday)/(rampDuration/dt))*[0:rampDuration/dt] + RRday; end */ /* if( 1 ) { */ /* fprintf( stderr, " Day-night ramp...\n" ) ;*/ RampDuration = (Rand0To1( ) * 50.0 + 40.0 ) * 60.0 ; for( Index = 0 ; Index <= RampDuration / cDT ; Index ++ ) RR[ SleepIndex + Index ] = ( ( RRnight - RRday ) / ( RampDuration / cDT ) ) * Index + RRday ; /* } */ RampDuration = (Rand0To1( ) * 15.0 + 15.0 ) * 60.0 ; for( Index = 0 ; Index <= RampDuration / cDT ; Index ++ ) RR[ WakeIndex + Index ] = ( ( RRday - RRnight ) / ( RampDuration / cDT ) ) * Index + RRnight ; /* %cycle of 50 to 70 mins of stable sleep RR (RRnight) followed by 15 to 45 mins unstable sleep RR (REM sleep) %reduces mean RR during REM sleep by between 5 to 20 % if (1) numCycles = fix((wakeTime-sleepTime)/(1.5*60*60)); for (cycle = 1:numCycles) remDuration(cycle) = (rand*30+15)*60; stableDuration(cycle) = (rand*20+50)*60; if (cycle == 1) remTime = sleepTime+stableDuration(cycle); else remTime = remTime+stableDuration(cycle)+remDuration(cycle-1); end remRR = RRnight*(1-(rand*0.15+0.05)); remInd(cycle).vect = find(t > remTime & t < remTime+remDuration(cycle)); RR(remInd(cycle).vect) = ones(size(remInd(cycle).vect))*remRR; end end */ /* if( 1 ) { */ // fprintf( stderr, " REM cycles...\n" ) ; /* NumRemCycles = (int) floor( ( WakeTime - SleepTime ) / ( 1.5 * 60.0 * 60.0 ) ) ; for( Cycle = 1 ; Cycle <= NumRemCycles ; Cycle ++ ) { RemDuration[ Cycle ] = ( Rand0To1( ) * 30.0 + 15.0 ) * 60.0 ; StableDuration[ Cycle ] = ( Rand0To1( ) * 40.0 + 40.0 ) * 60.0 ; if( Cycle == 1 ) RemTime = SleepTime + StableDuration[ Cycle ] ; else RemTime = RemTime + StableDuration[ Cycle ] + RemDuration[ Cycle - 1 ] ; RemRR = RRnight * ( 1 - ( Rand0To1( ) * 0.03 +0.02 ) ) ; RemFrom[ Cycle ] = TimeToIndex( RemTime ) ; RemTo[ Cycle ] = TimeToIndex( RemTime + RemDuration[ Cycle ] ) ; SetArray( RR, RemFrom[ Cycle ], RemTo[ Cycle ], RemRR ) ; }*/ /* } */ /* %2 to 5 sudden falls in mean RR of 10 to 20 % between start of recording and sleep %with durations of 3 to 30 minutes to simulate various physical activities! if (1) numFalls = fix(rand*3)+2; for (fall = 1:numFalls) duration = (rand*27+3)*60; fallTime = rand*sleepTime; fallRR = RRday*(1-(rand*0.1+0.1)); fallInd1(fall).vect = find(t > fallTime & t < fallTime+duration); RR(fallInd1(fall).vect) = ones(size(fallInd1(fall).vect))*fallRR; end end */ /* if( 1 ) { */ // fprintf( stderr, " Falls 1...\n" ) ; NumFalls1 = (int) floor( Rand0To1( ) * 3.0 ) + 2.0 ; for( Fall = 1 ; Fall < NumFalls1 ; Fall ++ ) { Duration = ( Rand0To1( ) * 27.0 + 3.0 ) * 60.0 ; FallTime = Rand0To1( ) * SleepTime ; FallRR = RRday * ( 1.0 - ( Rand0To1( ) * 0.1 + 0.1 ) ) ; FallFrom1[ Fall ] = TimeToIndex( FallTime ) ; FallTo1[ Fall ] = TimeToIndex( FallTime + Duration ) ; SetArray( RR, FallFrom1[ Fall ], FallTo1[ Fall ], FallRR ) ; } /* } */ /* %1 to 2 sudden falls in mean RR of 10 to 20 % between awaking and end of recording %with durations of 3 to 30 minutes to simulate various physical activities! if (1) numFalls = fix(rand)+1; for (fall = 1:numFalls) duration = (rand*27+3)*60; fallTime = rand*(24*60*60-wakeTime)+wakeTime; fallRR = RRday*(1-(rand*0.1+0.1)); fallInd2(fall).vect = find(t > fallTime & t < fallTime+duration); RR(fallInd2(fall).vect) = ones(size(fallInd2(fall).vect))*fallRR; end end */ /* if( 1 ) { */ // fprintf( stderr, " Falls 2...\n" ) ; NumFalls2 = (int) floor( Rand0To1( ) ) + 1.0 ; for( Fall = 1 ; Fall < NumFalls2 ; Fall ++ ) { Duration = ( Rand0To1( ) * 27.0 + 3.0 ) * 60.0 ; FallTime = Rand0To1( ) * ( 24.0 * 60.0 * 60.0 - WakeTime ) + WakeTime ; FallRR = RRday * ( 1.0 - ( Rand0To1( ) * 0.1 + 0.1 ) ) ; FallFrom2[ Fall ] = TimeToIndex( FallTime ) ; FallTo2[ Fall ] = TimeToIndex( FallTime + Duration ) ; SetArray( RR, FallFrom2[ Fall ], FallTo2[ Fall ], FallRR ) ; } /* } */ /* %15 to 25 spontaneous arousels of 5 to 20 s duration during sleep %where mean sleep RR falls to by 10 to 20 % if (1) numArouse = fix(rand*10)+15; for (arouse = 1:numArouse) duration = rand*15+5; arouseTime = rand*(wakeTime-sleepTime)+sleepTime; arouseRR = RRnight*(1-(rand*0.1+0.1)); arouseInd(fall).vect = find(t > arouseTime & t < arouseTime+duration); RR(arouseInd(fall).vect) = ones(size(arouseInd(fall).vect))*arouseRR; end end */ /* if( 1 ) { */ // fprintf( stderr, " Arousals...\n" ) ; NumArouse = floor( Rand0To1( ) * 10.0 ) + 15.0 ; for( Arouse = 1 ; Arouse <= NumArouse ; Arouse ++ ) { Duration = Rand0To1( ) * 15.0 + 5.0 ; ArouseTime = Rand0To1( ) * ( WakeTime - SleepTime ) + SleepTime ; ArouseRR = RRnight * ( 1.0 - ( Rand0To1( ) * 0.1 + 0.1 ) ) ; ArouseFrom[ Arouse ] = TimeToIndex( ArouseTime ) ; ArouseTo[ Arouse ] = TimeToIndex( ArouseTime + Duration ) ; SetArray( RR, ArouseFrom[ Arouse ], ArouseTo[ Arouse ], ArouseRR ) ; } /* } */ /* if (1) %Baseline wander - low frequency 2 to 5 hrs freqWander = 1/((rand*3+2)*60*60); sleepIndex = fix(sleepTime/dt); amp = 0.1; RR(1:sleepIndex) = RR(1:sleepIndex)+amp*sin(2*pi*freqWander*t(1:sleepIndex)); end */ /* if( 1 ) { */ // fprintf( stderr, " Baseline wander...\n" ) ; // FreqWander = 1.0 / ( ( Rand0To1( ) * 3.0 + 2.0 ) * 24 * 60.0 * 60.0 ) ; FreqWander = 1.0/(13.0*60.0*60.0); SleepIndex = floor( SleepTime / cDT ) ; amp = 0.1 ; for( Index = 1 ; Index <= SleepIndex ; Index ++ ) RR[ Index ] += amp * (sin( 2.0 * cPI * FreqWander * IndexToTime( Index )+(cPI/6.5) ) +sin( 2.0 * cPI * (FreqWander-(1.0/(5.0*60.0*60.0))) * IndexToTime( Index ) ) +sin( 2.0 * cPI * (FreqWander-(1.0/(8.0*60.0*60.0))) * IndexToTime( Index ) )); /* } */ /* if (1) %measurement artifact %0 or 5 artifacts during period of activity with high (1.5 to 2) rr for (fall = 1:size(fallInd1)) numArtifacts = fix(rand*5); for (artifact = 1:numArtifacts) randInd = fix(rand*(fallInd1(fall).vect(end)-fallInd1(fall).vect(1))+1); ind = fallInd1(fall).vect(randInd); RR(ind) = rand*0.5+1.5; end end end */ for( Fall = 1 ; Fall <= NumFalls1 ; Fall ++ ) { // fprintf( stderr, " Artifacts...\n" ) ; NumArtifacts = floor( Rand0To1( ) * 5.0 ) ; for( Artifact = 1 ; Artifact <= NumArtifacts ; Artifact ++ ) { RandInd = floor( Rand0To1( ) * ( FallTo1[ Fall ] - FallFrom1[ Fall ] ) + 1.0 ) ; RR[ (int)RandInd ] = Rand0To1( ) * 0.5 + 1.5 ; } } /* This is where Emma's respiration bit starts */ /* if( 1) {*/ // fprintf( stderr, " Emma's Resp bit...\n" ) ; /*initialise variables */ SetArray( Resp, 1, cSamplesPerDay - 1, 0 ) ; i=1 ; FallIndex1 = 1 ; FallIndex2 = 1 ; /* RemIndex = 1 ;*/ ArouseIndex = 1 ; /*ordering the random activities and arousal*/ Sort( FallFrom1, NumFalls1, OrderedFallFrom1 ) ; Sort( FallTo1, NumFalls1, OrderedFallTo1 ) ; Sort( ArouseFrom, NumArouse, OrderedArouseFrom ) ; Sort( ArouseTo, NumArouse, OrderedArouseTo ) ; Sort( FallFrom2, NumFalls2, OrderedFallFrom2 ) ; Sort( FallTo2, NumFalls2, OrderedFallTo2 ) ; /* while (1) */ while( 1 ) { /* calculating what is happening at each time instant*/ IsSleep = ( ( T[ i ] > SleepTime ) && ( T[ i ] < WakeTime ) ) ; IsFall1 = ( i <= OrderedFallTo1[ FallIndex1 ] ) && ( i >= OrderedFallFrom1[ FallIndex1 ] ) ; IsFall2 = ( i <= OrderedFallTo2[ FallIndex2 ] ) && ( i >= OrderedFallFrom2[ FallIndex2 ] ) ; IsArousal = ( i <= OrderedArouseTo[ ArouseIndex ] ) && ( i >= OrderedArouseFrom[ ArouseIndex ] ) ; /* IsREM = ( i <= RemTo[ RemIndex ] ) && ( i >= RemFrom[ RemIndex ] ) ;*/ if( !IsSleep && IsFall1 ) { /*this is day (before sleep) and active %so breathing rate increases so RR changes are modulated by a faster less varying frequency %the positive part of the cycle %calculating the amp, freq and number of points the half cycle passes through*/ // fprintf( stderr, " Emma's Resp bit...\n" ) ; choosefreqamp( cDayActiveAmpMaxPos, cDayActiveAmpMinPos, cDayActiveAmpVarPos, previousamppos, cDayActiveFrqMaxPos, cDayActiveFrqMinPos, cDayActiveFrqVarPos, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; /*appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = amppos * sin( 2.0 * cPI * freqpos * T[ Index + 1 ] ) ; i += numberoftimeelement ; /*fall out of the loop if the array is full */ if( i >= cSamplesPerDay ) break ; /*remembering previous values so changes can be related to previous values */ previousamppos = amppos ; /*negative part of the cycle %calculating the amp, freq and number of points the half cycle passes through */ choosefreqamp( cDayActiveAmpMaxNeg, cDayActiveAmpMinNeg, cDayActiveAmpVarNeg, previousampneg, cDayActiveFrqMaxNeg, cDayActiveFrqMinNeg, cDayActiveFrqVarNeg, cSamplesPerDay, i, &numberoftimeelement, &freqneg, &neg ) ; /*%appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = -ampneg * sin( 2.0 * cPI * freqneg *T[ Index + 1 ] ) ; i += numberoftimeelement ; /*%fall out of the loop if the array is full*/ if( i >= cSamplesPerDay ) break ; /*remembering previous values so changes can be related to previous values */ previousampneg = ampneg ; /*got to the end of the active time look at next one*/ if( ( i >= OrderedFallTo1[ FallIndex1 ] ) && ( FallIndex1 != NumFalls1 ) ) FallIndex1 ++ ; } /* end if */ else if( !IsSleep && IsFall2 ) { /*this is day (after wake up) and active %so breathing rate increases so RR changes are modulated by a faster less varying frequency %the positive part of the cycle %calculating the amp, freq and number of points the half cycle passes through*/ fprintf( stderr, " Emma's Resp bit2...\n" ) ; choosefreqamp( cDayActiveAmpMaxPos, cDayActiveAmpMinPos, cDayActiveAmpVarPos, previousamppos, cDayActiveFrqMaxPos, cDayActiveFrqMinPos, cDayActiveFrqVarPos, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; /*appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = amppos * sin( 2.0 * cPI * freqpos * T[ Index + 1 ] ) ; i += numberoftimeelement ; /*fall out of the loop if the array is full */ if( i >= cSamplesPerDay ) break ; /*remembering previous values so changes can be related to previous values */ previousamppos = amppos ; /*negative part of the cycle %calculating the amp, freq and number of points the half cycle passes through */ choosefreqamp( cDayActiveAmpMaxNeg, cDayActiveAmpMinNeg, cDayActiveAmpVarNeg, previousampneg, cDayActiveFrqMaxNeg, cDayActiveFrqMinNeg, cDayActiveFrqVarNeg, cSamplesPerDay, i, &numberoftimeelement, &freqneg, &neg ) ; /*%appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = -ampneg * sin( 2.0 * cPI * freqneg *T[ Index + 1 ] ) ; i += numberoftimeelement ; /*%fall out of the loop if the array is full*/ if( i >= cSamplesPerDay ) break ; /*remembering previous values so changes can be related to previous values */ previousampneg = ampneg ; /*got to the end of the active time look at next one*/ if( ( i >= OrderedFallTo2[ FallIndex2 ] ) && ( FallIndex2 != NumFalls2 ) ) FallIndex2 ++ ; } /* end if */ else if( !IsSleep ) { /*%this is day and not active %large variation in breathing frequency so RR interval changes are modulated by a large number of frequencies %the positive part of the cycle %calculating the amp, freq and number of points the half cycle passes through*/ choosefreqamp( cDayPassiveAmpMaxPos, cDayPassiveAmpMinPos , cDayPassiveAmpVarPos, previousamppos, cDayPassiveFrqMaxPos, cDayPassiveFrqMinPos, cDayPassiveFrqVarPos, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; /*appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = amppos * sin( 2.0 * cPI * freqpos * T[ Index + 1 ] ) ; i += numberoftimeelement ; /*%fall out of the loop if the array is full*/ if( i >= cSamplesPerDay ) break ; /*%remembering previous values so changes can be related to previous values*/ previousamppos = amppos ; /*negative part of the cycle %calculating the amp, freq and number of points the half cycle passes through*/ choosefreqamp( cDayPassiveAmpMaxNeg, cDayPassiveAmpMinNeg, cDayPassiveAmpVarNeg, previousampneg, cDayPassiveFrqMaxNeg, cDayPassiveFrqMinNeg, cDayPassiveFrqVarNeg, cSamplesPerDay, i, &numberoftimeelement, &freqneg, &neg ) ; /*appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = -ampneg * sin( 2.0 * cPI * freqneg * T[ Index + 1 ] ) ; i += numberoftimeelement ; /*%fall out of the loop if the array is full*/ if( i >= cSamplesPerDay ) break ; /*%remembering previous values so changes can be related to previous values*/ previousampneg = ampneg ; } /* end else if */ else if( IsSleep && IsArousal ) { /*%this is night and aroused %when aroused breathing frequency increases so RR changes are modulated by a faster frequency %the positive part of the cycle %calculating the amp, freq and number of points the half cycle passes through*/ choosefreqamp( cNightActiveAmpMaxPos, cNightActiveAmpMinPos, cNightActiveAmpVarPos, previousamppos, cNightActiveFrqMaxPos, cNightActiveFrqMinPos, cNightActiveFrqVarPos, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; /*%appending samples to the end of the resp signal*/ for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = amppos * sin( 2.0 * cPI * freqpos * T[ Index + 1 ] ) ; i += numberoftimeelement ; /*%fall out of the loop if the array is full*/ if( i >= cSamplesPerDay ) break ; /*%remembering previous values so changes can be related to previous values*/ previousamppos = amppos ; /* %appending samples to the end of the resp signal resp(i:i+numberoftimeelement-1)=-ampneg*sin(2*pi*freqneg*t(1:numberoftimeelement)); %incrementing i i=i+numberoftimeelement; %fall out of the loop if the array is full if (i >= length(t)) break; end %remembering previous values so changes can be related to previous values previousampneg=ampneg; */ /*negative part of the cycle %calculating the amp, freq and number of points the half cycle passes through*/ choosefreqamp( cNightActiveAmpMaxNeg, cNightActiveAmpMinNeg, cNightActiveAmpVarNeg, previousampneg, cNightActiveFrqMaxNeg, cNightActiveFrqMinNeg, cNightActiveFrqVarNeg, cSamplesPerDay, i, &numberoftimeelement, &freqneg, &neg ) ; for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = -ampneg * sin( 2.0 * cPI * freqneg * T[ Index + 1 ] ) ; i += numberoftimeelement ; if( i >= cSamplesPerDay ) break ; previousampneg = ampneg ; /* %got to the end of the aroused time look at next one if (t(i)>=t(orderarouseend(ArouseIndex)) & ArouseIndex~=numArouse) ArouseIndex=ArouseIndex+1; end */ if( ( i >= OrderedArouseTo[ ArouseIndex ] ) && ( ArouseIndex != NumArouse ) ) ArouseIndex ++ ; } /* end else if */ /* elseif ((t(i)>sleepTime) & (t(i)=t(orderremstart(RemIndex))) %this is night and REM %when aroused breathing frequency increases so RR changes are modulated by a faster frequency %the positive part of the cycle %calculating the amp, freq and number of points the half cycle passes through [numberoftimeelement,freqpos,amppos]=choosefreqamp(0.2,0.1,0.15,previousamppos,0.2666,0.0666,1,dt,length(t),i); %appending samples to the end of the resp signal resp(i:i+numberoftimeelement-1)= amppos*sin(2*pi*freqpos*t(1:numberoftimeelement)); %incrementing i i=i+numberoftimeelement; %fall out of the loop if the array is full if (i >= length(t)) break; end %remembering previous values so changes can be related to previous values previousamppos=amppos; */ /* else if( IsSleep && IsREM ) { choosefreqamp( cNightActiveAmpMaxPos, cNightActiveAmpMinPos, cNightActiveAmpVarPos, previousamppos, cNightActiveFrqMaxPos, cNightActiveFrqMinPos, cNightActiveFrqVarPos, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = amppos * sin( 2.0 * cPI * freqpos * T[ Index + 1 ] ) ; i += numberoftimeelement ; if( i >= cSamplesPerDay ) break ; previousamppos = amppos ;*/ /* %negative part of the cycle %calculating the amp, freq and number of points the half cycle passes through [numberoftimeelement,freqneg,ampneg]=choosefreqamp(0.2,0.1,0.15,previousampneg,0.2666,0.0666,1,dt,length(t),i); %appending samples to the end of the resp signal resp(i:i+numberoftimeelement-1)=-ampneg*sin(2*pi*freqneg*t(1:numberoftimeelement)); %incrementing i i=i+numberoftimeelement; %fall out of the loop if the array is full if (i >= length(t)) break; end %remembering previous values so changes can be related to previous values previousampneg=ampneg; */ /* choosefreqamp( cNightActiveAmpMaxNeg, cNightActiveAmpMinNeg, cNightActiveAmpVarNeg, previousamppos, cNightActiveFrqMaxNeg, cNightActiveFrqMinNeg, cNightActiveFrqVarNeg, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = -ampneg * sin( 2.0 * cPI * freqneg * T[ Index + 1 ] ) ; i += numberoftimeelement ; if( i >= cSamplesPerDay ) break ; previousampneg = ampneg ;*/ /* %got to the end of the aroused time look at next one if (t(i)>=t(orderremend(RemIndex)) & RemIndex~=numRem) RemIndex=RemIndex+1; end */ /* if( ( i >= RemTo[ RemIndex ] ) && ( RemIndex != NumRemCycles ) ) RemIndex ++ ; }*/ /* end else if */ /* elseif ((t(i)>sleepTime) & (t(i)= length(t)) break; end %remembering previous values so changes can be related to previous values previousamppos=amppos; */ else if( IsSleep ) { choosefreqamp( cNightPassiveAmpMaxPos, cNightPassiveAmpMinPos, cNightPassiveAmpVarPos, previousamppos, cNightPassiveFrqMaxPos, cNightPassiveFrqMinPos, cNightPassiveFrqVarPos, cSamplesPerDay, i, &numberoftimeelement, &freqpos, &pos ) ; for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = amppos * sin( 2.0 * cPI * freqpos * T[ Index + 1 ] ) ; i += numberoftimeelement ; if( i >= cSamplesPerDay ) break ; previousamppos = amppos ; /* %negative part of the cycle %calculating the amp, freq and number of points the half cycle passes through [numberoftimeelement,freqneg,ampneg]=choosefreqamp(0.3,0.05,0.15,previousampneg,0.2,0.0666,0.15,dt,length(t),i); %appending samples to the end of the resp signal resp(i:i+numberoftimeelement-1)=-ampneg*sin(2*pi*freqneg*t(1:numberoftimeelement)); %incrementing i i=i+numberoftimeelement; %fall out of the loop if the array is full if (i >= length(t)) break; end %remembering previous values so changes can be related to previous values previousampneg=ampneg; */ choosefreqamp( cNightPassiveAmpMaxNeg, cNightPassiveAmpMinNeg, cNightPassiveAmpVarNeg, previousampneg, cNightPassiveFrqMaxNeg, cNightPassiveFrqMinNeg, cNightPassiveFrqVarNeg, cSamplesPerDay, i, &numberoftimeelement, &freqneg, &neg ) ; for( Index = 0 ; Index < numberoftimeelement ; Index ++ ) Resp[ Index + i ] = -ampneg * sin( 2.0 * cPI * freqneg * T[ Index + 1 ] ) ; i += numberoftimeelement ; if( i >= cSamplesPerDay ) break ; previousampneg = ampneg ; } /* end else if */ } /* end while( 1 ) */ for( Index = 0 ; Index < cSamplesPerDay ; Index ++ ) RR[ Index ] += Resp[ Index ] ; /* } end if( 1 ) */ } /* function rrsim(dummy) */ /* Calls RRSim1Day as many times as necessary */ void RRSim( void ) { int Day ; int StartIndex ; int Index ; double RRday ; double RRnight ; double amp ; double phase ; double deltaphase ; double freq1 ; /* double ThisFreq ;*/ /* %Define mean night and day RR %day RR between 0.6 & 0.9 s %night RR between 0.8 & 1.2 s %difference between night & day RR greater than 0.2 s while (1) RRday = rand*0.3+0.6; RRnight = rand*0.4+0.8; if ((RRnight - RRday) > 0.2) break; end end */ // fprintf( stderr, "\nSimulating background...\n" ) ; while( 1 ) { RRday = Rand0To1( ) * 0.3 + 0.8 ; RRnight = Rand0To1( ) * 0.4 + 0.8 ; if ( ( RRnight - RRday ) > 0.2 ) break ; } /* if (1) %0.1 Hz effect %amplitude 5 to 10 % of mean daytime RR %frequency 0.05 to 0.2 Hz - multiple components to smooth peak out freq1 = (rand*0.15)+0.05; amp = RRday*(rand*0.008+0.008); freq=[freq1-0.005:0.001:freq1+0.005]; lfeffect = zeros(size(t)); for (freqNo = 1:length(freq)) lfeffect = lfeffect+amp*sin(2*pi*freq(freqNo)*t); end RR = RR+lfeffect; subplot(4,1,2),plot([0:dt/(60*60):24],RR); set(gca,'xtick',[0:6:24]); xlim([0 24]); end */ /* if( 0 ) { fprintf( stderr, " Phil's 0.1 Hz effect...\n" ) ; freq1 = ( Rand0To1( ) * 0.15 ) + 0.05 ; for( Index = 0 ; Index < cSamplesTotal ; Index ++ ) nLFeffect[ Index ] = 0.0 ; for( ThisFreq = freq1 - 0.005 ; ThisFreq <= freq1 + 0.005 ; ThisFreq +=0.001 ) { amp = RRday * ( Rand0To1( ) * 0.008 + 0.008 ) ; phase = Rand0To1( ) * 2.0 * cPI ; for( Index = 0 ; Index < cSamplesPerDay ; Index ++ ) nLFeffect[ Index ] = nLFeffect[ Index ] + amp * sin( 2.0 * cPI * ThisFreq * IndexToTime( Index ) + phase ); } for( Index = 0 ; Index < cSamplesTotal ; Index ++ ) nRR[ Index ] += nLFeffect[ Index ] ; } */ /* My effort at a 0.1 Hz effect */ /* if( 1 ) { */ // fprintf( stderr, " MD's 0.1 Hz effect...\n" ) ; amp = RRday * ( Rand0To1( ) * 0.008 + 0.008 ) ; freq1 = ( Rand0To1( ) * 0.15 ) + 0.05 ; for( Index = 0 ; Index < cSamplesTotal ; Index ++ ) nLFeffect[ Index ] = 0.0 ; phase = 0.0 ; for( Index = 0 ; Index < cSamplesTotal ; Index ++ ) { if( Rand0To1( ) < 0.01 ) freq1 += ( Rand0To1( ) - 0.5 ) * 0.01 ; if( freq1 == 0.2 ) freq1 = 0.17 ; if( freq1 == 0.05 ) freq1 = 0.08 ; deltaphase = 2.0 * cPI * freq1 * cDT ; amp += Rand0To1( ) * 0.002 - 0.001 ; amp = Min( amp, 0.016 ) ; amp = Max( amp, 0.008 ) ; phase += deltaphase * ( 0.95 + Rand0To1( ) * 0.1 ) ; nLFeffect[ Index ] = amp * sin( phase ); } for( Index = 0 ; Index < cSamplesTotal ; Index ++ ) nRR[ Index ] += nLFeffect[ Index ] ; /* } */ /* if (1) %1/f effect (pink noise) %amplitude 1 to 2 % of mean daytime RR whitenoise=randn(size(t)); %poles and zero from Mr. Tom Bruhns of HP %for white to pink noise filter for audio signals %poles=[0.9986823 0.9914651 0.9580812 0.8090598 0.2896591]'; %zeros=[0.9963594 0.9808756 0.9097290 0.6128445 -0.0324723]'; %get the numerator and denominator of the filter %[b,a]=zp2tf(zeros,poles,1); a = [ 1.0000 -4.0469 6.3705 -4.8224 1.7212 -0.2223]; b = [ 1.0000 -3.4673 4.4317 -2.4428 0.4608 0.0177]; %filter to get pink noise pinknoise=filter(b,a,whitenoise); amp = RRday*(rand*0.025+0.025); RR = RR+amp*pinknoise; end */ /* if( 1 ) { */ // fprintf( stderr, " White noise...\n" ) ; for( Index = 0 ; Index < cSamplesTotal ; Index ++ ) nWhitenoise[ Index ] = RandGauss( ) ; // Filter( whitenoise, cSamplesPerDay, CoeffA, CoeffB, 5, Pinknoise ) ; // fprintf( stderr, " Pink noise...\n" ) ; FirstOrderLowPassFilter( nWhitenoise, cSamplesTotal, 0.1, nPinknoise ) ; amp = RRday * ( Rand0To1( ) * 0.025 + 0.025 ) ; for( Index = 0 ; Index < cSamplesPerDay ; Index ++ ) nRR[ Index ] += amp * nPinknoise[ Index ] ; /* } */ for( Day = 0 ; Day < cMaxDays ; Day ++ ) { // fprintf( stderr, "\nSimulating day %d...\n", Day + 1 ) ; StartIndex = Day * cSamplesPerDay ; RRSim1Day( RRday, RRnight, nT + StartIndex, nRR + StartIndex, nResp + StartIndex ) ; } } /* This function is called once, before generate() is called. See main() for a description of its arguments. */ void initialize(long seed, long tmax) { int Index=0 ; /* FILE *OutputFile ; */ /* tmax not used here, added to aviod compiler warning */ Index = tmax < Index ? tmax : Index; /* Allocate memory for the arrays */ /* Abandon hope if you can't */ nT = malloc( cMaxAlloc * sizeof( double ) ) ; if( nT != NULL ) nRR = malloc( cMaxAlloc * sizeof( double ) ) ; if( nRR != NULL ) nLFeffect = malloc( cMaxAlloc * sizeof( double ) ) ; if( nLFeffect != NULL ) nWhitenoise = malloc( cMaxAlloc * sizeof( double ) ) ; if( nWhitenoise != NULL ) nPinknoise = malloc( cMaxAlloc * sizeof( double ) ) ; if( nPinknoise != NULL ) nResp = malloc( cMaxAlloc * sizeof( double ) ) ; if( nResp == NULL ) { fprintf( stderr, "Couldn't allocate enough memory to continue." ) ; return ; } srand((unsigned int) seed ) ; RRSim( ) ; /* Write the intermediate data to the file RR.num in the working directory */ /* fprintf( stderr, "\nWriting to RR.num...\n" ) ; OutputFile = fopen( "RR.num", "wt" ) ; if( OutputFile == NULL ) return ; for( Index = 1 ; Index < cSamplesTotal ; Index ++ ) fprintf( OutputFile, "%8.3lf %8.3lf %8.3lf %8.3lf %8.3lf %8.3lf\n", nT[ Index ], nLFeffect[ Index ], nWhitenoise[ Index ], nPinknoise[ Index ], nResp[ Index ], nRR[ Index ] ) ; fclose( OutputFile ) ;*/ } /* This function is called once per RR interval. It should return the length of the next (simulated) RR interval in seconds. The example code generates samples of a noisy sine wave. */ /* float generate(void) { float rr; static float t; rr = 0.8 + 0.05*sin(t/5.5) + ((float)rand() - RAND_MAX)/(RAND_MAX*100.0); t += rr; return (rr); } */ /* Generate a series of RR intervals */ /* Start with the first element in the RR array */ /* This is the first RR interval - 'ThisRR' */ /* Following this heartbeat, time has moved on by 'ThisRR' seconds */ /* Which tells us where to find the next RR interval from */ float generate( void ) { static double CurrentTime = 0.0 ; double ThisRR ; int Index ; Index = TimeToIndex( CurrentTime ) + 1 ; if( Index < cSamplesTotal ) ThisRR = nRR[ Index ] ; CurrentTime += ThisRR ; return ThisRR ; } int main(int argc, char **argv) { float t = 0.0; /* sum of intervals since the beginning of the simulation, in seconds */ long ts = 0; /* t, in sample intervals */ long tsp = 0; /* previous value of ts */ long tmax = 24*60*60; /* 24 hours, in seconds */ long seed; /* a 32-bit random variable that can be used to initialize the generator */ long atol(); if (argc < 2) { fprintf(stderr, "usage: %s seed [tmax]\n", argv[0]); exit(1); } seed = atol(argv[1]); if (argc > 2) tmax = atol(argv[2]); initialize(seed, tmax); /*#if 0*/ while ((t += generate()) < tmax) { /* add RR interval to running time */ /* calculate and output a quantized RR interval */ ts = (long)(SPS*t + 0.5); printf("%5.3f\n", (ts - tsp)/((float)SPS)); tsp = ts; } /*#endif */ free( nT ) ; free( nRR ) ; free( nLFeffect ) ; free( nWhitenoise ) ; free( nPinknoise ) ; free( nResp ) ; exit(0); }