在算法信号中寻找周期性

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 代码 sn-p 在 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

上面有 3 条重要信息：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

不适用于自相关检验。

Answer 2

更新。

@tom10 对自相关进行了全面调查，并解释了为什么自相关中的第一个凸点可以给出周期信号的周期。

我尝试了 FFT 和自相关这两种方法。他们的结果同意，尽管我更喜欢 FFT 而不是自相关，因为它更直接地为您提供了周期。

使用自相关时，我们只需确定第一个峰值的坐标。手动检查自相关图将显示您是否有“正确”峰值，因为您可以注意到周期（尽管对于 7 以上的素数，这变得不太清楚）。我相信你也可以制定一个简单的算法来计算“正确”的峰值。也许有人可以详细说明一些可以完成这项工作的简单算法？

例如，请参见下面的自相关序列图。 代码：

 1   # Plotting sequences satisfying, x_{i+1} = p-1 - (p*i-1 mod x_i)
 2   # with p prime and x_0 = 1, next to their autocorrelation.
 3   
 4   from __future__ import print_function
 5   import numpy as np
 6   from matplotlib import pyplot  as plt
 7   
 8   # The length of the sequences.
 9   seq_length = 10000
 10  
 11  upperbound_primes = 12 
 12  
 13  # Computing a list of prime numbers up to n
 14  def primes(n):
 15   sieve = [True] * n
 16   for i in xrange(3,int(n**0.5)+1,2):
 17     if sieve[i]:
 18         sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
 19   return [2] + [i for i in xrange(3,n,2) if sieve[i]]
 20  
 21  # The list of prime numbers up to upperbound_primes
 22  p = primes(upperbound_primes)
 23  
 24  # The amount of primes numbers
 25  no_primes = len(p)
 26  
 27  # Generate the sequence for the prime number p
 28  def sequence(p):
 29    x = np.empty(seq_length)
 30    x[0] = 1
 31    for i in range(1,seq_length):
 32      x[i] = p - 1 - (p * (i-1) - 1) % x[i-1]
 33    return x
 34  
 35  # List with the sequences.
 36  seq = [sequence(i) for i in p]  
 37  
 38  # Autocorrelation function.
 39  def autocor(x):
 40    result = np.correlate(x,x,mode='full')
 41    return result[result.size/2:]
 42  
 43  fig = plt.figure("The sequences next to their autocorrelation")
 44  plt.suptitle("The sequences next to their autocorrelation")
 45  
 46  # Proper spacing between subplots.
 47  fig.subplots_adjust(hspace=1.2)
 48  
 49  # Set up pyplot to use TeX.
 50  plt.rc('text',usetex=True)
 51  plt.rc('font',family='serif')
 52  
 53  # Maximize plot window by command.
 54  mng = plt.get_current_fig_manager()
 55  mng.resize(*mng.window.maxsize())
 56  
 57  k = 0 
 58  for s in  seq:
 59    k = k + 1
 60    fig.add_subplot(no_primes,2,2*(k-1)+1)
 61    plt.title("Sequence of the prime %d" % p[k-1])
 62    plt.plot(s)
 63    plt.xlabel(r"Index $i$")
 64    plt.ylabel(r"Sequence number $x_i$")
 65    plt.xlim(0,100)
 66    
 67    # Constrain the number of ticks on the y-axis, for clarity.
 68    plt.locator_params(axis='y',nbins=4)
 69  
 70    fig.add_subplot(no_primes,2,2*k)
 71    plt.title(r"Autocorrelation of the sequence $^{%d}x$" % p[k-1])
 72    plt.plot(autocor(s))
 73    plt.xlabel(r"Index $i$")
 74    plt.xticks
 75    plt.ylabel("Autocorrelation")
 76    
 77    # Proper scaling of the y-axis.
 78    ymin = autocor(s)[1]-int(autocor(s)[1]/10)
 79    ymax = autocor(s)[1]+int(autocor(s)[1]/10)
 80    plt.ylim(ymin,ymax)
 81    plt.xlim(0,500)
 82    
 83    plt.locator_params(axis='y',nbins=4)
 84  
 85    # Use scientific notation when 0< t < 1 or t > 10
 86    plt.ticklabel_format(style='sci',axis='y',scilimits=(0,1))
 87  
 88  plt.show()

使用 FFT 时，我们对序列进行傅里叶变换并寻找第一个峰值。第一个峰值的坐标给出了代表我们信号的最粗略的频率。这将给出我们的周期，因为最粗略的频率是我们的序列（理想情况下）振荡的频率。

请参阅下面的傅立叶变换旁边的序列图。

