【Delphi】如何在三轴加速器的频谱分析中使用FFT(快速傅里叶变换)算法
关于傅里叶变换的作用，网上说的太过学术化，且都在说原理，以及如何编码实现，可能很多人有个模糊印象，在人工智能，图像识别，运动分析，机器学习等中，频谱分析成为了必备的手段，可将离散信号量转换为数字信息进行归类分析。
今天这里讲的不是如何实现，而是如何使用傅里叶变换
但频谱分析中，涉及到的信号处理知识对大部分软件开发的人来说，太过于晦涩难懂，傅里叶变换，拉普拉斯，卷积，模相，实数，虚数，复数，三角函数等等，已经能让软件工程师望而却步，造成懂知识的人无法开发，懂开发的人无法分析，而同时具备两种技能的人，除了专业的研究生博士，剩下的少之又少。
站在软件开发人员的角度，能否抛开这些晦涩难懂概念，进行纯软件方面的信号处理？这个问题可能没有答案，最好能够通过实践来证明。而且若是抛开这些概念，让那些有兴趣的开发人员学习到信号处理，频谱分析，先不管可不可能，光想想，其作用也是有的，至少“高深”的东西能够被用上了。
三轴加速器步数计算分析
下面就以三轴加速器（X，Y，Z三个加速器）的实际应用开发为例，来讲一讲（限于个人水平，可能比较粗略）。
时下运动APP也有不少，不少APP开发人员就会接触到陀螺仪，加速器等，抛开计步的准确性不讲，我们这里来看看这些计步分析的原理，如何利用这些传感器进行行为分析，首先需要获取传感器的数据（如何获取就不在这里讲了，有Android，iOS，独立蓝牙/WIFI设备等），将获取的数据输出到平面坐标图表中：
如上图，即为时域坐标曲线图，X,Y,Z三种颜色分别代表三轴加速器的三个方向加速度随时间的变化曲线。坐标的横轴为时间轴，纵轴为加速度值（该值由加速度传感器的电压值转换而来，具体由硬件开发者决定）。
我们的采集频率为50Hz，即横轴（时间轴）每一个点为50分之1秒，可以看到前4秒（0-200）没有太明显走动，4秒后到12秒走动速度较慢，12秒走路稍微加快。在这里不得不感叹人脑简直是宇宙超级无敌，光看图表就能够直接知道行为，无法想象未来谁能够让人工智能超越人脑。
在软件开发中，我们要如何通过代码分析得到上面“人脑的分析结果“？ 要直接上代码不是我的作风，百度可以搜到一篇腾讯的文章《利用三轴加速器的计步测算方法》，此文章的分析方法可以说大致能够理解，但对于读者可能帮助不是特别大，因为涉及到的计算方法却语焉不详，比如，在没有进行滤波前，上图的X轴波峰波谷并不平滑，可以看出每个大波峰上又有2个小波谷，如果以此来计算，得出来的步数误差非常大，如何解决误差？
可能一些做过信号处理的人会说，做一个“均值滚动滤波+低通滤波+高通滤波”就差不多了，我只能说对于复杂的“周期性”信号分析是可行的，对于三轴加速器的步数计算来说，实际中并不存在什么周期性信号，人的上坡下坡走路非直线，时快时慢等动作都会加大分析难度，往往理论说起来容易，实现起来难，会发现实际情况和理论不太一样。
先不说这么多了，拿到信号原始数据，马上进行傅里叶变换，Delphi可以使用FPC一个开源的库AlgLib，最新版好像不开源了且分为免费版和收费版，但codetyphon收录了以前的开源版本可以用。
alglib中关于快速傅里叶变换是在u_fft.pp文件中，delphi只把后缀pp改成pas就可以使用了。使用复数FFTC1D和实数FFTR1D函数进行变换，这里仅做单轴分析，所以使用实数FFTR1D函数就好了，复数傅里叶变化应用在三维波动曲线的分析中，这里不用就不关心了（针对运动曲线来讲，至于无线电信号的复合器或者多维矩阵数据分析的情况就多了），2个函数中的N即为采集的数据量，当然实际中我们有时候采集了几十秒，而N的最佳范围是在128-2048范围（既充分又计算量不高且做代码运算方便），所以我们每10秒（50Hz即500个数据，最好是直接扩展到512个。2的n次方）进行一次傅里叶变换分析即可。
题外话：
专业人士看到这里可能会特别说明FFT仅仅是离散傅里叶变换(DFT)，由于计算机是由时间片定时执行一次的，所以采集的数据也是每隔一段时间采集一次，当然间隔时间越短越贴近连续采集，总的来讲，就是非连续的，离散的数据，计算机也只能做离散傅里叶变换分析了，电子电路设备光学设备等才能做连续傅里叶变换。
至于为何使用实数傅里叶变换，因为我们采集的数据是一个一个采集的（三轴实际上是分开独立来分析的），所以采集的数据只是一维（时间维度不计在内），使用实数傅里叶变换再适合不过了。
FFTR1D函数的调用方法
比如采集到的数据为[3.4, 4.1, 6.7, 3.1...],则如下代码（仅参考用）：
var
A : TComplex1DArray;
R : TReal1DArray;
I, N : Integer;
begin
N := 50*10 // 10秒（建议使用512，2的N次方）
SetLength(R, N);
for I:=0 to N-1 do
begin // 时域的复合幅度值
R[I]:=Data[nIdx];//按顺序填上采集的数据[3.4, 4.1, 6.7, 3.1...]
end;
FFTR1D(R, N, A); // 使用实数傅里叶变换方法（速度理论上比复数快一倍）
//////////////////////////
// 到这里已经完成傅里叶变换，如何使用变换后的数据？？
end;
复数数组A即为傅里叶变换结果，但是这个结果有啥用？这里先引用百科上的一段话来帮助理解后面如何利用傅里叶变换结果的代码：
假设FFT之后某点n用复数a+bi表示， 那么这个复数的模就是An=根号a*a+b*b， 相位就是Pn=atan2(b,a)。 根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。 对于n=1点的信号，是直流分量，幅度即为A1/N。 由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。
经过变换后的复数数组A，即是变换结果，但并不是分析的数据，而是一个中间数据，由该数据可以得到模值和相位，根据傅里叶变换原理，模值计算如下：
// 对A数组每个点的复数求模，即实部和虚部平方和再开根号。 nAValue := Math.Power(Abs(A[I].X*A[I].X)+Abs(A[I].Y*A[I].Y), 0.5);
至于相位则是根据公式Pn=atan2(b,a)计算，但是频谱分析比较有用的是模值（主要是滤波，这里区别于物理硬件上的滤波-采用共振/折射等方法），所以其他不管。
前面讲到使用10秒的数据进行计算，这样得到的傅里叶变换结果复数也有500个，由于傅里叶的对称性，我们只取0-250的结果来计算模值即可，因为251-500是对称的结果。
计算模值后再输出平面坐标图表（频域图），如下：
当然幅度也有作用，幅度的计算说明如下（来自百度百科）：
假设原始信号中某个频率的峰值（幅度）为V，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是V的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。
计算幅度 = 模值*2/N（第一个点= 模值/N），后再输出平面坐标图表（频域图），如下：
如上图，网络上各种傅里叶变换文章中的标志性的结果图就出来了（请忽略图中的时间刻度与文章的10秒500个数据等相关性，因为我用的图是随便一段数据来的。）
看到上面的傅里叶变换结果图，会不会觉得悲剧才开始发生。没错，结果是出来了，但我们该如何分析呢？
如何分析傅里叶变换结果
首先傅里叶变换后的结果要看得懂，直白的讲，变换后的结果反应了各个频率的模值，模值越高代表该频率产生的“作用”越大，也就代表越切合实际存在该频率的信号在起作用，而图中的纵坐标就是模值（或幅度），横坐标就是频率，但频率是整数的，实际频率与前面的采集频率和数据量有关，关系如下：
若采样频率是50Hz，采集50个数据说明我们能够分辨出1Hz的频率，可以简单认为分辨率为1，而500个点代表频率的分辨率是10，也说明这个傅里叶变换结果能够分辨出0.1Hz的频率，在这里代表每一个横坐标的整数点对应的频率刻度为0.1Hz。
所以看到上图15-20范围内的X加速器（蓝色线）的模值很高，代表1.5-2Hz的频率最切合实际，说明在这一部分的分析结果人的走路速度是1秒钟1.5-2步（走得比较慢），而我们的业余半马运动员步频在180步/分钟左右，专业的190步/分钟左右。
继续话题，实际中我们开发的软件是需要自动化处理，不可能跟人脑一样高级，一看图就知道结果，当然如果靠人脑，傅里叶变换都不需要就可以分析了，所以别想那么多了，在人工智能无法代替人脑的时代，软件的信号处理，到这里才刚刚开始。
如何进行下一步分析呢，我们可以取一个范围，如1-5Hz内的模值数据，其他的滤波处理掉，同时根据采样频率的二分之一为带宽(奈奎斯特频率)之间的关系，50Hz的采样频率必须把高于25Hz的频段滤掉（幸好人行走最高5Hz，非要跑到6Hz以上的疯子就不管了吧），得到简化后的模值复数数组。
另外，小幅度波动也可以滤掉，至于幅度值多少可以去掉，可以自行统计判断，可能不同加速度出来的值也不一样，简单的按照最大幅度的百分比算，低于10%的全部滤掉。不过为了后面可能需要的积分计算准确性更高，保留着也是可以的。
利用傅里叶反变换进行滤波
（未完待续）
后续还有如何均值化，方波化等等，最终得到非常贴近实际的步数结果
