python - 在算法信号中寻找周期性

Question

在检验关于以下递归关系的猜想

,

它声称数字序列具有某种周期性，我编写了一个 python 程序来计算序列并将它们打印在表格中。

 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  

现在我想研究我计算的这些序列的周期性。在网上四处寻找后，我发现自己有两个选择:

• 对数据进行自相关并寻找第一个峰值。这应该给出周期的近似值。
• 对数据执行 FFT。这显示了数字的频率。我看不出这如何能提供有关数字序列周期性的任何有用信息。

最后几行显示了我使用自相关的尝试，灵感来自 How can I use numpy.correlate to do autocorrelation? 的已接受答案.

它给出了下面的情节

很明显，我们看到所有素数的递减序列。

当使用以下简化的 python 代码片段在 sin 函数上测试相同的方法时

 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()

我得到了类似的结果，它给出了以下正弦函数图

例如，我怎样才能读出正弦函数情况下的周期性？

无论如何，我不理解导致峰值的自相关机制，这些峰值提供了信号周期性的信息。有人可以详细说明吗？在这种情况下，您如何正确使用自相关？

另外，我在实现自相关时做错了什么？

欢迎就确定数字序列中的周期性的替代方法提出建议。

Answer 1

这里有很多问题，所以我将开始描述自相关如何从“3”的情况产生周期，即第一张图片的第二个子图。

对于素数 3，序列是(在不太一致的开始之后)1,2,1,2,1,2,1,2,...。为了计算自相关，数组基本上是相对于自身平移的，所有对齐的元素相乘，所有这些结果相加。所以它看起来像这样，对于一些测试用例，其中 A 是自相关:

 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

从上面可以看出三个关键信息:1.自相关 A 在相似元素排列和相乘时具有更大的值，这里在每个其他步骤。 2. 自相关的指数对应于相对偏移。 3.当对整个数组进行自相关时，如此处所示，总是会出现一个向下的斜坡，因为在每个连续的移位中加在一起产生值的点数都会减少。

所以在这里您可以看到为什么“质数 3”在您的图表中有一个周期性的 20% 的凸起:因为当它们对齐时求和的项是 1+4，而当它们不对齐时是 2+2，即，5 对 4。这就是您在阅读期间时要寻找的颠簸。也就是说，这里显示周期为 2，因为这是您的第一个颠簸的索引。 (另外，顺便说一句，在上面我只计算了对的数量，看看这个已知的周期性如何导致你在自相关中看到的结果，也就是说，人们通常不想考虑对的数量.)

在这些计算中，如果您在进行自相关之前先减去均值，则相对于底数的凸起值将会增加。如果您使用末端经过修剪的数组进行计算，则可以删除斜坡，因此始终存在相同的重叠；这通常是有道理的，因为通常人们正在寻找比整个样本波长短得多的周期性(因为需要多次振荡才能定义良好的振荡周期)。

对于正弦波的自相关，基本答案是周期显示为第一个凸点。除了应用时间轴外，我重新绘制了情节。在这些事情中使用实时轴总是最清楚的，所以我稍微更改了您的代码以包含它。 (此外，我用一个适当的矢量化 numpy 表达式替换了列表理解来计算正弦波，但这在这里并不重要。我还明确定义了 f(x) 中的频率，只是为了更清楚地说明发生了什么——作为 1 的隐含频率令人困惑。)

要点在于，由于自相关是通过沿轴一次移动一个点来计算的，因此自相关的轴就是时间轴。所以我将其绘制为轴，然后可以读取其中的周期。这里我放大看清楚了(下面是代码):

# 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，这是预期的。

综上所述，根据我的经验，自相关通常适用于物理信号，但适用于算法信号则不太可靠。例如，如果周期性信号跳过一个步骤，这很容易被忽略，算法可能会发生这种情况，但振动物体的可能性较小。您可能认为计算算法信号的周期应该是微不足道的，但稍微搜索一下就会发现事实并非如此，甚至很难定义周期的含义。例如系列:

1 2 1 2 1 2 0 1 2 1 2 1 2

不适用于自相关检验。