代码：

 1   # Plotting sequences satisfying, x_{i+1} = p-1 - (p*i-1 mod x_i)
 2   # with p prime and x_0 = 1, next to their Fourier transforms.
 3   
 4   from __future__ import print_function
 5   import numpy as np
 6   from matplotlib import pyplot  as plt
 7   
 8   # The length of the sequences.
 9   seq_length = 10000
 10  
 11  upperbound_primes = 12 
 12  
 13  # Computing a list of prime numbers up to n
 14  def primes(n):
 15   sieve = [True] * n
 16   for i in xrange(3,int(n**0.5)+1,2):
 17     if sieve[i]:
 18         sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
 19   return [2] + [i for i in xrange(3,n,2) if sieve[i]]
 20  
 21  # The list of prime numbers up to upperbound_primes
 22  p = primes(upperbound_primes)
 23  
 24  # The amount of primes numbers
 25  no_primes = len(p)
 26  
 27  # Generate the sequence for the prime number p
 28  def sequence(p):
 29    x = np.empty(seq_length)
 30    x[0] = 1
 31    for i in range(1,seq_length):
 32      x[i] = p - 1 - (p * (i-1) - 1) % x[i-1]
 33    return x
 34  
 35  # List with the sequences.
 36  seq = [sequence(i) for i in p]  
 37  
 38  fig = plt.figure("The sequences next to their FFT")
 39  plt.suptitle("The sequences next to their FFT")
 40  
 41  # Proper spacing between subplots.
 42  fig.subplots_adjust(hspace=1.2)
 43  
 44  # Set up pyplot to use TeX.
 45  plt.rc('text',usetex=True)
 46  plt.rc('font',family='serif')
 47  
 48  
 49  # Maximize plot window by command.
 50  mng = plt.get_current_fig_manager()
 51  mng.resize(*mng.window.maxsize())
 52  
 53  k = 0 
 54  for s in  seq:
 55    f = np.fft.rfft(s)
 56    f[0] = 0
 57    freq  = np.fft.rfftfreq(seq_length)
 58    k = k + 1
 59    fig.add_subplot(no_primes,2,2*(k-1)+1)
 60    plt.title("Sequence of the prime %d" % p[k-1])
 61    plt.plot(s)
 62    plt.xlabel(r"Index $i$")
 63    plt.ylabel(r"Sequence number $x_i$")
 64    plt.xlim(0,100)
 65    
 66    # Constrain the number of ticks on the y-axis, for clarity.
 67    plt.locator_params(nbins=4)
 68    
 69    fig.add_subplot(no_primes,2,2*k)
 70    plt.title(r"FFT of the sequence $^{%d}x$" % p[k-1])
 71    plt.plot(freq,abs(f))
 72    plt.xlabel("Frequency")
 73    plt.ylabel("Amplitude")
 74    plt.locator_params(nbins=4)
 75    
 76    # Use scientific notation when 0 < t < 0 or t > 10
 77    plt.ticklabel_format(style='sci',axis='y',scilimits=(0,1))
 78  
 79  plt.show()

要了解为什么 FFT 方法比自相关方法更方便，请注意我们有一个确定周期的明确算法：找到傅立叶变换的第一个峰值。对于足够数量的样本，这总是有效的。

见下表，由FFT方法得到，与自相关方法一致。

 prime   frequency   period
 2       0.00        1000.00
 3       0.50        2.00
 5       0.08        12.00
 7       0.02        59.88
 11      0.00        1000.00

以下代码实现了该算法，打印了一个表格，指定每个素数序列的频率和周期。

 1   # Print a table of periods, determined by the FFT method,
 2   # of sequences satisfying, 
 3   # x_{i+1} = p-1 - (p*i-1 mod x_i) with p prime and x_0 = 1.
 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 = 10000
 11  
 12  upperbound_primes = 12 
 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  # Function that finds the first peak.
 40  # Assumption: seq_length >> 10 so the Fourier transformed
 41  #        signal is sufficiently smooth. 
 42  def firstpeak(x):
 43    for i in range(10,len(x)-1):
 44      if x[i+1] < x[i]:
 45        return i
 46    return len(x)-1
 47  
 48  k = 0 
 49  for s in  seq:
 50    f = np.fft.rfft(s)
 51    freq  = np.fft.rfftfreq(seq_length)
 52    k = k + 1
 53    if k == 1:
 54      print("prime \t frequency \t period")
 55    print(p[k-1],'\t %.2f' % float(freq[firstpeak(abs(f))]), \
 56      '\t\t %.2f' % float(1/freq[firstpeak(abs(f))]))

我在上述所有代码中使用了 10000 个样本（seq_length）。随着我们增加样本数量，可以看到周期收敛到某个整数值（使用 FFT 方法）。

在我看来，FFT 方法是确定算法信号周期的理想工具，它仅受设备可以处理的样本数量的限制。