1. 激活谷首页
  2. 激活谷笔记

做fft变换后0频非常大_为什么fft变换后0频率振幅大

• 激活谷笔记

做fft变换后0频非常大_为什么fft变换后0频率振幅大DFT, FFT与音频处理本文原来只是想把知乎当我的云笔记本记录一下FFT的，备着万一以后要搞下图像处理或者图形学里的基于FFT的海洋模拟。但是又想搞点差异化，于是脑子发热想写点音频处理，没想到坑也不浅。且因为我因为本科坑爹的课程安排导致我甚至还没上过信号处理相关的课，所以本文很可能有错误的地方，

DFT, FFT与音频处理
　　本文原来只是想把知乎当我的云笔记本记录一下FFT的，备着万一以后要搞下图像处理或者图形学里的基于FFT的海洋模拟。但是又想搞点差异化，于是脑子发热想写点音频处理，没想到坑也不浅。且因为我因为本科坑爹的课程安排导致我甚至还没上过信号处理相关的课，所以本文很可能有错误的地方，欢迎轻怼orz，图侵删

　　傅里叶变换

　　我们可能或多或少听说过，一个在长度为 P 的定义域里的可积(integrable)实数函数 f(x) 可以展开成无限个三角函数的和，也就是傅里叶级数 [1](Fourier series)：

　　f(x)= \frac{a_0}{2}+\sum^N_{n=1} \left( a_n\cos \left(\frac{2\pi nx}{P} \right) + b_n \sin \left(\frac{2\pi nx}{P} \right) \right)\\

　　其中序列 a_n, b_n 就是各个三角函数分量的权重。我们把这一大坨三角函数加权求和就可以逼近我们的原函数 f(x) ，就像泰勒级数那样。 另外呢，傅里叶级数有其他的形式，就是”振幅-相位”(amplitude-phase)形式：

　　 f(x) = \frac{A_0}{2} + \sum^N_{N=1}A_n \cos \left(\frac{2\pi nx}{P} - \phi_n \right)\\

　　其中振幅序列 A_n= \sqrt{a^2_n +b^2_n} ，相位 \phi_n=\text{atan2}(b_n,a_n)。这种形式可能是在做数字信号处理里面，我们比较关心的一种。 那么好了，现在我们有了一个已知的函数 f(x) ，要怎么求得这些三角函数的“权重”序列呢？在 N \rightarrow \infty 的时候，我们就用傅里叶变换就好啦：

　　 F(\xi)=\int^\infty_{-\infty}f(x)e^{-2\pi ix\xi} dx\\

　　不过还要提一下，前面的级数形式看起来还很正常，突然傅里叶变换就在 e 的指数上加了虚数看起来这么恐怖？？其实可以直接用欧拉公式展开，并把傅里叶变换的这个形式和前文提到的其他傅里叶级数的形式联系起来：

　　e^{ix}=cosx+isinx\\

　　这些三角函数的“权重”，就是频域里面的系数。我们每确定一次频率参数 \xi 并计算一次傅里叶变换的积分，就可以得到一个“权重”。做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图: 傅里叶变换

　　离散傅里叶变换(DFT)

　　DFT与频谱

　　下一个问题，在计算机的世界里，那必然不可能用无限个权重，输入的函数(信号)也在绝大部分时候不是连续的。于是我们有了离散傅里叶变换(DFT, Discrete Fourier Transform)，给定一个离散的实数序列 x(n) ，我们可以用DFT得到一个离散的频谱(spectrum) X(k)，其中频谱第 k 个点计算公式为：

　　 X(k)=\sum^{N-1}_{n=0} x(n) e^{-j \frac{2\pi}{N} k}\\

　　虽然理论上，你输入的是一堆复数(Complex)，输出也是一堆复数。但在做音频处理这种数字信号处理的场景下，输入是实数信号，你要对每个输入信号采样先加上一个0i，再拿去做DFT。 那这个DFT得到的频谱要怎么理解/解析/使用呢？音频处理的话，我们肯定是倾向于用振幅-相位形式的傅里叶级数来翻译啦。假设我们的频谱的复数值X(k)=a_k+b_kj，那么当前频谱点对应的谐波分量的振幅和相位为：

　　\begin{aligned} \\ A(k)&=\sqrt{a_k^2+b_k^2} \\ \phi(k)&= \text{atan2}(b_k,a_k) \\ \end{aligned} \\

　　有了上面的公式，我们可以把离散信号DFT的输出结果解析成振幅谱和相位谱，分别给出了每个(余弦)三角函数的振幅(权重)及相位。DFT及频谱的一些性质

