Loop matrix

#pragma rtGlobals=3 // Use modern global access method and strict wave access. //------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ //----------------------------------------------------------------------------------------- SELECT BC waves to be sorted ----------------------------------------------------------------------------------------- //------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ Function WavesSorted_BC_Test() String date_utc_time //Name of the Date&Time wave that we want to analyze string BC1 string BC2 prompt date_utc_time, "Date&Time wave to be sorted", popup WaveList("*", ";", "") prompt BC1, "BC1 wave to be sorted", popup WaveList("*", ";", "") prompt BC2, "BC2 wave to be sorted", popup WaveList("*", ";", "") DoPrompt "Waves sorted", date_utc_time, BC1, BC2 SortedByTimeOfDay_BC_Test($date_utc_time, $BC1, $BC2) End //-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- //-------------------- Function to sort input Date&Time wave by the time of day. NaNs (blanks) go to the end -------------------- //-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- Function/S SortedByTimeOfDay_BC_Test(date_utc_time, BC1, BC2) wave date_utc_time //Date&Time wave to be sorted wave BC1, BC2 //BC waves to be sorted variable i duplicate/O date_utc_time, TimeOfDay_Sorted duplicate/O BC1, BC1_Sorted duplicate/O BC2, BC2_Sorted TimeOfDay_Sorted = mod(date_utc_time, 24*60*60) // TimeOfDay now contains just time of day (remove date) Sort {TimeOfDay_Sorted,date_utc_time}, TimeOfDay_Sorted, BC1_Sorted, BC2_Sorted // Sort by TimeOfDay, then by date CalcConcStats_BC_Test(TimeOfDay_Sorted, BC1_Sorted, BC2_Sorted, 60) End //---------------------------------------------------------------------------------------------------------------- //-------------------- Function to estimate the statistics of the waves -------------------- //---------------------------------------------------------------------------------------------------------------- Function CalcConcStats_BC_Test(sortedTimeOfDay, sortedBC1, sortedBC2, averagingIntervalInMinutes) wave sortedTimeOfDay // Time-of-day wave wave sortedBC1, sortedBC2 // Concentrations wave variable averagingIntervalInMinutes // Interval over which to average in minutes; e.g., 60 for 1 hour variable numBins = (24 * 60) / averagingIntervalInMinutes variable averagingIntervalInSeconds = averagingIntervalInMinutes * 60 //Number of columns of Statistical tables variable numStatsColumns = 12 // Variables to find the first and last point of each bin variable numInputPoints = numpnts(sortedTimeOfDay) variable i, bin, timeOfDay Variable BC1, BC2 // To find the median and standard deviation for each bin variable first variable last variable count variable median //-------------------------------------------------------------------------------------------// //-------------------- Make output waves of BC waves --------------------// //-------------------------------------------------------------------------------------------// //-------------------------------------------------- //-------------------- BC1 -------------------- //-------------------------------------------------- Make /O /D /N=(numBins,numStatsColumns) BC1Stats = 0 // 2D wave SetDimLabel 1, 0, Count, BC1Stats // Column 0 is count of concentrations in bin SetDimLabel 1, 1, First, BC1Stats // Column 1 is first input point in bin SetDimLabel 1, 2, Last, BC1Stats // Column 2 is last input point in bin SetDimLabel 1, 3, Min, BC1Stats // Column 3 is median concentration for bin SetDimLabel 1, 4, Max, BC1Stats // Column 4 is median concentration for bin SetDimLabel 1, 5, Mean, BC1Stats // Column 5 is mean concentration for bin SetDimLabel 1, 6, StdDev, BC1Stats // Column 6 is standard deviation for bin SetDimLabel 1, 7, Median, BC1Stats // Column 7 is median concentration for bin SetDimLabel 1, 8, Q25, BC1Stats // Column 8 is Quartil25 for bin SetDimLabel 1, 9, Q75, BC1Stats // Column 9 is Quartil75 for bin SetDimLabel 1, 10, upperInnerFence1, BC1Stats // Column 10 is for bin SetDimLabel 1, 11, lowerInnerFence1, BC1Stats // Column 11 is for bin // Set X scaling such that the statistics for each bin fall in the middle of bin when graphed SetScale/P x, averagingIntervalInSeconds-3600, averagingIntervalInSeconds, "dat", BC1Stats // Preset columns of output BC1Stats[][%first] = -1 BC1Stats[][%last] = -1 // Find the first and last point of each bin, and count for (i=0; i<numInputPoints; i+=1) // Bin and sum up the concentrations timeOfDay = sortedTimeOfDay[i] bin = trunc(timeOfDay / averagingIntervalInSeconds) if (BC1Stats[bin][%first] == -1) BC1Stats[bin][%first] = i endif BC1Stats[bin][%last] = i if (numtype(BC1) == 0) // Exclude NaNs and INFs BC1Stats[bin][%count] += 1 // Increment bin count endif count = BC1Stats[bin][%count] endfor // Find the minimum for each bin BC1Stats[][%min] = WaveMin(sortedBC1, BC1Stats[p][%first], BC1Stats[p][%last]) // Find the maximum for each bin BC1Stats[][%max] = WaveMax(sortedBC1, BC1Stats[p][%first], BC1Stats[p][%last]) // Find the mean, median and standard deviation for each bin for(bin=0; bin<numBins; bin+=1) first = BC1Stats[bin][%first] // First input point for this bin last = BC1Stats[bin][%last] // Last input point for this bin Duplicate /R=[first,last] /FREE sortedBC1, temp1 // Needed because StatsMedian works on the whole wave WaveTransform zapNaNs, temp1 // Needed because StatsMedian does not support NaNs median = StatsMedian(temp1) BC1Stats[bin][%median] = median // WaveStats /M=2 /Q temp1 WaveStats/Q temp1 BC1Stats[bin][%mean] = V_avg BC1Stats[bin][%stddev] = V_sdev StatsQuantiles/Q temp1 BC1Stats[bin][%Q25] = V_Q25 BC1Stats[bin][%Q75] = V_Q75 Variable H1=V_Q75 - V_Q25 Variable step1=1.5*H1 Variable upperInnerFence1=V_Q75 + step1 Variable upperOuterFence1=V_Q75 + 2*step1 Variable lowerInnerFence1=V_Q25 - step1 Variable lowerOuterFence1=V_Q25 - 2*step1 BC1Stats[bin][%upperInnerFence1] = upperInnerFence1 BC1Stats[bin][%lowerInnerFence1] = lowerInnerFence1 endfor //-------------------------------------------------- //-------------------- BC2 -------------------- //-------------------------------------------------- Make /O /D /N=(numBins,numStatsColumns) BC2Stats = 0 // 2D wave SetDimLabel 1, 0, Count, BC2Stats // Column 0 is count of concentrations in bin SetDimLabel 1, 1, First, BC2Stats // Column 1 is first input point in bin SetDimLabel 1, 2, Last, BC2Stats // Column 2 is last input point in bin SetDimLabel 1, 3, Min, BC2Stats // Column 3 is median concentration for bin SetDimLabel 1, 4, Max, BC2Stats // Column 4 is median concentration for bin SetDimLabel 1, 5, Mean, BC2Stats // Column 5 is mean concentration for bin SetDimLabel 1, 6, StdDev, BC2Stats // Column 6 is standard deviation for bin SetDimLabel 1, 7, Median, BC2Stats // Column 7 is median concentration for bin SetDimLabel 1, 8, Q25, BC2Stats // Column 8 is Quartil25 for bin SetDimLabel 1, 9, Q75, BC2Stats // Column 9 is Quartil75 for bin SetDimLabel 1, 10, upperInnerFence2, BC2Stats // Column 10 is for bin SetDimLabel 1, 11, lowerInnerFence2, BC2Stats // Column 10 is for bin // Set X scaling such that the statistics for each bin fall in the middle of bin when graphed SetScale/P x, averagingIntervalInSeconds-3600, averagingIntervalInSeconds, "dat", BC2Stats // Preset columns of output BC2Stats[][%first] = -1 BC2Stats[][%last] = -1 // Find the first and last point of each bin for (i=0; i<numInputPoints; i+=1) // Bin and sum up the concentrations timeOfDay = sortedTimeOfDay[i] bin = trunc(timeOfDay / averagingIntervalInSeconds) if (BC2Stats[bin][%first] == -1) BC2Stats[bin][%first] = i endif BC2Stats[bin][%last] = i if (numtype(BC2) == 0) // Exclude NaNs and INFs BC2Stats[bin][%count] += 1 // Increment bin count endif count = BC2Stats[bin][%count] endfor // Find the minimum for each bin BC2Stats[][%min] = WaveMin(sortedBC2, BC2Stats[p][%first], BC2Stats[p][%last]) // Find the maximum for each bin BC2Stats[][%max] = WaveMax(sortedBC2, BC2Stats[p][%first], BC2Stats[p][%last]) // Find the median and standard deviation for each bin for(bin=0; bin<numBins; bin+=1) first = BC2Stats[bin][%first] // First input point for this bin last = BC2Stats[bin][%last] // Last input point for this bin Duplicate /R=[first,last] /FREE sortedBC2, temp2 // Needed because StatsMedian works on the whole wave WaveTransform zapNaNs, temp2 // Needed because StatsMedian does not support NaNs median = StatsMedian(temp2) BC2Stats[bin][%median] = median // WaveStats /M=2 /Q temp2 WaveStats/Q temp2 BC2Stats[bin][%mean] = V_avg BC2Stats[bin][%stddev] = V_sdev StatsQuantiles/Q temp1 BC2Stats[bin][%Q25] = V_Q25 BC2Stats[bin][%Q75] = V_Q75 Variable H2=V_Q75 - V_Q25 Variable step2=1.5*H2 Variable upperInnerFence2=V_Q75 + step2 Variable upperOuterFence2=V_Q75 + 2*step2 Variable lowerInnerFence2=V_Q25 - step2 Variable lowerOuterFence2=V_Q25 - 2*step2 BC2Stats[bin][%upperInnerFence2] = upperInnerFence2 BC2Stats[bin][%lowerInnerFence2] = lowerInnerFence2 endfor End

