http://blog.sina.com.cn/s/blog_6163bdeb0102ebi5.html
最近做有關加速度的數據處理，需要把加速度積分成位移，網上找了找相關資料，發現做這個並不多，把最近做的總結一下吧！
積分操作主要有兩種方法：時域積分和頻域積分，積分中常見的問題就是會產生二次趨勢。關於積分的方法，在國外一個論壇上有人提出了如下說法，供參考。
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函數實現積分，它是封裝之后的函數，可以實現時域積分和頻域積分，其代碼如下
% 積分操作由加速度求位移，可選時域積分和頻域積分
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
其中時域積分的子函數如下
% 時域內梯形積分
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
頻域積分的子函數如下（此代碼是一個老外編的，在頻域內實現積分和微分操作）
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