查找算法信号中的周期性

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.这显示了数字的频率.我不知道这如何能提供有关数字序列周期性的任何有用信息.

最后一行显示了我使用自相关的尝试,受到了如何使用numpy.correlate进行自相关的公认答案的启发？.

它给出了以下图表

显然,我们看到所有质数的数字序列递减.

使用以下简单的python-code片段在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 vs 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

不适用于自相关测试。