问题描述
我有一个周期信号，我想找到周期.
由于存在边界效应，因此我首先切出边界并通过查看第一个和最后一个最小值来保留N个周期.
然后，我计算FFT.
代码:
import numpy as np
from matplotlib import pyplot as plt
# The list of a periodic something
L = [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]
# Round because there is a slight variation around actually equals values: 2.762, 2.761 or 1.508, 1.507
L = [round(elt, 1) for elt in L]
minima = min(L)
min_id = L.index(minima)
start = L.index(minima)
stop = L[::-1].index(minima)
L = L[start:len(L)-stop]
fft = np.fft.fft(np.asarray(L))/len(L)
fft = fft[range(int(len(L)/2))]
plt.plot(abs(fft))
我知道列表的两点之间有多少时间(即采样频率，在这种情况下为190 Hz).我认为fft应该使我的峰值等于周期中点的数量，从而给我点的数量和周期.但是，这根本不是我观察到的输出:
我目前的猜测是，在0处的尖峰对应于我的信号的均值，大约在7处的这个小尖峰应该是我的周期(尽管重复模式仅包含5个点).
我做错了什么?谢谢！
一旦信号的DC部分被删除，该函数就可以自身进行卷积以捕获周期.实际上，卷积将在该周期的每个倍数处出现峰值. FFT可以用于计算卷积.
fft = np.fft.rfft(L, norm="ortho")
def abs2(x):
return x.real**2 + x.imag**2
selfconvol=np.fft.irfft(abs2(fft), norm="ortho")
第一个输出不是很好，因为图像的大小不是周期的倍数.
如Nils Werner所注意到的，可以应用一个窗口来限制频谱泄漏的影响.或者，可以使用该期间的第一个粗略估算来中继信号，并且可以重复该过程，就像我在.对周期信号进行自卷积从那里开始，将时间段归结为找到第一个最大值.这是一种可以完成的方法:
import numpy as np
import scipy.signal
from matplotlib import pyplot as plt
L = 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 Werner
L -= 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**2
selfconvol=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 max
fmax=np.max(selfconvol[1:len(Ltrunk)/4])
fmin=np.min(selfconvol[1:len(Ltrunk)/4])
xstartmin=1
while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4:
xstartmin=xstartmin+1
xstartmax=xstartmin
while selfconvol[xstartmax]<fmin+0.7*(fmax-fmin) and xstartmax< len(Ltrunk)//4:
xstartmax=xstartmax+1
xstartmin=xstartmax
while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4:
xstartmin=xstartmin+1
period=np.argmax(selfconvol[xstartmax:xstartmin])+xstartmax
print "The period is ",period
I have a periodic signal I would like to find the period.
Since there is border effect, I first cut out the border and keep N periods by looking at the first and last minima.
Then, I compute the FFT.
Code:
import numpy as np
from matplotlib import pyplot as plt
# The list of a periodic something
L = [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]
# Round because there is a slight variation around actually equals values: 2.762, 2.761 or 1.508, 1.507
L = [round(elt, 1) for elt in L]
minima = min(L)
min_id = L.index(minima)
start = L.index(minima)
stop = L[::-1].index(minima)
L = L[start:len(L)-stop]
fft = np.fft.fft(np.asarray(L))/len(L)
fft = fft[range(int(len(L)/2))]
plt.plot(abs(fft))
I know how much time I have between 2 points of my list (i.e. the sampling frequency, in this case 190 Hz). I thought that the fft should give me a spike at the value corresponding to the number of point in a period, , thus giving me the number of point and the period.Yet, that is not at all the output I observed:
My current guess is that the spike at 0 corresponds to the mean of my signal and that this little spike around 7 should have been my period (although, the repeating pattern only includes 5 points).
What am I doing wrong? Thanks!
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**2
selfconvol=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 np
import scipy.signal
from matplotlib import pyplot as plt
L = 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 Werner
L -= 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**2
selfconvol=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 max
fmax=np.max(selfconvol[1:len(Ltrunk)/4])
fmin=np.min(selfconvol[1:len(Ltrunk)/4])
xstartmin=1
while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4:
xstartmin=xstartmin+1
xstartmax=xstartmin
while selfconvol[xstartmax]<fmin+0.7*(fmax-fmin) and xstartmax< len(Ltrunk)//4:
xstartmax=xstartmax+1
xstartmin=xstartmax
while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4:
xstartmin=xstartmin+1
period=np.argmax(selfconvol[xstartmax:xstartmin])+xstartmax
print "The period is ",period
这篇关于从FFT中找出信号的周期的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！