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