加速度频域积分的实现及其局限性分析

首先，相信大家都尝试过直接在时域中通过加速度传感器积分得到位移。在加速度精度不高或者加速度数据不经处理的情况下，积分得到的位移量会一直有一个累计误差，而且会越来越大，这时有人就会把目光移到频域中，在频域中对加速度进行积分会怎样呢？会不会有出乎意料的效果呢？

什么是频域积分

单片机或者传感器采样得到的点都是离散的，在时域中，对于离散点的积分就是求和。

而频域积分需要先把时域的数组通过快速傅里叶变换（FFT）转到频域，再通过傅里叶变换的积分性质，在频域中实现积分与滤波0后，通过IFFT傅里叶逆变换返回时域，此时得到的数组即为时域数组的积分结果。傅里叶变换的原理此处不再赘述，大概的作用就是把一个时域上的信号波形分解成若干个单一频率信号的叠加，通俗点讲就是通过无穷多个幅度、频率不同的正弦波的叠加去近似原信号的波形，而傅里叶变换得到的频谱就是这些不同频率正弦波形对应的幅值。这里丢两张动图方便理解。

而频域的积分需要用到傅里叶变换的积分性质：

一次积分

二次积分

这两个公式将时域积分放到了到频域做，去掉了积分符号，看起来比较复杂，其实如果以不变量一个变量地看都比较简单。

如何实现频域积分

用程序语言描述傅里叶变换的积分性质，看似复杂，其实只需要把这个公式每一项单独求出即可。

下面步骤是参考了《matlab在振动信号处理中的应用》（王济著）中第五章频域积分。

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。ω数组计算公式如下，其中n为积分次数：

一次积分的ω数组的与二次积分的ω数组如图所示：

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语言代码点击下面链接

基于STM32F4的加速度频域二次积分振动位移C语言算法

基于STM32F407与JY901模块的加速度频域积分实现（完整工程）

clear;
clc;
close all hidden;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fno = 'output.txt';
fid = fopen('1.txt','r');
x = fscanf(fid,'%f',[1,inf]);
status = fclose(fid);sf = 1000;            %采样频率
fmin = 1;           %最小截止频率
fmax = 200;           %最大截止频率
c = 9810.5;           %单位变换系数
it = 2;               %积分次数sx = '时间';          %横坐标标注
sy1 = 'y1';           %纵坐标标注
sy2 = 'y2';           %纵坐标标注n=length(x);                    %待分析数组长度
t=0:1/sf:(n-1)/sf;              %建立时间向量
nfft=2^nextpow2(n);             %大于并接近n的2的幂次方的FFT长度
y=fft(x,nfft);                  %FFT变换df=sf/nfft;                     %计算频率间隔（Hz/s）
ni=round(fmin/df+0.5);          %高通频率截止点=最小截止频率/时间步长t+1
na=round(fmax/df+0.5);          %低通频率截止点=最大截止频率/时间步长t+1
dw=2*pi*df;                     %圆频率间隔（rad/s）
w1=0:dw:2*pi*(0.5*sf-df);       %正离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;    %负离散圆频率向量
w=[w1,w2];                      %连接w1,w2
w=w.^it;                        %以积分次数为指数，建立圆频率变量向量
subplot(2,2,2);
plot(abs(y));a=zeros(1,nfft);                %进行积分的频域变换
a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1);
subplot(2,2,4);                 %绘制积分前的时程曲线图形
plot(abs(a));if it==2y=-a;                       %进行二次积分的相位变换
elsey = imag(a) - 1i*real(a);   %进行一次积分的相位变换
enda=zeros(1,na);                  
a(ni:na)=y(ni:na);              %消除指定正频带外的频率成分
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);%消除指定负频带外的频率成分y=ifft(a,nfft);                 %IFFT变换
y=real(y(1:n))*c;               %取逆变换的实部n个元素并乘以单位变换系数为积分结果
subplot(2,2,1);                 %绘制积分前的时程曲线图形
plot(t,x);
xlabel('时间(s)');
ylabel('加速度(g)');
grid on;
subplot(2,2,3);                 %绘制积分前的时程曲线图形
plot(t,y);
xlabel('时间(s)');
if it==2ylabel('距离(m)');
elseylabel('速度(m/s)');
end
grid on;
fid=fopen(fno,'w');             %输出积分后的数据
for k=1:nfprintf(fid,'%f/n',y(k));
end
status=fclose(fid);

局限性

但是上面的计算思路仅局限于振动的位移计算，也就是说总位移为0，当总位移不为0时，我们的积分结果最终还是会等于0，下面看一下仿真结果

图一是原始数据，图三为图一的积分，积分结果不为0，图二为FFT后去除最低频信号（即直流部分）后再IFFT的波形，可以看到整个图形都向下偏移了一段距离，图四为图二的积分结果，最终为0，这与实际积分结果严重不符。因为根据傅里叶变换公式 $F(\omega )=\int_{-\propto }^{+\propto }f(t)e^{-i\omega t}dt$ 当ω=0时， $F(0)=\int_{-\propto }^{+\propto }f(t)dt$ ，可见F(0)本身就等于时域图形的积分，而滤低频时将其直接置为0，最后的积分结果当然为0，也就是说上述频域积分仅适用于总积分为0的情况。

若滤除直流部分，最终积分为0，若不滤除直流分量部分，又无法消除加速度积分后的累计误差，因此互相矛盾。所以频域积分不适用于总位移不为0的加速度积分。

非零位移

关于非零位移的问题，首先我们要清楚加速度计累积误差的来源，加速度计都带有一定的线性偏移，即需要通过kx+b的方法进行校正，而其中的b正是引起加速度积分累积误差的罪魁祸首，因此需要对加速度进行线性补偿，可通过Matlab中的lsqcurvefit，非线性拟合文中提到的方法确定k和b，如果在要求不是很严格的场合，只需要对b进行补偿即可，方法是在确定系统处于静止状态时（通过其他方式确定的系统状态），首先令速度为零，然后将一段时间内的加速度取平均即可作为补偿值b，每隔适当时间校正一次b这样就能保证速度基本无累积误差，然后用同样的方法计算位移。

但是上述方法仍然有比较大的误差存在，仅加速度计显然无法精确求得位移，因此需要增加其他传感器。对于在平面上运动的系统，建议增加GPS，如果是在垂直方向上运动的系统则可以增加气压计，然后通过卡尔曼滤波对两组数据进行数据融合，能得到比较准确的位移、速度和加速的。具体方法可以参考这篇博文：二阶卡尔曼滤波计算加速度、速度及高度