加速度位移积分.doc

1、加速度积分位移Matlab2013-02-04 05:30:00|分类： MATLAB 应用 |举报 |字号订阅最近做有关加速度的数据处理，需要把加速度积分成位移，网上找了找相关资 料，发现做这个并不多，把最近做的总结一下吧！积分操作主要有两种方法：时域积分和频域积分，积分中常见的问题就是会 产生二次趋势。关于积分的方法，在国外一个论坛上有人提出了如下说法，供Double in tegrati on of raw accelerati on data is a pretty poor estimate for displaceme nt. The reason is that at each

2、integration, you arcompounding the noise in the data . If you are dead set on working in thetime-domain , the best results come from the followi ng steps.1. Remove the mean from your sample (now have zero-mea n sample)2. Integrate once to get velocity using some rule (trapezoidal, etc.)3. Remove the

3、 mean from the velocity4. Integrate again to get displacement.5. Remove the mean. Note, if you plot this, you will see drift over time.6. To elim in ate (some to most) of the drift (tre nd), use a least squares fit (high degree depe nding on data) to determ ine polyno mial coefficie nts.7. Remove th

4、e least squares polyno mial fun ctio n 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. data2. Take the Fourier tran sform (FFT) of the accel. data.3. Convert the tran sf

5、ormed accel. data to displaceme nt data by dividi ng each eleme nt by -omega A 2, where omega is the freque ncy band.4. Now take the in verse FFT to get back to the time-doma in and scale your result.This will give you a much better estimate of displaceme nt.说到底就是频域积分要比时域积分效果更好，实际测试也发现如此。原因可能是时域积分时积

6、分一次就要去趋势，去趋势就会降低信号的能量，所以最后得到的结果常常比真实幅值要小。下面做一些测试，对一个正弦信号的二次微分做两次积分，正弦频率为 50Hz , 采样频率 1000Hz , 恢复效果如下时域积分位移速度500vwwwrmm-5000O.102Q5C&0.70.8091频域积分位移原始信号恢复信号?rnnrmYmm00.40.6C.&0708091速度500原始信号 |1immmwmmn50000.10.20 305C&0 70.80 91欢迎下载2可见恢复信号都很好（对于50Hz 是这样的效果）分析两种方法的频率特性曲线如下时域积分位移积分的颓率特性由线1.1频

7、域积分欢迎下载3位移祝弁的频率特性曲言夷1u n _ 5601400OOO6-O 20O250 30O可以看到频域积分得到信号更好，时域积分随着信号频率的升高恢复的正弦幅值会降低。对于包含两个正弦波的信号，频域积分正常恢复信号，时域积分恢复的高频信息有误差；对于有噪声的正弦信号，噪声会使积分结果产生大的趋势项（不是简单的二次趋势），如下图欢迎下载4对此可以用滤波的方法将大的趋势项去掉。测试的代码如下%测试积分对正弦信号的作用clcclearclose all% 原始正弦信号ts = 0.001;fs = 1/ts;t = 0:ts:1000*ts;f = 50;dis = sin(2*pi*f

8、*t); %位移vel = 2*pi*f*cos(2*pi*f*t); %速度acc = -(2*pi*f). A 2.*sin(2*pi*f*t); %加速度%多个正弦波的测试% fl = 400;% disl = sin(2*pi*f1*t); %位移% veil = 2*pi*f1.*cos(2*pi*f1*t); %速度% accl = -(2*pi*f1).A2.*sin(2*pi*f1*t); %加速度% dis = dis + disl;% vel = vel + vell;欢迎下载5% acc = acc + accl;%结：频域积分正常恢复信号，时域积分恢复加入的高频信息有误

