See Sethares and Staley, “Periodicity Transforms”  IEEE Trans. Signal Processing, 1998.
projectP[s,p] finds the closest p-periodic vector to the vector s by projecting s onto the periodic subspace P_p
periodNorm[x] calculates the norm of the vector x using the inner product defined on the p-periodic subspaces
argMax[ ] and argMin[ ] find the index of the Max and Min elements of a sequence
projectRangeP[r,range] carries out projections of the sequence r onto all the p in the range
The four basic Transforms are:
small2Large[ ]
mBest[ ] and the imporved version mBestGam[ ]
bestFrequency[ ]
bestCorrelation[ ]

projectP[s_,p_]:=Module[{baseElem=ConstantArray[0,Length[s]],ones=ConstantArray[1,Length[s]],ind,u},Do[ind=Select[PadRight[u=UnitVector[p,shift],Length[s],u] Range[Length[s]],#>0&];
baseElem[[ind]]=Mean[s[[ind]]] ones[[ind]];,{shift,1,p}];
baseElem];
projectRangePGam[r_,range_,gam_]:=Module[{theP},Table[{j,periodNorm[theP=projectP[r,j]]/(j^gam),theP},{j,range}]];
periodNorm[x_]:=Norm[x]/Sqrt[Length[x]];
argMax[x_]:=Flatten[Position[x,Max[x]]];
argMin[x_]:=Flatten[Position[x,Min[x]]];
small2Large[s_,thresh_,numCheck_]:=Module[{s2lList={},ss=s,sNorm=periodNorm[s],baseElem,normImp,r,ind},Do[baseElem=projectP[ss,period];
r=ss-baseElem;
normImp=(periodNorm[ss]-periodNorm[r])/sNorm;
If[normImp>thresh,ss=r;
AppendTo[s2lList,{period,normImp,baseElem}];];,{period,1,numCheck}];
ind=Reverse[Ordering[N[Transpose[s2lList][[2]]]]];
Transpose[s2lList[[ind]]]];
mBestGam[x_,m_,len_,gam_]:=Module[{r=x,mBestList=ConstantArray[{{},{},{}},m],biggest,bigP,listChange,minMBestInd,thisP,allProjections,maxNorm,maxP,xBigQ,xSmallQ,nBigQ,nSmallQ,minQ,ptTwice,i,ind},Do[allProjections=projectRangePGam[r,Range[len],gam];
biggest=First[argMax[allProjections[[All,2]]]];
mBestList[[i]]={biggest,allProjections[[biggest,2]],bigP=allProjections[[biggest,3]]};r=r-bigP;,{i,1,m}];
listChange=1;
While[listChange==1,i=1;
While[i≤m,listChange=0;
minMBestInd=First[argMin[mBestList[[All,2]]]];
thisP=mBestList[[i,1]];
If[Length[Divisors[thisP]]>1,factors=Most[Rest[Divisors[thisP]]],factors={}];
If[factors≠{},allProjections=projectRangePGam[xqi=mBestList[[i,3]],factors,gam];
maxNorm=First[argMax[allProjections[[All,2]]]];
maxP=allProjections[[maxNorm,1]];
xBigQ=allProjections[[maxNorm,3]];
xSmallQ=xqi-xBigQ;
nBigQ=periodNorm[xBigQ]/maxP^gam;
nSmallQ=periodNorm[xSmallQ]/maxP^gam;
minQ=Min[mBestList[[All,2]]];
ptTwice=Position[mBestList[[All,1]],maxP];
If[nSmallQ+nBigQ>minQ+mBestList[[i,2]]&&nSmallQ>minQ&&nBigQ>minQ&&ptTwice=={},listChange=1;
mBestList[[minMBestInd]]={maxP,nBigQ,xBigQ};
mBestList[[i]]={thisP,nSmallQ,xSmallQ};];];
i=i+1;];];
mBestList[[All,2]]=Length[x]^gam mBestList[[All,2]]/periodNorm[x];
ind=Reverse[Ordering[Transpose[mBestList][[2]]]];
Transpose[mBestList[[ind]]]];
bestFrequency[x_,m_,harm_]:=Module[{mFreqList=ConstantArray[{{},{},{}},m],r=x,lenX=Length[x],ssf,i=1,fft,absFFT,maxF,p,bas,ind,rTest,indMin,pC},ssf=Range[0,lenX/2,1]/lenX;ssf[[1]]=1;
ssf[[2]]=2/lenX;
While[i≤m,fft=Fourier[r,FourierParameters->{1,-1}];
absFFT=Abs[fft[[1;;Round[lenX/2]]]];
maxF=argMax[absFFT];
pC=First[Round[1/ssf[[maxF]]]];
p=Sort[Flatten[{{pC Range[harm]-1},{pC Range[harm]},{pC Range[harm]+1}}]];
bas=Table[projectP[r,p[[j]]],{j,1,Length[p]}];
rTest=Table[r-bas[[j]],{j,1,Length[p]}];
indMin=First[argMin[Table[periodNorm[rTest[[i]]],{i,1,Length[p]}]]];
mFreqList[[i]]={p[[indMin]],periodNorm[bas[[indMin]]],bas[[indMin]]};
r=rTest[[indMin]];
i=i+1;];
ind=Reverse[Ordering[Transpose[mFreqList][[2]]]];
Transpose[mFreqList[[ind]]]];
bestCorrelation[x_,ratPower_,numIter_,longPer_]:=Module[{lenX=Length[x],bestCorrList={},numt=0,oldNorm,s,ind,xNorm=periodNorm[x]},oldNorm=xNorm;
s=x-Mean[x];
While[numt<numIter,numt=numt+1;
allCorr=Table[Abs[Total[s[[Range[shift,lenX,k]]]]],{k,1,longPer},{shift,1,k}];
{maxP,indS}=First[Position[allCorr,Max[allCorr]]];
bas=projectP[s,maxP];
s=s-bas;
normTest=(oldNorm-periodNorm[s])/xNorm;
If[normTest>ratPower,AppendTo[bestCorrList,{maxP,normTest,bas}];
oldNorm=periodNorm[s];,numt=numIter];];
ind=Reverse[Ordering[N[Transpose[bestCorrList][[2]]]]];
Transpose[bestCorrList[[ind]]]];
(*alternative definitions allow use of default values*)

small2Large[s_,thresh_]:=small2Large[s,thresh,Length[s]/5];
small2Large[s_]:=small2Large[s,0.1,Length[s]/5];
mBestGam[x_,m_,len_]:=mBestGam[x,m,len,0.5];
mBestGam[x_,m_]:=mBestGam[x,m,Length[x]/5,0.5];
mBestGam[x_]:=mBestGam[x,5,Length[x]/5,0.5];
mBest[x_,m_,len_]:=mBestGam[x,m,len,0.];
mBest[x_,m_]:=mBestGam[x,m,Length[x]/5,0.];
mBest[x_]:=mBest[x,5,Length[x]/5,0.];
bestFrequency[x_,m_]:=bestFrequency[x,m,1];
bestFrequency[x_]:=bestFrequency[x,5,1];
bestCorrelation[x_,ratPower_,numIter_]:=bestCorrelation[x,ratPower,numIter,Length[x]/5];
bestCorrelation[x_,ratPower_]:=bestCorrelation[x,ratPower,10,Length[x]/5];
bestCorrelation[x_]:=bestCorrelation[x,0.01,10,Length[x]/5];