　　1. 实数序列DFT频谱的对称性：因为一般情况下我们的输入序列 x(n) 是实数序列，则其频谱 X(k) 满足共轭对称性，即 X(k)=X^*(N-k)。这带来的一个结果是，振幅谱 A(k) 关于 k=N/2 圆周偶对称，相位谱 \phi(k) 关于 k=N/2 圆周奇对称。这意味着，我们可以做分析的时候会把频谱砍掉一半，一般情况下抛弃后半段。如果砍半之前的振幅谱叫双边振幅谱。做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图: DFT输出振幅谱的偶对称性

　　2. 折叠频率：假设输入信号的采样频率为 f_s，则这个 k=N/2 处对应的频率为折叠频率。其实折叠频率就是Nyquist Frequency，其值为 f_s/2。

　　3. 频谱点对应的谐波频率：振幅谱/相位谱每个点都有其对应的频率，N/2 个点均匀间隔地分布在 [0,f_s/2] 里面，相当于每个点都掌管一小段频率范围，每小段范围大小是 f_s/N 。每个频率小区域我们称为frequency bin。第k个frequency bin的频率 f_k=kf_s/N。

　　4. 输入信号补0：假设原始输入序列的点数为 N，在实际操作的时候我们可以通过在后面疯狂补0，把序列扩充到长度为 2N,3N,4N….，这样子DFT输出的频谱点数会增加，提高视觉上的分辨率，但遗憾的是频谱分辨能力并不会因此增加，补0前区分不开的频率分量，补0后还是区分不开。补0后只能让振幅谱和相位谱平滑连续一点，不需要自己做插值了。频谱泄漏

　　理想的连续Fourier Transform的频谱采样精度是无限，在计算机的世界里是没有这等好事的，DFT/FFT由于不能达到无穷的精度，有时每个频谱点对应的频率，不一定就是对应着输入信号的频率。例如我用1024 point FFT分析采样率20480Hz的信号，其源信号Nyquist Frequency是10240Hz，"有效频谱点"数为1024/2=512，则频谱每个点对应20Hz, 40 Hz, 60Hz, … 512 * 20Hz。如果我的输入信号的频率是20Hz的倍数还好，还能完美捕捉。但假设我有个输入信号不是20的倍数赫兹，那么这个信号不能被FFT完美的捕捉到，会有很多能量泄露并散布在整个频谱中，造成频谱泄漏(spectral leakage)。 下面是我写的FFT的一些test cases，输入信号为200Hz, 250Hz, 440Hz的正弦波叠加在一起，FFT window =1024 point，输出超过一定强度阈值的频谱点，应该可以体现Spectral Leakage：做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图1: 源信号采样率等于10000Hz做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图: 源信号采样率等于10240Hz

　　图1的frequency bin真是完美避开了源信号的分量频率，导致还挺严重的频谱泄漏。而图2的bin就真的可以全部精确覆盖到输入谐波的频率，就没有严重的泄漏了（但还是有点浮点数误差的）。这些会在主频率隔壁产生很多小山坡，降低频谱分辨率

　　快速傅里叶变换(FFT)

　　DIT Base2-FFT

　　其实按照公式，DFT的实现就两个for几行代码的事（当然你可能得先实现一个复数类 class Complex)：

　　那为什么后来搞出来个FFT呢？还不是因为计算速度和时间复杂度嘛，DFT的时间复杂度是 O(n^2)。一般我们记 W_N=e^{-j 2π/N} ，则DFT可以写成：

　　X(k)=∑_{n=0}^{N-1}x(n) W_N^{nk}\\

　　但是充分开发利用一下W_N的周期性和对称性，我们能节省掉很多计算量，最终推导出时间复杂度为 O(NlogN) 的快速傅里叶变换(FFT, Fast Fourier Transform)。具体的推导这里不展开，其他很多文章都还讲得挺详细的=.=。不过说回来，哪来那么多对称性和周期性呢？因为 W_N 是一个模长为1的复数，而复数可以看做一个复平面的点，也可以看做是二维笛卡尔坐标系里面的旋转变换。W_N 相当于是旋转了一小格 ( 2π/N 的弧度)，W_N^k 就是旋转了 k 小格，那么那些周期性和对称性就给了我们很多化简和做优化的余地。 假设输入的信号长度为 N=2^D,D \in Z^*，然后我们原输入序列 x(n) 分成奇数和偶数序号两组：

　　\begin{aligned} x_{even}(r) &= x(2r) \\ x_{odd}(r) &= x(2r + 1) \\ \end{aligned}\\

