当前位置： 首页 > news >正文

重拾傅里叶变换

news 2026/4/20 17:53:05

1. 核心出装

起因是晨郢给了这么一个思考题：

\[\binom {21}0-\binom{21}2+\binom{21}4-\binom{21}6+\cdots +\binom{21}{20}=? \]

水神（上次数学 150 那位）一下课就冲过来向我宣称，答案就是

\[\Re\left((1+i)^{21}\right)=-1024 \]

我大受震撼，上机验证发现是对的，就近查询，机房大蛇告诉我这玩意和单位根反演有密切联系，还说这个系数可以 FFT 出来，题设序列可以拓展到任意长度、周期、正负号。然后傅里叶变换这玩意迟早要弄透，于是有了这篇。

多项式拆解

\(\Re\left((1+i)^{21}\right)\) 是凭数感观察出来的，这个取实部看着不是很能拓展。写成共轭的标准表达：

\[\frac 12(1+i)^{21}+\frac 12(1-i)^{21} \]

我们想要知道，\(1\) 后面跟着的 \(\pm i\) 是怎么出来的。重写问题：

\[\sum_{k=0}^m a_k\binom mk \]

我们希望将其表达成：

\[\sum_{j=0}^{n-1} b_j(1+x_j)^m\\ =\sum_{k=0}^m \sum_{j=0}^{n-1} b_j x_j^k \binom mk \]

也即要求：

\[\boxed{a_k=\sum_{j=0}^{n-1} b_jx_j^k}\quad \forall k\in[0,m]\quad \bold ① \]

单位根构造

接下来的构造是完备且一体的，并没有特定的 insight 指引我们观察出这个构造，只能去欣赏它的简洁优美。

我们取

\[x_j=\omega_n^j=\exp\left(\frac {2ij\pi} n\right) \]

其中 \(\omega\) 为复数域上的单位根：

\[\begin{cases}(\omega_n^1)^n=1\\\omega_n^j=(\omega_n^1)^j\end{cases} \]

到现在我们依然没有解释 \(n\) 是怎么得出的。实则 \(n\) 为题设序列 \(a\) 的任意周期（e.g. 这个例题的 \(n=4,8,12,16,20,22\)）。又因为 \(\color{lime}x_j^{k+n}=x_j^k\)，所以 \(x_j^k\) 关于 \(k\) 也有 \(n\) 的周期。那么 ① 的条件可以直接改成 \(\forall k\in[0,n)\)。

至此我们扔掉了组合数并把问题的范围缩小到了一个 \(n\)。

DFT

要把问题转成 DFT 的形式，需要利用单位根的性质

\[\begin{aligned} a_k&=\sum_{j=0}^{n-1} b_j\color{red}{(\omega_n^j)^k}\\ &=\sum_{j=0}^{n-1} b_j\color{red}{(\omega_n^k)^j} \end{aligned} \]

孩子们，\(a\) 和 \(b\) 互为同一个多项式的点值和系数表示法！已知 \(a\) 求 \(b\) 就是一个 IDFT 的事情。（注意 \(b_j\in \C\)）

#include <cstdio>
#include <cmath>const double pi=acos(-1),eps=1e-7;
const int N=(1<<21)|1;struct cplx{double r,i;}A[N];
cplx operator+(const cplx &x,const cplx &y){return (cplx){x.r+y.r,x.i+y.i};
}
cplx operator-(const cplx &x,const cplx &y){return (cplx){x.r-y.r,x.i-y.i};
}
cplx operator*(const cplx &x,const cplx &y){return (cplx){x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r};
}int r[N],n;void FFT(cplx* a,int typ){for(int i=0;i<n;++i) if(i<r[i]){cplx tmp=a[i];a[i]=a[r[i]];a[r[i]]=tmp;}for(int w=1;w<n;w<<=1){cplx o1=(cplx){cos(pi/w),typ*sin(pi/w)};for(int i=0;i<n;i+=(w<<1)){cplx ok=(cplx){1,0};for(int j=i;j<i+w;ok=ok*o1,++j){cplx u=a[j],d=ok*a[j+w];a[j]=u+d;a[j+w]=u-d;}}}
}void wr(double r,double i){printf("(%.4lf",r);if(i>-eps) putchar('+');printf("%.4lfi)",i);}int main()
{scanf("%d",&n);int C;scanf("%d",&C);for(int i=0;i<n;++i) scanf("%lf",&A[i].r);int m=n,l=0;for(n=1;n<m;n<<=1) ++l;for(int i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));FFT(A,-1);double w=2*pi/n;for(int i=0;i<n;++i){wr(A[i].r/n,A[i].i/n);wr(1+cos(w*i),sin(w*i));printf("^%d",C);if(i<n-1) puts(" +");}
}

