为了找到波形中的频率幅度,必须通过将波乘以它们正在搜索的频率来进行探测在两个不同的阶段(正弦和余弦)进行平均。相位是通过两者之间的关系找到的,其代码如下所示:
//simple pseudocode
var wave = [...]; //an array of floats representing amplitude of wave
var numSamples = wave.length;
var spectrum = [1,2,3,4,5,6...] //all frequencies being tested for.
function getMagnitudesOfSpectrum() {
var magnitudesOut = [];
var phasesOut = [];
for(freq in spectrum) {
var magnitudeSin = 0;
var magnitudeCos = 0;
for(sample in numSamples) {
magnitudeSin += amplitudeSinAt(sample, freq) * wave[sample];
magnitudeCos += amplitudeCosAt(sample, freq) * wave[sample];
}
magnitudesOut[freq] = (magnitudeSin + magnitudeCos)/numSamples;
phasesOut[freq] = //based off magnitudeSin and magnitudeCos
}
return magnitudesOut and phasesOut;
}
为了非常快地对很多频率执行此操作, FFT使用许多技巧。
使FFT比DFT快得多的技巧有哪些?
PS我曾尝试在网络上查看完整的FFT算法,但所有的窍门往往都被浓缩为一段漂亮的代码,而没有太多解释。在了解全部内容之前,我首先需要对这些有效更改的每个概念进行一些介绍。
谢谢。
#1 楼
$ N $点DFT的简单实现基本上是与$ N x N矩阵相乘。这导致$ \ mathcal {O}(N ^ 2)$的复杂性。最常见的快速傅立叶变换(FFT)算法之一是基数2 Cooley-Tukey抽取-时间FFT算法。这是基本的分而治之的方法。
首先将“旋转因子”定义为:
$$ W_N \ triangleq e ^ {-j \ frac {2 \ pi} {N}} $$
其中$ j \ triangleq \ sqrt {-1} $是虚数单位,则$ x [n] $的DFT $ X [k] $由
$$ X [k] = \ sum_ {n = 0} ^ {N-1} x [n] \,W_N ^ {kn} \,。$$
如果$ N $是偶数(并且$ \ tfrac {N} {2} $是整数) ,则可以将总和分为以下两个总和
$$ X [k] = \ sum_ {n = 0} ^ {N / 2-1} x [2n] W_N ^ {2kn} + \ sum_ {n = 0} ^ {N / 2-1} x [2n + 1] W_N ^ {k(2n + 1)} $$
其中第一个求和处理$ x [n]的偶数样本$,第二个具有$ x [n] $的奇数样本。定义$ x_e [n] \ triangleq x [2n] $和$ x_o [n] \ triangleq x [2n + 1] $并使用以下事实:
$ W_N ^ {k (2n + 1)} = W_N ^ {2kn} W_N ^ k $,并且
$ W_N ^ {2kn} = W_ {N / 2} ^ {kn} $,
此可以重写为
$$
\ begin {align}
X [k]&= \ sum_ {n = 0} ^ {N / 2-1} x_e [n] W_ {N / 2} ^ {kn} + W_N ^ k \ sum_ {n = 0} ^ {N / 2-1} x_o [n] W_ {N / 2} ^ {kn} \\
& = X_e [k] + W_N ^ k X_o [k]
\ end {align}
$$
其中$ X_e [k] $和$ X_o [k] $是$ \对$ x [n] $的偶数和奇数样本分别进行tfrac {N} {2} $点DFT变换。因此,我们只是将单个$ N $点DFT转换为两个较小的$ \ tfrac {N} {2} $点DFT。这会减少计算成本,因为
$$ 2 \ left(\ frac {N} {2} \ right)^ 2 + N
然后我们可以在这两个较小的DFT上重复相同的过程。这种分而治之的方法允许达到$ \ mathcal {O}(N \ log N)$的复杂度,这比我们通过朴素DFT获得的$ \ mathcal {O}(N ^ 2)$更好。实现(如leftaroundabout的答案所充分说明)。
评论
$ \ begingroup $
您愿意列出每个变量代表什么吗?我对此还比较陌生,因此W,j,X(),N和k尚未为我定义。
$ \ endgroup $
– Seph Reed
17年6月8日在20:44
$ \ begingroup $
$ W $已经在我的答案中定义了。我试图更好地定义其他一些符号。 $ k $表示频域中的索引,$ n $表示时域。
$ \ endgroup $
– anpar
17年6月8日在20:53
#2 楼
http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629DFT,尺寸16
复杂性的区别就很明显了,不是吗?
这就是我对FFT的理解。
/>
首先,我总是将傅立叶变换主要考虑为连续函数的变换,即双射映射$ \ operatorname {FT}:\ mathcal {L} ^ 2(\ mathbb {R})\ to \ mathcal {L} ^ 2(\ mathbb {R})$。从这个角度来看,很明显,实际上并没有必要进入“最深层次”并遍历单个元素,因为“单个元素”是实线上的单点,其中有无数个无限大的点。 。
那么,这种转换又如何定义得很好呢?好吧,至关重要的是,它不能在通用函数空间$ \ mathbb {R} \ to \ mathbb {C} $上运行,而只能在(Lebesgue-,square-)可积函数的空间上运行。现在,这种可集成性并不是一个很强的属性(比可微性等要弱得多),但是它确实要求该功能“在本地可用大量信息来区分”。这种描述是通过短时傅立叶变换的系数给出的。†最简单的情况是您的函数是连续的,并且将其划分为很小的区域,因此每个区域基本上都是恒定的。然后,每个STFT具有最强的第零项。如果忽略(无论如何衰减)其他系数,则每个域只是一个数据点。在所有这些短时LF极限系数中,您可以进行离散傅立叶变换。实际上,这就是对实测数据进行任何FT时所做的事情!
但是,测得的数据不一定与基本物理量相对应。例如,当您测量某些光强度时,实际上只是在测量电磁波的幅度,该电磁波的频率本身太高而无法使用ADC进行采样。但是显然,尽管光波的频率是疯狂的,但您也可以便宜地计算出采样光强度信号的DFT。
这可以理解为FFT便宜的最重要原因:
不必试图从最高水平看到各个振荡周期。取而代之的是,仅转换已经在本地进行了预处理的高级信息。
这还不是全部。 FFT的优点是仍然可以为您提供完整DFT所能提供的所有信息。即采样准确的光束电磁波时,您还将获得所有信息。可以通过转换光电二极管信号来实现吗? –您能从中测量出确切的光频率吗?
好吧,答案是不,您不能。也就是说,除非您应用其他技巧。首先,您至少需要在短时间内大致测量频率。好吧,用光谱仪是可能的。但这只能达到$ \ Delta \ nu = 1 / {\ Delta t} $的精度,这是一个典型的不确定性关系‡。
通过总体上较长的时间跨度,我们也应该能够缩小频率不确定度。如果您不仅可以在本地测量粗略的频率,而且还可以测量波的相位,则确实可以实现。您知道,如果一秒钟后再看一眼,则1000 Hz信号将具有完全相同的相位。而1000.5 Hz信号虽然在短范围内无法区分,但一秒钟后将具有反相。
幸运的是,该相位信息可以很好地存储在单个复数中。这就是FFT的工作方式!它以许多小的局部转换开始。这些价格便宜-一方面,显然是因为它们仅使用少量数据,但是其次,因为他们知道,由于时间跨度短,它们无论如何都无法精确地解析频率-因此,即使您做很多这样的转变。
但是,它们确实也记录了相位,因此您可以在顶层使频率分辨率更加精确。所需的转换还是很便宜的,因为它本身不会打扰任何高频振荡,而只会打扰经过预处理的低频数据。在这一点上有点循环。让我们称其为递归就可以了...
‡这种关系不是量子力学,但是海森堡不确定性实际上具有相同的根本原因。
评论
$ \ begingroup $
关于这个问题的精美图画。 :-)
$ \ endgroup $
–罗伯特·布里斯托-约翰逊
17年9月9日在18:12
$ \ begingroup $
您是否不喜欢在任何地方都重复出现且从未在任何地方进行过实际解释的图:)
$ \ endgroup $
–user541686
17年6月9日在23:03
$ \ begingroup $
我刚刚看了anpar的答案,就明白了这张照片。
$ \ endgroup $
–JDługosz
17年6月11日在16:48
#3 楼
这是添加到Robert好的答案中的图片,展示了操作的“重用”,在本例中为8点DFT。图中的“旋转因子”使用符号$ W_N ^ {nk} $表示,该符号等于$ e ^ {j2 \ pi \ frac {nk} {N}} $注意所示的路径和下面的方程式显示了频率仓X(1)的结果,如罗伯特方程式所示。
虚线与实线没有什么不同,只是为了明确求和联接的位置。
#4 楼
本质上,在直接根据求和计算朴素DFT时:$$ X [k] = \ sum \ limits_ {n = 0} ^ {N-1} x [n] \,e ^ {j 2 \ pi \ frac {nk} {N}} $$
存在用于旋转因子$ e $ N $的表查找^ {j 2 \ pi \ frac {nk} {N }},$ N $复数乘法和$ N-1 $加法。那只是为了$ X [k] $的一个值和$ k $的一个实例。那么朴素的DFT会丢弃所有中间数据,并再次对所有数据进行$ X [k + 1] $。
,因此FFT保留了一些中间数据。
FFT还将稍微利用旋转因子,以便相同的因子可用于数据的中间组合。
#5 楼
我是一个有视觉感的人。我更喜欢将FFT想象成矩阵技巧,而不是求和技巧。从高层次进行解释:
朴素的DFT独立地计算每个输出样本,并在每次计算中使用每个输入样本(经典N²算法)。
常见的FFT使用DFT定义中的对称性和模式在“层”(对数N层)中进行计算,每个层的每个样本具有恒定时间要求,因此创建了N log N算法。 br />更多细节:
可视化这些对称性的一种方法是将DFT视为1×N矩阵输入乘以所有复杂指数的NxN矩阵。让我们从“基数2”的情况开始。我们将拆分矩阵的偶数和奇数行(对应于偶数和奇数输入样本),并将其视为两个单独的矩阵乘法,它们相加在一起可获得相同的最终结果。
现在看一下这些矩阵:在第一个矩阵中,左半部分与右半部分相同。另一方面,右半部分是左半部分x -1。这意味着我们只需要真正使用这些矩阵的左半部分进行乘法,并通过乘以1或-1廉价地创建右半部分。接下来,观察第二个矩阵与第一个矩阵的不同之处,因为每列中的因子相同,因此我们可以将其分解为因子并将其乘以输入,因此现在偶数和奇数样本都使用相同的矩阵,但是需要乘数第一。最后一步是观察得到的N / 2×N / 2矩阵与N / 2 DFT矩阵相同,我们可以一次又一次地进行操作,直到达到1×1矩阵,其中DFT是恒等式。
要泛化除基数2之外的其他内容,您可以查看每三行的拆分和三列的列,或者每四列的第一个,等等。 -pad,FFT和截断,但这超出了此答案的范围。
请参见:http://whoiskylefinn.com/MatrixFFT.html
评论
$ \ begingroup $
主FFT,各种FFT。使用零垫不是唯一的选择。抱歉,我发现零填充已被过度使用。一个小问题,我不明白您所说的“每个样本具有恒定时间要求的每一层”是什么意思,如果您能解释一下,那就太好了。
$ \ endgroup $
–邪恶
17年6月9日在22:30
$ \ begingroup $
对不起,我并不是说零填充就是这种方式,只是想指出更多的阅读内容。 “层”表示递归或从N DFT到2 N / 2 DFT的转换,每个样本的时间恒定,表示此步骤为O(N)。
$ \ endgroup $
– kylefinn
17年6月9日在23:20
$ \ begingroup $
到目前为止,在所有描述中,这似乎是使复杂的问题变得简单的最接近的描述。但是,它缺少的最大的东西就是这些矩阵的一个例子。你碰巧有一个吗?
$ \ endgroup $
– Seph Reed
17年6月10日在18:30
$ \ begingroup $
上传了此内容,应该会有所帮助:whoiskylefinn.com/MatrixFFT.html
$ \ endgroup $
– kylefinn
17年6月13日在1:51
#6 楼
DFT进行了蛮力的N ^ 2矩阵乘法。FFT巧妙地利用了矩阵的属性(对矩阵乘法进行去归一化)以降低计算成本。
让我们首先来看一个小的DFT:
W = fft(eye(4));
x = rand(4,1)+ 1j * rand(4 ,1);
X_ref = fft(x);
X = W * x;
太好了,我们可以通过填充FFT函数的矩阵,通过小的4x4(复杂)矩阵乘法来代替对FFTW库的MATLAB调用。那么这个矩阵是什么样的呢?
N = 4,
Wn = exp(-1j * 2 * pi / N),
f =((0:N-1)'*(0:N-1))
f =
0 0 0 0
0 1 2 3
0 2 4 6
0 3 6 9
W = Wn 。^ f
W =
1 1 1 1
1 -i -1 i
1 -1 1 -1
1 i -1 -i
每个元素都是+ 1,-1,+ 1j或-1j。显然,这意味着我们可以避免完全复杂的乘法。此外,第一列是相同的,这意味着我们将x的第一个元素一遍又一遍地乘以相同的因子。
结果证明,Kronecker张量积,“旋转因子”和置换矩阵根据翻转的二进制表示形式更改索引的方法既紧凑,又提供了另一种观点,即如何将FFT作为一组稀疏矩阵运算进行计算。
下面的几行内容是一个简单的频率抽取( DIF)基数2前向FFT。虽然这些步骤似乎很麻烦,但可以方便地重用正向/反向FFT,radix4 / split-radix或及时抽取,同时可以很好地表示在现实世界中倾向于实现就地FFT,我相信。
N = 4;
x = randn(N,1)+ 1j * randn(N,1);
T1 = exp(-1j * 2 * pi *([zeros(1,N / 2),0:(N / 2-1)])。'/ N),
M0 = kron(eye(2),fft(eye(2))),
M1 = kron(fft(eye(2)),eye(2)),
X = bitrevorder(x。'* M1 * diag(T1)* M0),
X_ref = fft(x)
CF Van Loan在这方面有一本很棒的书。
#7 楼
如果您想从智慧之火中喝酒,我建议:“快速变换-算法,分析,应用程序”
Douglas F. Elliott,K. Ramamohan Rao
它涵盖了FFT,Hartley,Winograd和应用程序。
一个优点是,它说明了FFT是如何以位反转顺序进行一组稀疏矩阵分解的。
评论
“ DFT”不是指算法:它是指数学运算。 “ FFT”是指用于计算该操作的一类方法。只是想指出,在您的代码示例中使用sudo可能会造成混淆,因为这是计算机世界中众所周知的命令。您可能是指psuedocode。
@nwfeather他可能的意思是“伪代码”。