once the DC part of the signal is removed, the function can be convoluted with
itself to catch the period. Indeed, the convolution will feature peaks at each
multiple of the period. The FFT can be applied to compute the convolution.
fft = np.fft.rfft(L, norm="ortho")def abs2(x): return x.real**2 + x.imag**2selfconvol=np.fft.irfft(abs2(fft), norm="ortho")
The first output is not that good because the size of the image is not a
multiple of the period.
As noticed by Nils Werner, a window can be applied to limit the effect of
spectral leakage. As an alternative, the first crude estimate of the period
can be used to trunk the signal and the procedure can be repeated as I
answered in How do I scale an FFT-based cross-correlation such that its peak
is equal to Pearson’s rho.
From there, getting the period boils down to finding the first maximum. Here
is a way it could be done:
import numpy as npimport scipy.signalfrom matplotlib import pyplot as pltL = np.array([2.762, 2.762, 1.508, 2.758, 2.765, 2.765, 2.761, 1.507, 2.757, 2.757, 2.764, 2.764, 1.512, 2.76, 2.766, 2.766, 2.763, 1.51, 2.759, 2.759, 2.765, 2.765, 1.514, 2.761, 2.758, 2.758, 2.764, 1.513, 2.76, 2.76, 2.757, 2.757, 1.508, 2.763, 2.759, 2.759, 2.766, 1.517, 4.012])L = np.round(L, 1)# Remove DC component, as proposed by Nils WernerL -= np.mean(L)# Window signal#L *= scipy.signal.windows.hann(len(L))fft = np.fft.rfft(L, norm="ortho")def abs2(x): return x.real**2 + x.imag**2selfconvol=np.fft.irfft(abs2(fft), norm="ortho")selfconvol=selfconvol/selfconvol[0]plt.figure()plt.plot(selfconvol)plt.savefig('first.jpg')plt.show()# let's get a max, assuming a least 4 periods...multipleofperiod=np.argmax(selfconvol[1:len(L)/4])Ltrunk=L[0:(len(L)//multipleofperiod)*multipleofperiod]fft = np.fft.rfft(Ltrunk, norm="ortho")selfconvol=np.fft.irfft(abs2(fft), norm="ortho")selfconvol=selfconvol/selfconvol[0]plt.figure()plt.plot(selfconvol)plt.savefig('second.jpg')plt.show()#get ranges for first min, second maxfmax=np.max(selfconvol[1:len(Ltrunk)/4])fmin=np.min(selfconvol[1:len(Ltrunk)/4])xstartmin=1while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4: xstartmin=xstartmin+1xstartmax=xstartminwhile selfconvol[xstartmax]<fmin+0.7*(fmax-fmin) and xstartmax< len(Ltrunk)//4: xstartmax=xstartmax+1xstartmin=xstartmaxwhile selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4: xstartmin=xstartmin+1period=np.argmax(selfconvol[xstartmax:xstartmin])+xstartmaxprint "The period is ",period