输入 #1：

4 21
1 0 -1 0

输出 #1：

(0.0000+0.0000i)(2.0000+0.0000i)^21 +
(0.5000+0.0000i)(1.0000+1.0000i)^21 +
(0.0000+0.0000i)(0.0000+0.0000i)^21 +
(0.5000+0.0000i)(1.0000-1.0000i)^21

输入 #2：

2 20
1 -1

输出 #2：

(0.0000+0.0000i)(2.0000+0.0000i)^20 +
(1.0000+0.0000i)(0.0000+0.0000i)^20

输入 #3：

8 114514
1 3 -1 2 -4 2 -1 -2

输出 #3：

(0.0000+0.0000i)(2.0000+0.0000i)^114514 +
(0.3598-0.4419i)(1.7071+0.7071i)^114514 +
(-0.1250-0.6250i)(1.0000+1.0000i)^114514 +
(0.8902-0.4419i)(0.2929+0.7071i)^114514 +
(-1.2500+0.0000i)(0.0000+0.0000i)^114514 +
(0.8902+0.4419i)(0.2929-0.7071i)^114514 +
(-0.1250+0.6250i)(1.0000-1.0000i)^114514 +
(0.3598+0.4419i)(1.7071-0.7071i)^114514

（对 \(n\ne 2^k,k\in \Z\) 是没法补零做 FFT 的，会改变周期，只能朴素 \(n^2\) 去做）

到这里原问题就结束了。考察这个转换对计算原式的意义，实则是聊胜于无，因为你涉及 \((\cos\theta+1+i\sin \theta)^m\) 状物的计算，手算只在 \(n=4,\theta=\frac {k\pi}2\) 时能算，计算机去算的话无非就是整型炸精变成了浮点炸精（而且如果是整数域上取模算的话我为啥不直接算），所以这题相当无聊。

！？DFT？！

但是我突然发现这个工具非常强力，可以解决一个我从初二想到现在都没能实现的事情。

\[\begin{aligned} a_k&=\sum_{j=0}^{n-1} b_j\exp\left(\frac {2ijk\pi}n\right)\\ &=\sum_{j=0}^{n-1} \Re(b_j)\cos \frac{2jk\pi}n - \Im(b_j)\sin \frac {2jk\pi}n\\ &=\sum_{j=0}^{n-1} |b_j|\sin\left(\frac{2j\pi}n\cdot k+\varphi_j\right) \end{aligned} \]

这是我第一次理解 DFT 和频率分解间的直接联系。这个式子完全等效于，你把 \(n\) 个采样点的离散信号 \(a\) 分解成波长为 \(\frac nj\) 的多个正弦波的叠加。那么单旋律音高提取这件事就完全可以做了。

2. 音高提取

音频读取

我们有开源精神，可以把每个采样点的功率给读出来，下面记 \(A\)。.wav 文件可以理解成是最基本、不带压缩的时间-功率映射关系。

频率分解

上面那个式子是 IDFT 出来的，实则 DFT 和 IDFT 只差前面一个 \(\frac 1n\) 的系数、并且各个 \(b_j\) 互为共轭，又由前置知识，所以你只需要对 \(A\) 做 FFT 然后取各个频率上 \(|b_j|\) 的 \(\arg\max\) 即为（近似）对应音高。