9、差%加噪声测试acc = acc + (2*pi*f). A2*0.2*ra ndn (size(acc);%结：噪声会使积分结果产生大的趋势项figureax(1) = subplot(311);plot(t, dis), title( 位移 )ax(2) = subplot(312);plot(t, vel), title( 速度 )ax(3) = subplot(313);plot(t, acc), title( 加速度 )lin kaxes(ax, x);%由加速度信号积分算位移dis int, veli nt = In tFc n( acc, t, ts, 2);axes(ax(2)

10、; hold onplot(t, velint, r), legend(原始信号 ，恢复信号 )axes(ax(1); hold onplot(t, disint, r), legend(原始信号 ，恢复信号 )% 测试积分算子的频率特性n = 30;amp = zeros( n, 1);f = 5:30 40:10:480;figurefor i = 1:le ngth (f)fi = f(i);acc = -(2*pi*fi).A2.*sin(2*pi*fi*t); %加速度disi nt, veli nt = In tFcn(acc, t, ts, 2); %积分算位移amp(i) =

11、sqrt(sum(dis in t.A2)/sqrt(sum(dis.A2);plot(t, disi nt)draw now% pause endclosefigureplot(f, amp)title( 位移积分的频率特性曲线)xlabel(f)ylabel( 单位正弦波的积分位移幅值)以上代码中使用IntFcn 函数实现积分，它是封装之后的函数，可以实现时域积分和频域积分，其代码如下欢迎下载6%积分操作由加速度求位移，可选时域积分和频域积分fun cti on dis int, veli nt = In tFc n( acc, t, ts, flag)if flag = 1%时域积分di

12、sint, velint = IntFcn_Time(t, acc);vele nergy = sqrt(sum(veli nt.A 2);veli nt = detre nd(veli nt);velree nergy = sqrt(sum(veli nt.A2);veli nt = veli nt/velree nergy*vele nergy;dise nergy = sqrt(sum(dis in t.A2);dis int = detre nd(dis in t);disree nergy = sqrt(sum(dis in t.A2);disi nt = disi nt/disre

13、e nergy*dise nergy; %此操作是为了弥补去趋势时能量的损失%去除位移中的二次项p = polyfit(t, disint, 2);dis int = dis int - polyval(p, t);else%频域积分veli nt = iomega(acc, ts, 3, 2);veli nt = detre nd(veli nt);dis int = iomega(acc, ts, 3, 1);%去除位移中的二次项p = polyfit(t, disint, 2);dis int = dis int - polyval(p, t);endend其中时域积分的子函数如下%时域

14、内梯形积分fun cti on xn, vn = In tFc n_Time(t, an)vn = cumtrapz(t, an);vn = vn - repmat(mea n(v n), size( vn ,1), 1);xn = cumtrapz(t, vn);xn = xn - repmat(mea n(xn), size(x n,1), 1);end频域积分的子函数如下 ( 此代码是一个老外编的，在频域内实现积分和微分操作)fun cti on dataout = iomega(data in, dt, data in _type, dataout_type)% %欢迎下载7Aouan

15、baij QnbAN pue % 1U81U8J0UI Aouanbaij( 乙 jo 8/vod ixou)siuiod jo aqiunN euiiujaiaci % pus !(,C o 乙 1 e eq isniu adAQnoeep oj en|eA,)J0JJ8(g adAQnoeep) jiesp!(,Co 乙 1 e eq isniu adACujeiepoj en|eA,)J0JJ8(C od人厂 upiep) 上 o 乙 7eq;ie aje adAQnoeep pueadACujeiep iei|i o ns8|e|/| % % % % % % %U0lieJ8|800V

16、%-%AipopA - 乙 二 od 人厂inoeep % iu8iu8oe|ds!ci - L% U0lieJ8|800V -%AIIOOPA - 乙=od 人厂 uQiep %iu8iu8oe|ds!ci - L%( 8|diues od spuooas jo suun) luoiuaou! etui; = IP % % adAQnoeep odd jo eep iujoj8Ae/v indjo 二 noeep % od 人厂 uQiep odd jo eep iujoj8Ae/v dui = uieep %-% :S8|qeueA %LUJOj % pajisap eq; oui pap

17、aAuoo sip 屮 Qnoeep) indino iq 6ui;|nsai uQiuop %etui; eq; oui |oeq spaAuoo ueq; pue iuoj indino pajisap oui paAuoo % A|ise88J01U oi apjo ui uieiuop-Aouanbo#。屮 oui siujojsueji (uieiep) % eeplUJopABM ) o Aejje ue sa|ei duos oi|丄 saijas-saiuii uonejapooe %o gooQ 八uaiuaoeidsip eq;ie oi sauas-aiun uoneja

18、pooe %o gooQ 八 uaiuaoeidsip6uip8Auoo oj duos avilVIAl e siV93IAIOI % % % %N = 2M extpow2(max(size(data in);df = 1/(N*dt);Nyq = 1/(2*dt);% Save freque ncy arrayiomega_array = 1i*2*pi*(-Nyq : df : Nyq-df);iomega_exp = dataout_type - data in _type;% Pad datai n array with zeros (if n eeded) size1 = siz

19、e(data in ,1);size2 = size(data in, 2);if (N-size1 = 0 & N-size2 = 0)if size1 size2datai n = vertcat(datai n,zeros(N-size1,1);elsedata in = horzcat(datai n,zeros(1,N-size2);endend% Tran sform datai n in to freque ncy doma in via FFT and shift output (A)% so that zero-freque ncy amplitude is in the m