Polar averaging and errors

Averaging waves with polar data (list of angles) will give incorrect results in many cases (see image below). This function does a vector summation of X and Y components to determine the average angle, the standard deviation, or the standard error.

Function polarAvgErr(inWave,opStr,degrad)
	Wave inWave //single wave with list of angles
	
	//opStr = "avg", "sdev", or "sem"
	//degrad = "deg" or "rad" depending on the input data type 
	String opStr,degrad
	Variable i,size,avgX,avgY,avg,sdev,sem,errX,errY,value,multiplier
	
	size = DimSize(inWave,0)
	Make/FREE/N=(size) xTotal,yTotal
		
	//input is degrees or radians
	strswitch(degrad)
		case "deg":
			xTotal = cos(inWave * pi/180)
			yTotal = sin(inWave * pi/180)
			multiplier = 180/pi
			break
		case "rad":
			xTotal = cos(inWave)
			yTotal = sin(inWave)
			multiplier = 1
			break
	endswitch

	//polar average
	avgX = mean(xTotal)
	avgY = mean(yTotal)
	avg = atan2(avgY,avgX)*multiplier
	
	If(avg < 0)
		avg += 360
	EndIf
	
	//polar standard deviation
	For(i=0;i<size;i+=1)		
		errX += (xTotal[i] - avgX)^2/(size-1)
		errY += (yTotal[i] - avgY)^2/(size-1)
	EndFor
	
	errX = sqrt(errX)		
	errY = sqrt(errY)
	sdev = atan2(errY,errX)*multiplier
			
	//polar standard error of the mean
	sem = sdev/sqrt(size)
	
	//select return parameter
	strswitch(opStr)
		case "avg":
			value = avg
			break
		case "sdev":
			value = sdev
			break
		case "sem":
			value = sem
			break
	endswitch
	
	return value
End

 

You may want to look at StatsCircularMoments as it may save you some coding.

I find the description and polar plots a little confusing. You imply that the only input is a 1-D wave of polar angles (or you use only the angular part of a {theta, r} wave), but the data set has multiple 'r' values. Are you selecting a constant-r set of values from your presented data?

Another approach for full 2-D (x,y) or (r,theta)-->(x,y) data might be to first perform a linear ODR fit to an (x,y) line constrained to pass through the origin to find the best-fit average angle. Then find the difference angle for each point to get the set of d_theta's for the standard deviation calculation. I note that StatsCircularMoments seems to apply only to 1-D data at a fixed radius.

Forum

Support

Gallery

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More