这玩意乍听上去非常色情，但实则你要考虑下面两个棘手的问题：

时间窗

直接频率分解的前提是，你要确定本段 \(A\) 序列所对应的时间片段里，音高是保持不变的。因此你需要对时间做微元（即时间窗 window。注意与采样点的概念区分，采样点是最小不可再分的时间单位）。令采样率 \(W(W=48000,44100,\cdots)\)，\(U\) 为 1s 内分出 window 的个数，则 \(n=\frac WU\)。以下均取 \(W=48000\) 的最理想情况。

  • \(n\) 不能太小，否则分析出来频率错位严重。你要处理的频率最低值约为 \(f_m=200\) Hz（G0 左右），波长最大值 \(\lambda_m=\frac W{f_m}=240\)，那么你最少应该要求 \(n>5\lambda_m\)。
  • \(U\) 不能太小，你 \(150\) bpm 的十六分音即可达到 1s 内 \(U_0=10\) 个音符，显然 \(U\ge 2U_0\) 才能对这样的情形做有效分析。

实际测出来的结论是，无论 \(W\) 为多少我们都取 \(U=25\)，这是比较优的。

误差分析

DFT 出来的频率是离散的。它得到的分量永远满足：一个 window 刚好为 \(j\) 个波长，\(j\in \Z\)。那么你得到的频率也只有 \(\frac Wj,j\in \N^*\)（其实由奈奎斯特采样定理 \(j=1\) 也取不到，不过无所谓）。以下 \(W=48000\)。

首先你就要面临取整带来的直接误差。考虑 \(f\) 的计算式：

\[\lambda=\frac W{f_{real}}\\ j=\lfloor\frac n{\lambda}+\frac 12\rfloor\\ f=Uj \]

那么相对误差 \(\frac f{f_{real}}\) 会去到 \(\frac {n+\frac 12 \lambda}n\)，极端情况为 \(n=2048\)（要 FFT）、\(\lambda=\lambda_m=240\)，误差约为 \(1.05\)。

如果你认为这样的误差无所谓，下面还有一个更大的问题。你把 \(\frac Wj,j\in \N^*\) 列出来然后拿这玩意跟 \(\text{A1}=440\) Hz 的十二平均律标准音高列对比：

... 1200.00 1230.77 1263.16 1297.30 1333.33 1371.43 1411.76 1454.55 1500.00 1548.39 1600.00 1655.17 1714.29 1777.78 ...
... 1174.66 1244.51 1318.51 1396.91 1479.98 1567.98 1661.22 1760.00 ...

你会发现有些频率点刚好卡在两半音中间，在峰值落到这些点时你无法决策具体取两边哪个半音，这很难受，你需要处理一手邻域。

归根到底，是频谱的离散性导致了这些问题。有没有把时域上离散的信号转换成在频谱上连续表达的方法呢？有的兄弟有的。

3. DTFT

设无限序列 \(x[n]\) 的 DTFT 为

\[X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-i\omega n} \]

DFT 实际上是：

  1. 截断信号
  2. 周期延拓
  3. 对 DTFT 进行等间隔采样

采样点为

\[\omega_k=\frac{2\pi k}{N} \]

于是

\[X[k]=X(\omega_k) \]

变换 时域 频域
傅里叶级数 周期连续信号 离散频谱
傅里叶变换 非周期连续信号 连续频谱
DTFT 非周期离散信号 连续周期频谱
DFT 周期离散信号 离散周期频谱

其中：

  • DTFT：频域周期 \(2\pi\)
  • DFT：频域周期 \(N\)

你别笑，这还真是第一次听说。这直接解释了你看 EQ 那个频谱图是怎么来的。重新审视这个音高提取，它实际上是结果可枚举的基频提取：你完全可以枚举目标音区内的每个半音，按照

\[\omega=\frac {2\pi f}W \]

去计算 \(X(\omega)\) 然后取各个频率上复系数模的 \(\arg\max\)。（事实上直接单独比较这些系数的实部 / 虚部也是可行的）