　　其中 r=-,1,...,\frac{N}{2}-1。 因为我们每次递归都可以把N点的计算分解成两个 N/2 点的计算，分解到最后一层，就直接从输入信号那取数据了。所以FFT靠化简和分治优化了DFT的时间复杂度。 但在实际操作中，出于各种因素的考虑，我们倾向于把这个递归形式的FFT完全展开。例如在GPU上用IFFT生成海洋表面的网格，在GPU上递归可就真的窒息了，所以flatten的FFT是很有优势的。展平的FFT可以使用蝶形算法来实现，算法流程也比较简单： 做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图：FFT蝶形变换流图

　　1. 用Reversed-bit index来交换元素顺序：假设我们的输入点数是 N=8=2^3，则指数 D=D=log_2⁡(N)=3 。我们把输入序列 x(n) 的序号用 D 比特整数表示，然后，我们把第 i 个元素和第 \text{reverseBit}(i)个元素调换位置。\text{reverseBit}(i) 就是把D bit整数的二进制位倒序排列，得到一个新的值。例如 6=110_2，要跟他对调位置的是 011_2=3。

　　2. D 次迭代的蝶形算法: 这个可以照着上图找规律，该乘的东西乘上，该加的东西也加上，然后抠进代码里。对于第 i 轮迭代，我们有 2^{D-i} 个蝶形单元，每个蝶形单元的“大小”为 2^i, i\in[1,D]。

　　上面的算法是按时间抽取(DIT)的Base2-FFT算法。还有另外一种Base2 FFT是按频率抽取(DIF)的Base2-FFT，这个算法打乱的是频域的元素，但这种方法的流图跟DIT-FFT不一样。

　　DFT, FFT的逆变换

　　离散傅里叶逆变换(IDFT, Inverse DFT) 是DFT的逆过程，IDFT的公式为：

　　x(n)=\frac{1}{N} \sum^{N-1}_{k=0}X(k)e^{i2\pi kn/N}=\frac{1}{N}\sum^{N-1}_{k=0}X(k)W^{-nk}_N\\

　　IDFT形式与正向DFT的基本形式很像，所以我们也不需要再次推导IFFT的蝶形变换了，直接把正向FFT的代码修改一下，加两个判断就可以变成IFFT。如果是IFFT算法，则把FFT中的W_N项的符号改一下，然后末尾多除一个N，就完事了。

　　音频处理

　　【**以下内容可能逐渐胡扯**】 音频是一维的离散数字信号。平时常用到的例如Adobe Audition等音频处理软件、FL Studio等编曲/混音软件里面就有巨量的音频处理算法。这些算法有基于时域的，也有基于频域的。我在平时用FL Studio做电子音乐的时候，会用到很多的混音特效插件，他们会在输入音频的基础上加点特效。这些专业的软件算法经过这么多年的发展，其底层实现肯定没有我想的那么trivial，但有一些特效算法肯定多多少少使用了Frequency Domain的方法，例如最基本的升降调(Pitch shifting)、均衡器(EQ, Equalizer)、滤波器(Filter)、Auto-tune、基本的修音防走调等。没了基于频域的方法，不少“歌手”音准修不了而被迫翻车，吴亦凡没了auto-tune也当不成加拿大电鳗（逃。

　　瞎搞的翻车现场

　　其实FFT离音频处理之间还是有点差距的，之前自己随便瞎搞了一下之后发现自己还是太拿衣服了orz。我先po几段自己瞎弄的pitch shifting的翻车音频，这些都是把音频生硬地分成N个frame做FFT，然后修改频谱之后IFFT得到输出。下场很惨，要不就是相位不连续疯狂爆音，要不就是FFT window大一点之后声音像鬼畜一样重复，要不就是FFT window太小导致频域精度不够，输出的和声非常不和谐………诸如此类，反正是没这么trivial的= =

　　【高能预警，小心音波冲击，可以提神】

　　STFT

　　所以如果我们要用频域算法来处理音频，正确的姿势是这样的：使用短时傅里叶变换(STFT, Short Time Fourier Transform)。定义短时傅里叶变换为[4]：

　　X(n,k)=\sum^\infty_{m=-\infty}x(m)h(n-m)e^{-j \frac{2\pi}{N} km}\\

　　其中 n 是时间，k 是频率， h(n) 是窗口函数。。这意思是我们要把原始的音频分成若干个FFT frame，每个frame都要乘上一个window function h(n)，例如可以使用hanning window[5]，用于“节选”一段离散信号序列：

　　h(n)=\frac{1}{2} \left(1-\cos\left( \frac{2\pi n}{N} \right)\right)\\做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图：hanning window

