【Delphi】如何在三轴加速器的频谱分析中使用FFT(快速傅里叶变换)算法

时间 2019-11-12

标签 Delphi 如何频谱分析使用 fft 快速傅里叶变换算法栏目 Delphi 繁體版

原文原文链接

关于傅里叶变换的做用，网上说的太过学术化，且都在说原理，以及如何编码实现，可能不少人有个模糊印象，在人工智能，图像识别，运动分析，机器学习等中，频谱分析成为了必备的手段，可将离散信号量转换为数字信息进行归类分析。数组

今天这里讲的不是如何实现，而是如何使用傅里叶变换

但频谱分析中，涉及到的信号处理知识对大部分软件开发的人来讲，太过于晦涩难懂，傅里叶变换，拉普拉斯，卷积，模相，实数，虚数，复数，三角函数等等，已经能让软件工程师望而却步，形成懂知识的人没法开发，懂开发的人没法分析，而同时具有两种技能的人，除了专业的研究生博士，剩下的少之又少。网络

站在软件开发人员的角度，可否抛开这些晦涩难懂概念，进行纯软件方面的信号处理？这个问题可能没有答案，最好可以经过实践来证实。并且如果抛开这些概念，让那些有兴趣的开发人员学习到信号处理，频谱分析，先无论可不可能，光想一想，其做用也是有的，至少“高深”的东西可以被用上了。机器学习

三轴加速器步数计算分析

下面就以三轴加速器（X，Y，Z三个加速器）的实际应用开发为例，来说一讲（限于我的水平，可能比较粗略）。函数

时下运动APP也有很多，很多APP开发人员就会接触到陀螺仪，加速器等，抛开计步的准确性不讲，咱们这里来看看这些计步分析的原理，如何利用这些传感器进行行为分析，首先须要获取传感器的数据（如何获取就不在这里讲了，有Android，iOS，独立蓝牙/WIFI设备等），将获取的数据输出到平面坐标图表中：学习

如上图，即为时域坐标曲线图，X,Y,Z三种颜色分别表明三轴加速器的三个方向加速度随时间的变化曲线。坐标的横轴为时间轴，纵轴为加速度值（该值由加速度传感器的电压值转换而来，具体由硬件开发者决定）。编码

咱们的采集频率为50Hz，即横轴（时间轴）每个点为50分之1秒，能够看到前4秒（0-200）没有太明显走动，4秒后到12秒走动速度较慢，12秒走路稍微加快。在这里不得不感叹人脑简直是宇宙超级无敌，光看图表就可以直接知道行为，没法想象将来谁可以让人工智能超越人脑。人工智能

在软件开发中，咱们要如何经过代码分析获得上面“人脑的分析结果“？要直接上代码不是个人做风，百度能够搜到一篇腾讯的文章《利用三轴加速器的计步测算方法》，此文章的分析方法能够说大体可以理解，但对于读者可能帮助不是特别大，由于涉及到的计算方法却语焉不详，好比，在没有进行滤波前，上图的X轴波峰波谷并不平滑，能够看出每一个大波峰上又有2个小波谷，若是以此来计算，得出来的步数偏差很是大，如何解决偏差？spa

可能一些作过信号处理的人会说，作一个“均值滚动滤波+低通滤波+高通滤波”就差很少了，我只能说对于复杂的“周期性”信号分析是可行的，对于三轴加速器的步数计算来讲，实际中并不存在什么周期性信号，人的上坡下坡走路非直线，时快时慢等动做都会加大分析难度，每每理论提及来容易，实现起来难，会发现实际状况和理论不太同样。
code

先不说这么多了，拿到信号原始数据，立刻进行傅里叶变换，Delphi可使用FPC一个开源的库AlgLib，最新版好像不开源了且分为免费版和收费版，但codetyphon收录了之前的开源版本能够用。blog

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%的所有滤掉。不过为了后面可能须要的积分计算准确性更高，保留着也是能够的。

利用傅里叶反变换进行滤波

（未完待续）

后续还有如何均值化，方波化等等，最终获得很是贴近实际的步数结果