在算法信号中找到周期性 [英] Finding periodicity in an algorithmic signal
问题描述
在测试关于以下递归关系的猜想
In testing a conjecture about the following recursive relation
，
它声称数字序列具有某种周期性，我编写了一个python程序来计算序列并将其打印在表格中.
which claims a periodicity of some kind for the sequence of numbers, I wrote a python program which computes the sequences and prints them in a table.
1 # Consider the recursive relation x_{i+1} = p-1 - (p*i-1 mod x_i)
2 # with p prime and x_0 = 1. What is the shortest period of the
3 # sequence?
4
5 from __future__ import print_function
6 import numpy as np
7 from matplotlib import pyplot as plt
8
9 # The length of the sequences.
10 seq_length = 100
11
12 upperbound_primes = 30
13
14 # Computing a list of prime numbers up to n
15 def primes(n):
16 sieve = [True] * n
17 for i in xrange(3,int(n**0.5)+1,2):
18 if sieve[i]:
19 sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
20 return [2] + [i for i in xrange(3,n,2) if sieve[i]]
21
22 # The list of prime numbers up to upperbound_primes
23 p = primes(upperbound_primes)
24
25 # The amount of primes numbers
26 no_primes = len(p)
27
28 # Generate the sequence for the prime number p
29 def sequence(p):
30 x = np.empty(seq_length)
31 x[0] = 1
32 for i in range(1,seq_length):
33 x[i] = p - 1 - (p * (i-1) - 1) % x[i-1]
34 return x
35
36 # List with the sequences.
37 seq = [sequence(i) for i in p]
38 """
39 # Print the sequences in a table where the upper row
40 # indicates the prime numbers.
41 for i in range(seq_length):
42 if not i:
43 for n in p:
44 print('\t',n,end='')
45 print('')
46 print(i+1,'\t',end='')
47 for j in range(no_primes):
48 print(seq[j][i],end='\t')
49 print('\n',end='')
50 """
51 def autocor(x):
52 result = np.correlate(x,x,mode='full')
53 return result[result.size/2:]
54
55
56 fig = plt.figure('Finding period in the sequences')
57 k = 0
58 for s in seq:
59 k = k + 1
60 fig.add_subplot(no_primes,1,k)
61 plt.title("Prime number %d" % p[k-1])
62 plt.plot(autocor(s))
63 plt.show()
64
现在，我想研究我计算出的这些序列中的周期性.在网上环顾四周后，我发现自己似乎有两种选择:
Now I want to investigate periodicities in these sequences that I computed. After looking around on the net I found myself two options it seems:
- 对数据进行自相关并寻找第一个峰.这应该给出周期的近似值.
- 对数据执行FFT.这显示了数字的频率.我看不到它如何提供有关数字序列的周期性的任何有用信息.
The last lines show my attempt of using autocorrelation, inspired by the accepted answer of How can I use numpy.correlate to do autocorrelation?.
它给出了以下情节
显然，我们看到所有质数的数字降序.
Clearly we see a descending sequence of numbers for all the primes.
使用以下简单的python代码片段在sin函数上测试相同的方法时
When testing the same method on a sin function with the following simplyfied python-code snippet
1 # Testing the autocorrelation of numpy
2
3 import numpy as np
4 from matplotlib import pyplot as plt
5
6 num_samples = 1000
7 t = np.arange(num_samples)
8 dt = 0.1
9
10 def autocor(x):
11 result = np.correlate(x,x,mode='full')
12 return result[result.size/2:]
13
14 def f(x):
15 return [np.sin(i * 2 * np.pi * dt) for i in range(num_samples)]
16
17 plt.plot(autocor(f(t)))
18 plt.show()
我得到了类似的结果，它给出了正弦函数的以下图
I get a similar result, it giving the following plot for the sine function
例如，在正弦函数情况下，如何读取周期性?
How could I read off the periodicity in the sine-function case, for example?
无论如何，我不了解自相关导致峰值的机制，这些峰值给出了信号的周期性信息.有人可以详细说明吗?在这种情况下，您如何正确使用自相关?
Anyhow, I do not understand the mechanism of the autocorrelation leading to peaks that give information of the periodicity of a signal. Can someone elaborate on that? How do you properly use autocorrelation in this context?
在实现自相关时，我还做错了什么?
Also what am I doing wrong in my implementation of the autocorrelation?
欢迎提出有关确定一系列数字中的周期性的其他方法的建议.
Suggestions on alternative methods of determining periodicity in a sequence of numbers are welcome.
推荐答案
这里有很多问题，所以我将开始描述自相关如何产生以"3"为例的周期，即您的第二个第一张图片的子图.
There are quite a few questions here, so I'll start be describing how an autocorrelation produces the period from the case of "3", ie, your second sub-plot of the first image.
对于素数3，序列是(不太一致的开始之后)1,2,1,2,1,2,1,2,....为了计算自相关，基本上将数组相对于其自身进行平移，将所有对齐的元素相乘，然后将所有这些结果相加.因此，在一些测试用例中，它看起来像这样，其中A是自相关:
For prime number 3, the sequence is (after a less consistent start) 1,2,1,2,1,2,1,2,.... To calculate the autocorrelation, the array is basically translated relative to itself, all the elements that align are multiplied, and all of these results are added. So it looks something like this, for a few test cases, where A is the autocorrelation:
0 1 2 3 4 5 6 7 ... 43 44 45 46 47 48 49 # indices 0
1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 # values 0
1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 # values 1
0 1 2 3 4 5 6 7 ... 43 44 45 46 47 48 49 # indices 1
1 4 1 4 1 4 1 4 4 1 4 1 4 1 4 # products
# above result A[0] = 5*25 5=1+4 25=# of pairs # A[0] = 125
0 1 2 3 4 5 6 7 ... 43 44 45 46 47 48 49 # indices 0
1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 # values 0
1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 # values 1
0 1 2 3 4 5 6 7 ... 43 44 45 46 47 48 49 # indices 1
2 2 2 2 2 2 2 2 2 2 2 2 2 2 # products
# above result A[1] = 4*24 4=2+2 24=# of pairs # A[1] = 96
0 1 2 3 4 5 6 7 ... 43 44 45 46 47 48 49 # indices 0
1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 # values 0
1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 # values 1
0 1 2 3 4 5 6 7 ... 43 44 45 46 47 48 49 # indices 1
1 4 1 4 1 4 1 4 4 1 4 1 4 # products
# above result A[2] = 5*23 5=4+1 23=# of pairs # A[2] = 115
上面有3条带回家的消息: 1..当相似元素排成一行并相乘时，自相关A具有较大的值，此处每隔一步. 2..自相关的索引对应于相对偏移. 3 .当在整个数组上进行自相关时，如此处所示，由于每次相加后每次相加产生的点数都会减少，因此总会有下降的趋势.
There are three take-home messages from the above: 1. the autocorrelation, A, has larger value when like elements are lined up and multiplied, here at every other step. 2. The index of the autocorrelation corresponds to the relative shift. 3. When doing the autocorrelation over the full arrays, as shown here, there's always a downward ramp since the number of points added together to produce the value are reduced at each successive shift.
因此，您可以在这里看到为什么"Prime number 3"中的图表会定期出现20％的颠簸:因为相加后的项在对齐时为1 + 4，而在未对齐时为2 + 2，即5对4.这就是您在阅读时段时要寻找的障碍.也就是说，这里显示的期间是2，因为那是您第一次碰撞的索引. (另外，请注意，在上面，我仅以对数的形式进行计算，以了解这种已知的周期如何导致您在自相关中看到的结果，也就是说，人们通常不希望考虑对数)
So here you can see why there's a periodic 20% bump in your graph from "Prime number 3": because the terms that are summed are 1+4 when they are aligned, vs 2+2 when they aren't, ie, 5 vs 4. It's this bump that you're looking for in reading off the period. That is, here is shows that the period is 2, since that is the index of your first bump. (Also, note btw, in the above I only do the calculation as number of pairs to see how this known periodicity leads to the result you see in the autocorrelation, that is, one doesn't in general want to think of number of pairs.)
在这些计算中，如果您在进行自相关之前先减去均值，则相对于基准的凹凸值将增加.如果使用端部经过修剪的数组进行计算，则可以删除坡道，因此总是有相同的重叠.这通常是有道理的，因为通常人们正在寻找一个比整个样本短得多的波长的周期性(因为要确定一个好的振荡周期需要很多次振荡).
In these calculations, the values of the bump relative to the base will be increased if you first subtract the mean before doing the autocorrelation. The ramp can be removed if you do the calculation using arrays with trimmed ends, so there's always the same overlap; this usually makes sense since usually one is looking for a periodicity of much shorter wavelength than the full sample (because it takes many oscillations to define a good period of oscillation).
对于正弦波的自相关，基本答案是该周期显示为第一个凸起.除了应用了时间轴外，我重做了该图.在这些方面，使用实时轴始终是最清晰的，因此我对您的代码进行了一些更改以使其包含在内. (此外，我用适当的矢量化numpy表达式代替了列表推导式来计算正弦波，但这并不重要.而且我还明确定义了f(x)的频率，只是为了更清楚地说明发生了什么事情-造成混淆的隐式频率为1.)
For the autocorrelation of the sine wave, the basic answer is that the period is shown as the first bump. I redid the plot except with the time axis applied. It's always clearest in these things to use a real time axis, so I changed your code a bit to include that. (Also, I replaced the list comprehension with a proper vectorized numpy expression for calculating the sin wave, but that's not important here. And I also explicitly defined the frequency in f(x), just to make it more clear what's going on -- as an implicitly frequency of 1 in confusing.)
主要要点是，由于自相关是通过一次沿一个轴移动一个点来计算的，因此自相关的轴就是时间轴.因此，我将其绘制为轴，然后可以读取其周期.在这里，我放大以清楚地看到它(下面是代码):
The main point is that since the autocorrelation is calculated by shifting along the axis one point at a time, the axis of the autocorrelation is just the time axis. So I plot that as the axis, and then can read the period off of that. Here I zoomed in to see it clearly (and the code is below):
# Testing the autocorrelation of numpy
import numpy as np
from matplotlib import pyplot as plt
num_samples = 1000
dt = 0.1
t = dt*np.arange(num_samples)
def autocor(x):
result = np.correlate(x,x,mode='full')
return result[result.size/2:]
def f(freq):
return np.sin(2*np.pi*freq*t)
plt.plot(t, autocor(f(.3)))
plt.xlabel("time (sec)")
plt.show()
也就是说，在上面，我将频率设置为0.3，该图显示的周期大约为3.3，这是预期的.
That is, in the above, I set the frequency to 0.3, and the graph shows the period as about 3.3, which is what's expected.
所有这些都表明，以我的经验，自相关通常适用于物理信号，但不适用于算法信号.例如，如果周期性信号跳过了某个步骤，这很容易抛出，这在算法中可能会发生，但在振动物体中发生的可能性较小.您可能会认为计算算法信号的周期应该是微不足道的，但是经过一番搜索就会发现它不是，甚至很难定义周期的含义.例如，系列:
All of this said, in my experience, the autocorrelation generally works well for physical signals but not so reliably for algorithmic signals. It's fairly easy to throw off, for example, if a periodic signal skips a step, which can happen with an algorithm, but is less likely with a vibrating object. You'd think that it should be trivial to calculate that period of an algorithmic signal, but a bit of searching around will show that it's not, and it's even difficult to define what's meant by period. For example, the series:
1 2 1 2 1 2 0 1 2 1 2 1 2
在自相关测试中不能很好地工作.
won't work well with the autocorrelation test.
这篇关于在算法信号中找到周期性的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
