加载中…加速度频域积分的实现及其局限性分析
(2020-05-11 12:34:29)
标签：
it健康美食文化杂谈 |
分类： 硬件交流 |
首先，相信大家都尝试过直接在时域中通过加速度传感器积分得到位移。在加速度精度不高或者加速度数据不经处理的情况下，积分得到的位移量会一直有一个累计误差，而且会越来越大，这时有人就会把目光移到频域中，在频域中对加速度进行积分会怎样呢？会不会有出乎意料的效果呢？
而频域积分需要先把时域的数组通过快速傅里叶变换（FFT）转到频域，再通过傅里叶变换的积分性质，在频域中实现积分与滤波0后，通过IFFT傅里叶逆变换返回时域，此时得到的数组即为时域数组的积分结果。傅里叶变换的原理此处不再赘述，大概的作用就是把一个时域上的信号波形分解成若干个单一频率信号的叠加，通俗点讲就是通过无穷多个幅度、频率不同的正弦波的叠加去近似原信号的波形，而傅里叶变换得到的频谱就是这些不同频率正弦波形对应的幅值
1. 获得时域数组
假设加速度的采样频率sf为1kHz，则需要先采集一段时间的数据得到时域数组，为FFT计算方便起见，nfft取2的幂次个采样点如1024个点；
2.计算FFT
对这1024个点做FFT快速傅里叶变换得到频谱，即上面公式中的F(ω)，将F(ω)取模长后我们会发现该图形是个关于500Hz对称的图形（第一个点为F(0)，不对称）
3. 计算ω数组
要获得ω数组首先需要算出频率间隔df，因为采样频sf率是1000Hz，而整个数组长度为1024，所以df=sf/nff，单位(Hz/s)，通过计算df把时域中采样的时间点和频域中的频率点一一对应起来，再说直白点就是时域中采样点的时间间隔是1/1000Hz=1ms，而频域中的的频率间隔是1000/1024=0.977Hz，我们FFT得到的频谱每个值之间间隔0.977Hz，1024个点对应0~1000Hz。ω数组也是对称的（或中心对称），分为正离散圆频率向量和负离散圆频率向量。先将刚才得到的频率间隔df转换为角频率间隔dω，dω=2π*df。
4. 虚数 j
得到ω数组后就是处理 j 了，j 是单位虚数，j的平方为-1。积分公式中F(ω)除了一个 j，而 j 在频域中表示相移，每乘一个 j 就逆时针旋转了90度，每除一个j，顺时针旋转90度。也就是说，时域的一次积分转到频域后需要顺时针旋转90度，而二次积分需要顺时针旋转180度。
5. 进行积分的频域变换
频域变换就是第二步中FFT的结果（即1024个复数）依次除以ω数组。
6. 进行积分的相位变换
第4步中提到一次积分顺时针旋转90度，而二次积分顺时针旋转180度。相位变换就是把第5步得到的复数数组每个复数都旋转，假设实部为real，虚部为imag，一次积分时虚部等于实部imag=real，实部等于负虚部real=-imag；二次积分时实部等于负实部real=-real，虚部等于负虚部imag=-imag
7. 滤波
此时如果不进行滤波直接跳到第8步就能得到时域的二次积分结果。如果想把其中的低频和高频部分滤除，则可以将第6步得到的复数数组中的开头几个复数和最后几个复数直接置零，假设需要滤除Fmin=10Hz以下的部分，先计算出10Hz大概是哪个点ni=Fim/df+1，10/0.977+1≈11，即把第6步得到的复数数组前11个复数都置位0，滤高频同理。
补充：以上方法其实就是用了矩形窗函数，如果想要达到不同的滤波效果，可以考虑其他窗函数如汉宁窗、海明窗、高斯窗等等。
8. IFFT返回时域
最后通过IFFT将处理完的复数数组转回时域，得到积分结果。
附上MATLAB代码，C语言代码点击下面链接
喜欢
0
赠金笔