積分操作主要有兩種方法：時域積分和頻域積分，積分中常見的問題就是會產生二次趨勢。關于積分的方法，在國外一個論壇上有人提出了如下說法，供參考。
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.
- Remove the mean from your sample (now have zero-mean sample)
- Integrate once to get velocity using some rule (trapezoidal, etc.)
- Remove the mean from the velocity
- Integrate again to get displacement.
- Remove the mean. Note, if you plot this, you will see drift over time.
- To eliminate (some to most) of the drift (trend), use a least squares fit (high degree depending on data) to determine polynomial coefficients.
- 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…
- Remove the mean from the accel. data
- Take the Fourier transform (FFT) of the accel. data.
- Convert the transformed accel. data to displacement data by dividing each element by -omega^2, where omega is the frequency band.
- 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是這樣的效果）。分析兩種方法的頻率特性曲線如下：
時域積分
頻域積分
對于包含兩個正弦波的信號，頻域積分正常恢復信號，時域積分恢復的高頻信息有誤差；對于有噪聲的正弦信號，噪聲會使積分結果產生大的趨勢項（不是簡單的二次趨勢），如下圖：
對此可以用濾波的方法將大的趨勢項去掉。測試的代碼如下：
<section class="" data-source="bj.96weixin.com">
<blockquote class="">% 測試積分對正弦信號的作用
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
end
close
figure
plot(f, amp)
title('位移積分的頻率特性曲線')
xlabel('f')
ylabel('單位正弦波的積分位移幅值')</blockquote>
</section>
以上代碼中使用IntFcn函數實現積分，它是封裝之后的函數，可以實現時域積分和頻域積分，其代碼如下：
<section class="" data-source="bj.96weixin.com">
<blockquote class="">% 積分操作由加速度求位移，可選時域積分和頻域積分
function [disint, velint] = IntFcn(acc, t, ts, flag)
if flag == 1
% 時域積分
[disint, velint] = IntFcn_Time(t, acc);
velenergy = sqrt(sum(velint.^2));
velint = detrend(velint);
velreenergy = sqrt(sum(velint.^2));
velint = velint/velreenergy*velenergy;
disenergy = sqrt(sum(disint.^2));
disint = detrend(disint);
disreenergy = sqrt(sum(disint.^2));
disint = disint/disreenergy*disenergy;
% 此操作是為了彌補去趨勢時能量的損失
% 去除位移中的二次項
p = polyfit(t, disint, 2);
disint = disint – polyval(p, t);
else
% 頻域積分
velint = ?iomega(acc, ts, 3, 2);
velint = detrend(velint);
disint = ?iomega(acc, ts, 3, 1);
% 去除位移中的二次項
p = polyfit(t, disint, 2);
disint = disint – polyval(p, t);
end
end</blockquote>
</section>
其中時域積分的子函數如下：
<section class="" data-source="bj.96weixin.com">
<blockquote class="">% 時域內梯形積分
function [xn, vn] = IntFcn_Time(t, an)
vn = cumtrapz(t, an);
vn = vn – repmat(mean(vn), size(vn,1), 1);
xn = cumtrapz(t, vn);
xn = xn – repmat(mean(xn), size(xn,1), 1);
end</blockquote>
</section>
頻域積分的子函數如下（此代碼是一個老外編的，在頻域內實現積分和微分操作）
<section class="" data-source="bj.96weixin.com">
<blockquote class="">function dataout = ?iomega(datain, dt, datain_type, dataout_type)
%%%%%%%%%%%%%%%%
%
% ? IOMEGA is a MATLAB script for converting displacement, velocity, or
% ? acceleration time-series to either displacement, velocity, or
% ? acceleration times-series. The script takes an array of waveform data
% ? (datain), transforms into the frequency-domain in order to more easily
% ? convert into desired output form, and then converts back into the time
% ? domain resulting in output (dataout) that is converted into the desired
% ? form.
%
% ? Variables:
% ? ———-
%
% ? datain ? ? ? = ? input waveform data of type datain_type
%
% ? dataout ? ? ?= ? output waveform data of type dataout_type
%
% ? dt ? ? ? ? ? = ? time increment (units of seconds per sample)
%
% ? ? ? ? ? ? ? ? ? ?1 – Displacement
% ? datain_type ?= ? 2 – Velocity
% ? ? ? ? ? ? ? ? ? ?3 – Acceleration
%
% ? ? ? ? ? ? ? ? ? ?1 – Displacement
% ? dataout_type = ? 2 – Velocity
% ? ? ? ? ? ? ? ? ? ?3 – Acceleration
%
%
%%%%%%%%%%%%%%%%
% ? Make sure that datain_type and dataout_type are either 1, 2 or 3
if (datain_type < 1 || datain_type > 3)
error('Value for datain_type must be a 1, 2 or 3');
elseif (dataout_type < 1 || dataout_type > 3)
error('Value for dataout_type must be a 1, 2 or 3');
end
% ? Determine Number of points (next power of 2), frequency increment
% ? and Nyquist frequency
N = 2^nextpow2(max(size(datain)));
df = 1/(N*dt);
Nyq = 1/(2*dt);
% ? Save frequency array
iomega_array = 1i*2*pi*(-Nyq : df : Nyq-df);
iomega_exp = dataout_type – datain_type;
% ? Pad datain array with zeros (if needed)
size1 = size(datain,1);
size2 = size(datain,2);
if (N-size1 ~= 0 && N-size2 ~= 0)
if size1 > size2
datain = vertcat(datain,zeros(N-size1,1));
else
datain = horzcat(datain,zeros(1,N-size2));
end
end
% ? Transform datain into frequency domain via FFT and shift output (A)
% ? so that zero-frequency amplitude is in the middle of the array
% ? (instead of the beginning)
A = fft(datain);
A = fftshift(A);
% ? Convert datain of type datain_type to type dataout_type
for j = 1 : N
if iomega_array(j) ~= 0
A(j) = A(j) * (iomega_array(j) ^ iomega_exp);
else
A(j) = complex(0.0,0.0);
end
end
% ? Shift new frequency-amplitude array back to MATLAB format and
% ? transform back into the time domain via the inverse FFT.
A = ifftshift(A);
datain = ifft(A);
% ? Remove zeros that were added to datain in order to pad to next
% ? biggerst power of 2 and return dataout.
if size1 > size2
dataout = real(datain(1:size1,size2));
else
dataout = real(datain(size1,1:size2));
end
return</blockquote>
</section>
杭州銳達微信二維碼