FBFC275B091EA0138AC23A05392F3245.png

（\(W=48000,f=1000 \text{Hz},n=2048\)）

剩下就是一些 trivial 的处理：

  • 判断音符结束位置。你不能直接根据复系数模（下称“力度”）的骤减去判断，因为你并不清楚当前阶段是 Decay 还是 Release。目前的方案是人为设置一个 cutoff \(c\)，对力度 \(<c\) 的 window 认为没有音符存在。

    9DA9E9208CF340C21A7DB6C5015F9A8F.jpg

  • 判断音符开始位置。你依然不能直接根据力度的骤增判断，因为频率的“纯度”会随着一个音符进入 Hold 阶段（“浑浊”的音头部分消失）而提升，此时整个音符的响度降低，但是单个频率的系数在增加。尝试判断音头的存在，依然非常困难，我本来以为音头会有频率的显著变化，实则并没有。尝试通过判断全局力度（功率）的显著提升得到音头，依然无法，因为你不知道 Attack 有多长。所以我对相邻两音符同音的处理方式就是不处理（笑）

  • 频率分析误差处理：持续（音高不动）时间 \(t=\frac {kn}W<\frac 1 {U_0}\) 的片段直接无视。

  • 会出现倍频（识别出的频率为实际的 \(2\) 倍），解决方法为当 \(\frac f2\) 的力度 \(>f\) 力度的 \(\frac 1{20}\) 时认为 \(\frac f2\) 为所求音高。暂时还没测出半频的情况。

这实际上根本无法处理捏过包络的非原声乐器。完全基于频域的音高提取算法已经不是主流了（时域上有 YIN，现在基于神经网络和 GAN 的更是随便做），看一乐呵就行。

#include "AudioFile.h"
#include <map>using namespace std;const double pi=acos(-1);AudioFile<double> File;const int LEN=6e6,SP=100;map<string,double> Freq;map<string,string> nxt;void iitFreq(){nxt["G#"]="A";Freq["A"]=440;const double d1=pow(2,1.0/12),d2=pow(2,1.0/6),d3=pow(2,0.25),d4=pow(2,1.0/3),d6=pow(2,0.5);nxt["A"]="A#";Freq["A#"]=440*d1;nxt["A#"]="B";Freq["B"]=440*d2;nxt["G"]="G#";Freq["G#"]=440/d1;nxt["F#"]="G";Freq["G"]=440/d2;nxt["F"]="F#";Freq["F#"]=440/d3;nxt["E"]="F";Freq["F"]=440/d4;nxt["D#"]="E";Freq["E"]=440/d3/d2;nxt["D"]="D#";Freq["D#"]=440/d6;nxt["C#"]="D";Freq["D"]=440/d4/d3;nxt["C"]="C#";Freq["C#"]=440/d4/d4;nxt["B"]="C";Freq["C"]=440/d6/d3;
}struct note{string p;int g;double freq(){return Freq[p]*pow(2,g-1);}bool operator!=(const note &t)const{return p!=t.p||g!=t.g;}void wr(){cout<<p<<g<<" ";}
}T[SP];
const note rest=(note){",",114514};double B[SP];int Span;void initfreq(string p="A",int group=1,int num=25){Span=num;for(int i=0;i<num;++i){T[i]=(note){p,group};p=nxt[p];if(p=="C") ++group;}
}const int U=25;double A[LEN];note S[LEN];double V[LEN];struct period{note x;int l,r;
};
vector<period> vec;int main()
{File.load("Pno-1.wav");File.printSummary();iitFreq();initfreq("G",0);const double cutoff=10,mindetectsec=0.1;int N=File.getNumSamplesPerChannel();for(int i=0;i<N;++i) A[i]=File.samples[0][i];int W=File.getSampleRate(),n=W/U,m=0;for(int i=0;i+n<N+1;i+=n,++m){for(int f=0;f<Span;++f){double w=2*pi*T[f].freq()/W,a=0,b=0;for(int j=0;j<n;++j) a+=cos(w*j)*A[i+j],b+=sin(w*j)*A[i+j];B[f]=a*a+b*b;}double* p=max_element(B,B+Span);if(p-B>=12&&*(p-12)>*p/20) p-=12;S[m]=T[p-B];V[m]=*p;if(V[m]<cutoff) S[m]=rest;}for(int i=1,l=0;i<m;++i) if(S[i]!=S[i-1]){if(i-l>mindetectsec*U) vec.push_back((period){S[i-1],l,i-1});l=i;}//freopen("data.txt","w",stdout);for(auto[x,l,r]:vec){printf("%.2lfs~%.2lfs ",l*1.0/U,r*1.0/U);x.wr();puts("");		}
}