　　为了使用Phase Vocoder算法达成更好的效果，我们需要在时间轴上超采样STFT[2]，让FFT frame之间都有重叠。按照Phase Vocoder算法的原理，这个相邻frame之间的重叠至少要达到75%（即需要达到至少四倍的super-sampling，原因后面会提到）。

　　为什么我们要用STFT呢而不是使用囊括整段音频的FFT呢？那是因为音乐或者普通讲话的音频的频谱都是随时间变化的，如果我们用一个FFT frame把整段音频给包住了，那IFFT之后，音频就没有任何频谱变化，所有律动和变化都被平均掉了。好的，那么我们使用STFT也可以，但是STFT还不能直接处理音频，有些问题还待解决：

　　1. 时域和频域的采样精度不可兼得：一般Short Time指的是多短呢？假设我们的音频是44100Hz的采样率，那么做audio processing的FFT window一般是1024 / 2048 / 4096个采样这样子。假设我们使用1024的FFT window，那样子我们的频域精度就不够了，因为频谱里面的相邻Bin的center frequency都相互隔了有40Hz左右。但是如果用很大的FFT window，又会导致时域的采样率不够了，关于这个问题，可以从之前用16384 window的翻车音频多少感受到一点。

　　2. 频谱泄漏：这个前文已经提到过，即某个frequency bin的能量泄漏到其他bin里面，这会降低频谱的分辨力

　　Phase Vocoder

　　音频处理里面有一个比较通用的范式，Analysis-Process-Synthesis，前期的信号分析阶段就不只是做一下FFT。我们可以使用Phase Vocoder算法来改善采样精度和频谱泄露带来的弊端，于是我们会在Analysis阶段使用FFT频谱的相位信息去”纠正”每个frequency bin对应的频率。或者说，phase vocoder让我们可以用原始频谱的frequency bin的频率、相位、振幅三个信息，推导出“用相位信息纠正后的真实频率”及其对应振幅这两个信息[2]。做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图：好运的时候，每个frame都能恰好捕捉到信号的整数倍周期做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图：但更多的情况，分界点有一定phase offset

　　其大致原理是，每个FFT frame开头的初始相位信息是在疯狂暗示我们他的实际频率信息的，所以才有了用相位信息推断“实际频率”的操作。另外，之前提到要至少在时间轴上4倍超采样STFT，就是为了消除这种phase offset到frequency deviation推断的二义性，4倍超采样相当于把最多 2\pi 的phase offset分的足够细。具体的demo可以参考：Pitch Shifting Using The Fourier Transform

　　及配套代码 smbPitchShift()。 而Sythesis阶段则是Analysis阶段操作的逆过程。 相比起Analyze和Synthesis阶段，在Process阶段对Analyze过的频谱进行处理，则算是简单的了，例如我们要做升降调，我们只需要给每个频谱的频率乘一个缩放系数，拉伸频谱。例如要升一个八度，则新频谱对应的频率通通x2，降一个八度则是x0.5，升降一个半音(semi-tone)则是乘or除2^{1/12}。而均衡器的实现也比较直观，就是用一条EQ增益包络曲线去乘以原来的频谱就可以了（？吧）。做fft变换后0频非常大_为什么fft变换后0频率振幅大做fft变换后0频非常大_为什么fft变换后0频率振幅大图：FL Studio 20里面的EQ2均衡器插件

　　结论：音频处理算法的水好像也挺深，我就浅尝辄止了哈。（虽然本来也只是想记录一下FFT的一些细节和概念，没想到扯了这么多）希望下篇不扯淡，回归一下图形学。。。。

　　参考资料

　　[1] https://en.wikipedia.org/wiki/Fourier_series [2] Time And Pitch Scaling in Audio Procesing, http://www.surina.net/article/time-and-pitch-scaling.html [3] Pitch Shifting using the Fourier Transform, http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/ [4] Phase Vocoder, https://ccrma.stanford.edu/~jos/sasp/Phase_Vocoder.html [5] The Phase Vocoder: A tutorial, https://pdfs.semanticscholar.org/31d9/e1cc5d87c2b84cde2d4527b15b644544380e.pdf?_ga=2.8276405.686440233.1591512780-458074043.1591512780 [6] https://ww2.mathworks.cn/help/signal/ref/hann.html

2024最新激活全家桶教程，稳定运行到2099年，请移步至置顶文章：https://sigusoft.com/99576.html

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请联系我们举报，一经查实，本站将立刻删除。 文章由激活谷谷主-小谷整理,转载请注明出处：https://sigusoft.com/95655.html

赞 (0)
0
上一篇 2024年 5月 24日
下一篇 2024年 5月 24日

相关推荐

关注微信