如何理解周期发现下的离散傅里叶变换的结果？答案

【问题标题】：How to understand the result of Discrete Fourier Transform under period finding?如何理解周期发现下的离散傅里叶变换的结果？
【发布时间】：2021-08-08 07:49:45
【问题描述】：

我正在学习如何使用离散傅里叶变换（DFT）找到关于a^x mod(N) 的周期，其中x 是正整数，a 是任意素数，N 是乘积两个质因数p 和q。

比如2^x mod(15)的周期是4，

>>> for x in range(8):
...     print(2**x % 15)
...

Output: 1 2 4 8 1 2 4 8
                ^-- the next period

DFT的结果如下，

（引自O'Reilly Programming Quantum Computers chapter 12）

4个单位间距的尖峰有4个，我认为后4个表示周期为4。

但是，当N 是 35 并且句点是 12 时

>>> for x in range(16):
...     print(2**x % 35)
...

Output: 1 2 4 8 16 32 29 23 11 22 9 18 1 2 4 8
                                       ^-- the next period

在这种情况下，有 8 个大于 100 的尖峰，其位置分别为 0、5、6、11、32、53、58、59。

位置序列是否暗示了幻数 12？以及如何从右图中理解“12 个均匀分布的尖峰”？

（引自O'Reilly Programming Quantum Computers chapter 12）

【问题讨论】：

标签： math fft dft

【解决方案1】：

查看How to compute Discrete Fourier Transform? 和所有子链接，尤其是How do I obtain the frequencies of each value in an FFT?。

如您所见，DFT 结果的i-th 元素（从0 到n-1 包括在内）表示Niquist 频率

f(i) = i * fsampling / n

而 DFT 结果仅使用那些正弦频率。因此，如果您的信号确实有不同的信号（即使频率或形状略有不同），就会出现混叠。

混叠正弦曲线在 DFT 输出中产生 2 个频率，一个较高的频率，一个较低的频率。

任何锋利的边缘都会被转换成许多频率（通常是像你上一个例子一样的连续频谱）

f(0) 没有频率，代表直流偏移。

除此之外，如果 DFT 的输入是实域，那么 DFT 结果是对称的，这意味着您只能使用 结果的前半部分作为第二部分只是镜像（不包括f(0)），这是有道理的，因为您不能在实际域数据中表示大于fsampling/2 的频率。

结论：

您无法获得 DFT 使用的信号频率，因为有无数种方法可以计算此类信号。 DFT 正在使用正弦波重构信号，而您的信号肯定不是正弦波，因此结果与您的想法不符。

通过为 DFT 正确选择 n 来将 niquist 频率与您的频率匹配，但是在不知道前面的频率的情况下您无法做到这一点...

有可能从它的 2 个别名中计算出奇异的正弦波频率，但是您的信号不是正弦波，因此无论如何这不适用于您的情况。

我会使用不同的方法来确定整数数字信号的频率：

计算信号直方图

数一数每个数有多少
测试可能的频率

您可以暴力破解所有可能的信号周期并测试后续周期是否相同，但是对于大数据来说这不是最佳选择...

我们可以使用直方图来加快速度。因此，如果您从直方图中查看频率 f 和周期 T 在大小为 n 的数据中的周期信号的计数 cnt(ix)，那么信号的周期应该是所有计数的公共分频器
```
T = n/f
k*f = GCD(all non zero cnt[i])
```
k 将 GCD 结果相除。但是，如果n 不是T 的精确倍数，或者信号中有噪声或轻微偏差，这将不起作用。然而，我们至少可以估计 GCD 并测试所有频率，这些频率仍然比蛮力更快。

所以对于每个计数（不考虑噪音）它应该遵守这个：
```
cnt(ix) = ~ n/(f*k)
k = { 1,2,3,4,...,n/f} 
```
所以：
```
f = ~ n/(cnt(ix)*k)
```
所以如果你收到这样的信号：
```
1,1,1,2,2,2,2,3,3,1,1,1,2,2,2,2,3,3,1
```
那么直方图将是cnt[]={0,7,8,4,0,0,0,0,...} 和n=19，因此每个使用的元素在每个n 的周期中计算f 会导致：
```
f(ix) = n/(cnt(ix)*k) 
f(1)  = 19/(7*k) = ~ 2.714/k
f(2)  = 19/(8*k) = ~ 2.375/k
f(3)  = 19/(4*k) = ~ 4.750/k
```
现在实际频率应该是结果的公共分频器 (CD)，因此将最大和最小计数向上和向下舍入（忽略噪声）导致以下选项：
```
f = CD(2,4) = 2
f = CD(3,4) = none
f = CD(2,5) = none
f = CD(3,5) = none
```
所以现在测试频率（幸运的是，在这种情况下它只有一个有效）每 19 个样本有 2 个周期，这意味着 T = ~ 9.5 所以测试向上和向下舍入...
```
signal(t+ 0)=1,1,1,2,2,2,2,3,3,1,1,1,2,2,2,2,3,3,1
signal(t+ 9)=1,1,1,2,2,2,2,3,3,1 // check 9 elements
signal(t+10)=1,1,2,2,2,2,3,3,1,?   // check 10 elements
```
如您所见，signal(t...t+9)==signal(t+9...t+9+9) 表示句点为T=9。

【讨论】：