Chi-Square periodogram

#pragma rtGlobals=1		// Use modern global access method.



//////////////////////////
/////	Chi-Square periodogram
////
////	chi_square_periodogram(w1, size, start_val, end_val)	
////	w1: wave to be analyzed. Its unit must be "min"
////	size : resolution of periodogram (in "min")
////	start_val, end_val:periodicity range from start_val to end_val (in "min") 
////						to be focused on for analysis
////
////	Result wave is made as "Qp"
////	The stronger is the rhythmiciy in a data set, the higher is the value of Qp.
////
////	References
////	1. Refinetti, Journal of Theoretical Biology 277, 571-581, 2004 
////		Qp=K*N*Sigma[h=1->P](Mh-M)^2 / Sigma[i=1->N](Xi-N)^2
////
////	2. Sokolve & Bushell, Journal of Theoretical Biology 72, 131-160, 1978
////	
////	made on 23 Jan 2009 by Y. Ootsuka
////	
////
Function/S chi_square_periodogram(w1, size, start_val, end_val)	
	wave w1
	variable size, start_val, end_val
	variable Period		// period  e.g. 10min, 20min, .... ,100min,..., 200min
	variable M			// mean of all N value
	variable N			// a number of total data points collected
	variable K			// a number of blocks. 
	variable P			//NumofBins			
	variable h,i,j, blockNum
	Variable Qp_denominator, Qp_nominator
	variable SumMh_M
	
	Make /o/n=(round((end_val-start_val)/size)) Qp //	probability distribution of chi-square with P-1 degress of freedom
	wavestats /Q w1
	Qp_denominator =V_sdev^2*(V_npnts-1)
	M=V_avg								//	mean of all values of  data 
	N=numpnts(w1)						//	TotalNumData
	
	Do	
		Period=start_val+size*j
		K=round(N*deltax(w1)/Period)	//	TotalNumofBlocks
		P=round(Period/deltax(w1))		//	NumofBins per period
		Make /o/n=(P) Mh	
		
		For (h=0;h<P;h+=1)
			For (blockNum=0;blockNum<K;blockNum+=1) //  cal Mh[h] from K blocks of period P
				if (blockNum==0)
				Mh[h]=w1[blockNum*P+h]
				else
				Mh[h]+=w1[blockNum*P+h]
				endif
			endfor
			Mh[h]/=(K-1)					//  Mean of Mh[h] from K blocks of period P
			if (h==0)
				SumMh_M=(Mh[h]-M)^2
			else
				SumMh_M+=(Mh[h]-M)^2
			endif
		endfor
		
		Qp_nominator=K*N*(SumMh_M-M)
		Qp[j]=Qp_nominator/Qp_denominator
		j+=1
	while(Period<=end_val || j<numpnts(Qp))
	
	SetScale/P x start_val,size,"min", Qp
	return nameofwave(Qp)
end