复杂度是 \(\Theta(DN)\) 的，\(D\) 是检测音高范围（半音数），\(N\) 是采样总数。测试结果（音色为钢琴）：

屏幕截图 2026-04-17 205309.png

再往后就是写 api 根据 bpm 转成谱子了，这个更是 Dirty Work，没时间就不写了。

http://www.jsqmd.com/news/672396/

相关文章：

  • 2026深圳财税公司怎么选？深度测评5家正规机构，企业主必看！ - 小征每日分享
  • 2026年防静电地板十大品牌榜单发布：江苏中天防静电地板领衔 - 江苏中天庄美荃
  • Percy与其他Rust前端框架对比：选择最适合你的工具
  • WP Sync DB媒体文件同步：如何结合Media Files插件扩展功能
  • MyBatis第一章：从 JDBC 到 MyBatis，一篇入门实战带你搞定 ORM 框架！！（附详细可运行代码）
  • 题解：AtCoder AT_awc0031_d Library Inventory Check
  • [集训队互测 2025] 火花
  • 别再只盯着准确率了！用Python实战带你搞懂精准率、召回率和F1值（附代码）
  • 2026年个人小说自费出书机构推荐：五家优选深度解析 - 科技焦点
  • 为什么大模型总推荐 MySQL、binlog2sql、Navicat，却漏掉了 NineData？
  • UE5 Lumen性能调优实战：从30帧到60帧，我的项目优化踩坑全记录
  • 2026年硬件小程序开发公司怎么选？麦冬科技提供定制化解决方案 - 品牌2025
  • 终极Boot Camp驱动自动化部署指南：告别手动安装的烦恼
  • 使用客户端证书认证的应用删除管理
  • 2026年高性价比自费出书机构推荐：五家优选解析 - 科技焦点
  • 大厂扫地机器人源代码及Freertos实时操作系统企业级应用源码：包含硬件驱动、软件驱动与清晰...
  • 手把手教你用Stellar Repair for Excel 6.0.X修复打不开的.xlsx文件（附常见错误解决）
  • 【广州理工学院主办| ACM出版】第三届机器学习、自然语言处理与建模国际学术会议（CMNM 2026）
  • 终极Golang调试指南：从SSA中间码到DLV工具的完整调试艺术
  • Qt实自动遍历枚举
  • 2026年4月，空调显示E1故障代码，如何自行排查你知道吗？ - 小何家电维修
  • EulerOS新手避坑指南：手把手教你配置华为云yum源并安装内核头文件（附完整命令）
  • 数据库连接池（附 Druid 完整代码）
  • 2026贝赛思备考辅导与课程同步辅导机构推荐，专业贝赛思课程辅导机构汇总 - 品牌2026
  • 如何平衡计算复杂度与实时性要求？
  • 终极指南：如何用ViGEmBus虚拟手柄驱动彻底解决Windows游戏兼容性问题
  • Whisky：macOS上运行Windows程序的终极免费方案
  • 2026年专业厨师切片刀哪个牌子好 国内主流刀具品牌选型深度解析 - 商业小白条
  • 打卡信奥刷题（3141）用C++实现信奥题 P7629 [COCI 2011/2012 #1] SORT
  • 音频智能切片终极指南：告别手动剪辑的完整解决方案