chozo

So you want to avoid writing the same code over and over again? Is there a specific meaning for CalcConcStats_BC_Test to take two waves BC1 and BC2 (I guess not from the code)? If not just rewrite the code to use dynamic wave names instead of hard-coded ones (i.e., BC1Stats etc. ... never a good idea but for to simplest cases), and then use a function to loop though a wave list repeatedly calling your stats function in the process. Here this should give you some idea (note the use of the '$', and you should probably read more about it by executing DisplayHelpTopic "Converting a String into a Reference Using $"):

Function DynamicExampe(inwave)
    wave inwave
    String statsname = nameofwave(inwave)+"_stats"
    Make/D/O/N=(round(DimSize(inwave,0)/2)) $statsname
    SetScale/P x,DimOffset(inwave,0),DimDelta(inwave,0)*2, $statsname
    Wave workwave = $statsname
    workwave = inwave[p*2] 
End

Btw. this function downsamples the incoming data to halve the points in the process.

July 2, 2015 at 04:04 am - Permalink

thomas_braun

For resampling data I prefer Resample.

July 2, 2015 at 04:12 am - Permalink

chozo

Just to clear things up: Of course the function has no real world application besides demonstrating the use of '$'. I thought it might be nice to have at least some effect beyond creating an empty wave. Btw I use Interpolate2 for 'resampling', because I rather dial in the output points.

July 2, 2015 at 05:51 am - Permalink

thomas_braun

@chozo: Sure, I didn't want to imply that this is how you would implement resampling.

July 2, 2015 at 03:20 pm - Permalink

Marta

Thank you chozo and thomas_braun for your help.

The waves come from an Aethalometer instrument, so I have seven black carbon waves (BC1, BC2, ..., BC7). The function CalcConcStats_BC_Test() calculates the statistics of each wave.
I am going to try what you wrote me, I hope it works!

Thanks again!

July 3, 2015 at 12:19 am - Permalink