time average waves

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3     // Use modern global access method and strict wave access.
// This function will take a pair of waves, are return a new pair of waves with
// the desired periodicity

function periodicMean(timeWave, valueWave, key, [mult,suffix])
    // timeWave is a real wave representing time in secoinds since the epoch
    // is a wave with y values of the same length as the time wave
    // key is sec, min, hour, day, and represents the time unit of hte output waves
    // multipler is an optiona multipler of the time unit
    // suffix is a user provided suffix for the output waves.   optioal if there is no muliplier
   
    // output is a pair of waves averaged to the time unit desired

    Wave timeWave // assumes time in seconds (since epoch 1/1/1904)
    Wave valueWave
    String key
    Variable mult
    String suffix
   
    variable multiplier

   
   
    if(ParamIsDefault(mult))
        multiplier = 1
    else
        multiplier = mult
        if (ParamIsDefault(suffix))
            print "Output wave needs a unique suffix parameter of the form: suffix=<string>"
            return -1
        endif
    endif
   
    if(ParamIsDefault(suffix))
        suffix = key
    endif
   
    String periodDict = "sec:1;min:60;hour:3600;day:86400"
    String newXName, newYName

    if (strsearch(periodDict, key, 0) == -1)
        print "Invalid period string, allowed values are sec, min, hour and day"
        return -1
    endif
   
    if (numpnts(timeWave) != numpnts(valueWave))
        print "wavelength mismatch"
        return -1
    endif
   
    // get the total duration of the time wave
    Variable startTime = timeWave[0]
    Variable endTime = timeWave[INF]
    Variable duration = endTime - StartTime
    Variable numSeconds = NumberByKey(key, periodDict)*multiplier
    Variable numPeriods = duration/numSeconds
    if (numPeriods < 1)
        print "the period is too large for this data set"
        return -1
    endif  
   
    Variable wholePeriods = ceil(numPeriods)
    Make /O/N=(wholePeriods)/D xWave, yWave

    variable i
    variable averageTime
    variable intervalStart, intervalEnd
    variable numpntsperinterval = numpnts(timeWave)/wholePeriods
    for (i=0; i < wholePeriods; i+=1 )
        intervalStart = i*numpntsperinterval
        intervalEnd = intervalStart + numpntsperinterval
        if (intervalEnd > numpnts(timeWave)-1)
            intervalEnd = numpnts(timeWave) - 1
        endif
       
        averageTime = (timeWave[intervalEnd] + timeWave[intervalStart])/2
        WaveStats/M=1/Q/R=[intervalStart,intervalEnd] valueWave
        xWave[i] = averageTime
        yWave[i] = V_avg
    endfor
   
    newXName = nameOfWave(timeWave) + "_" + suffix
    newYName = nameOfWave(valueWave) + "_" + suffix
    Rename xWave, $newXName
    Rename yWave, $newYName
end