加速度位移积分

加速度积分位移 Matlab

2013-02-04 05:30:00| 分类：MATLAB应用|举报|字号订阅

最近做有关加速度的数据处理，需要把加速度积分成位移，网上找了找相关资料，发现做这个并不多，把最近做的总结一下吧！

积分操作主要有两种方法：时域积分和频域积分，积分中常见的问题就是会产生二次趋势。关于积分的方法，在国外一个论坛上有人提出了如下说法，供参考。

Double integration of raw acceleration data is a pretty poor estimate for displacement. The reason is that at each integration, you are compounding the noise in the data. If you are dead set on working in the time-domain, the best results come from the following steps.

1. Remove the mean from your sample (now have zero-mean sample)

2. Integrate once to get velocity using some rule (trapezoidal, etc.)

3. Remove the mean from the velocity

4. Integrate again to get displacement.

5. Remove the mean. Note, if you plot this, you will see drift over time.

6. To eliminate (some to most) of the drift (trend), use a least squares fit (high degree depending on data) to determine polynomial coefficients.

7. Remove the least squares polynomial function from your data.

A much better way to get displacement from acceleration data is to work in the frequency domain. To do this, follow these steps...

1. Remove the mean from the accel. data

2. Take the Fourier transform (FFT) of the accel. data.

3. Convert the transformed accel. data to displacement data by dividing each element by -omega^2, where omega is the frequency band.

4. Now take the inverse FFT to get back to the time-domain and scale your result. This will give you a much better estimate of displacement.

说到底就是频域积分要比时域积分效果更好，实际测试也发现如此。原因可能是时域积分时积分一次就要去趋势，去趋势就会降低信号的能量，所以最后得到的结果常常比真实幅值要小。下面做一些测试，对一个正弦信号的二次微分做两次积分，正弦频率为50Hz，采样频率1000Hz，恢复效果如下

时域积分

频域积分

可见恢复信号都很好（对于50Hz是这样的效果）。

分析两种方法的频率特性曲线如下时域积分

频域积分

可以看到频域积分得到信号更好，时域积分随着信号频率的升高恢复的正弦幅值会降低。

对于包含两个正弦波的信号，频域积分正常恢复信号，时域积分恢复的高频信息有误差；对于有噪声的正弦信号，噪声会使积分结果产生大的趋势项（不是简单的二次趋势），如下图

对此可以用滤波的方法将大的趋势项去掉。

测试的代码如下

% 测试积分对正弦信号的作用

clc

clear

close all

%% 原始正弦信号

ts = 0.001;

fs = 1/ts;

t = 0:ts:1000*ts;

f = 50;

dis = sin(2*pi*f*t); % 位移

vel = 2*pi*f.*cos(2*pi*f*t); % 速度

acc = -(2*pi*f).^2.*sin(2*pi*f*t); % 加速度

% 多个正弦波的测试

% f1 = 400;

% dis1 = sin(2*pi*f1*t); % 位移

% vel1 = 2*pi*f1.*cos(2*pi*f1*t); % 速度

% acc1 = -(2*pi*f1).^2.*sin(2*pi*f1*t); % 加速度% dis = dis + dis1;

% vel = vel + vel1;

% acc = acc + acc1;

% 结：频域积分正常恢复信号，时域积分恢复加入的高频信息有误差

% 加噪声测试

acc = acc + (2*pi*f).^2*0.2*randn(size(acc));

% 结：噪声会使积分结果产生大的趋势项

figure

ax(1) = subplot(311);

plot(t, dis), title('位移')

ax(2) = subplot(312);

plot(t, vel), title('速度')

ax(3) = subplot(313);

plot(t, acc), title('加速度')

linkaxes(ax, 'x');

% 由加速度信号积分算位移

[disint, velint] = IntFcn(acc, t, ts, 2);

axes(ax(2)); hold on

plot(t, velint, 'r'), legend({'原始信号', '恢复信号'})

axes(ax(1)); hold on

plot(t, disint, 'r'), legend({'原始信号', '恢复信号'})

%% 测试积分算子的频率特性

n = 30;

amp = zeros(n, 1);

f = [5:30 40:10:480];

figure

for i = 1:length(f)

fi = f(i);

acc = -(2*pi*fi).^2.*sin(2*pi*fi*t); % 加速度

[disint, velint] = IntFcn(acc, t, ts, 2); % 积分算位移

amp(i) = sqrt(sum(disint.^2))/sqrt(sum(dis.^2));

plot(t, disint)

drawnow

% pause

end

close

figure

plot(f, amp)

title('位移积分的频率特性曲线')

xlabel('f')

ylabel('单位正弦波的积分位移幅值')

以上代码中使用IntFcn函数实现积分，它是封装之后的函数，可以实现时域积分和频域积分，其代码如下