设为首页收藏本站联系我们

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛 › 展示 ›振动理论› 信号处理 › 查看内容

matlab数值积分的实现：时域积分和频域积分

2016-12-28 11:34| 发布者: weixin| 查看: 1545| 评论: 7|来自: 了凡春秋的博客

摘要: 积分操作主要有两种方法：时域积分和频域积分，积分中常见的问题就是会产生二次趋势。关于积分的方法，在国外一个论坛上有人提出了如下说法，供参考。
积分操作主要有两种方法：时域积分和频域积分，积分中常见的问题就是会产生二次趋势。关于积分的方法，在国外一个论坛上有人提出了如下说法，供参考。
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，恢复效果如下：

时域积分
01.png
频域积分
02.png

可见恢复信号都很好（对于50Hz是这样的效果）。分析两种方法的频率特性曲线如下：

时域积分
03.png
频域积分
04.png


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

对于包含两个正弦波的信号，频域积分正常恢复信号，时域积分恢复的高频信息有误差；对于有噪声的正弦信号，噪声会使积分结果产生大的趋势项（不是简单的二次趋势），如下图：
05.png
对此可以用滤波的方法将大的趋势项去掉。测试的代码如下：
  1. % 测试积分对正弦信号的作用
  2. clc
  3. clear
  4. close all
  5. % 原始正弦信号
  6. ts = 0.001;
  7. fs = 1/ts;
  8. t = 0:ts:1000*ts;
  9. f = 50;
  10. dis = sin(2*pi*f*t);
  11. % 位移
  12. vel = 2*pi*f.*cos(2*pi*f*t);
  13. % 速度
  14. acc = -(2*pi*f).^2.*sin(2*pi*f*t);
  15. % 加速度
  16. % 多个正弦波的测试
  17. % f1 = 400;
  18. % dis1 = sin(2*pi*f1*t);
  19. % 位移
  20. % vel1 = 2*pi*f1.*cos(2*pi*f1*t);
  21. % 速度
  22. % acc1 = -(2*pi*f1).^2.*sin(2*pi*f1*t);
  23. % 加速度
  24. % dis = dis + dis1;
  25. % vel = vel + vel1;
  26. % acc = acc + acc1;
  27. % 结：频域积分正常恢复信号，时域积分恢复加入的高频信息有误差
  28. % 加噪声测试
  29. acc = acc + (2*pi*f).^2*0.2*randn(size(acc));
  30. % 结：噪声会使积分结果产生大的趋势项
  31. figure
  32. ax(1) = subplot(311);
  33. plot(t, dis), title('位移')
  34. ax(2) = subplot(312);
  35. plot(t, vel), title('速度')
  36. ax(3) = subplot(313);
  37. plot(t, acc), title('加速度')
  38. linkaxes(ax, 'x');
  39. % 由加速度信号积分算位移
  40. [disint, velint] = IntFcn(acc, t, ts, 2);
  41. axes(ax(2));hold on
  42. plot(t, velint, 'r'),
  43. legend({'原始信号', '恢复信号'})
  44. axes(ax(1));hold on
  45. plot(t, disint, 'r'),
  46. legend({'原始信号', '恢复信号'})

  48. % 测试积分算子的频率特性
  49. n = 30;
  50. amp = zeros(n, 1);
  51. f = [5:30 40:10:480];
  52. figure
  53. for i = 1:length(f)
  54.     fi = f(i);
  55.     acc = -(2*pi*fi).^2.*sin(2*pi*fi*t);
  56.     % 加速度
  57.     [disint, velint] = IntFcn(acc, t, ts, 2);
  58.     % 积分算位移
  59.     amp(i) = sqrt(sum(disint.^2))/sqrt(sum(dis.^2));
  60.     plot(t, disint)
  61.     drawnow
  62. end
  63. close
  64. figure
  65. plot(f, amp)
  66. title('位移积分的频率特性曲线')
  67. xlabel('f')
  68. ylabel('单位正弦波的积分位移幅值')
复制代码


以上代码中使用IntFcn函数实现积分，它是封装之后的函数，可以实现时域积分和频域积分，其代码如下：
% 积分操作由加速度求位移，可选时域积分和频域积分
  1. function [disint, velint] = IntFcn(acc, t, ts, flag)
  2. if flag == 1
  3.     % 时域积分
  4.     [disint, velint] = IntFcn_Time(t, acc);
  5.     velenergy = sqrt(sum(velint.^2));
  6.     velint = detrend(velint);
  7.     velreenergy = sqrt(sum(velint.^2));
  8.     velint = velint/velreenergy*velenergy;
  9.     disenergy = sqrt(sum(disint.^2));
  10.     disint = detrend(disint);
  11.     disreenergy = sqrt(sum(disint.^2));
  12.     disint = disint/disreenergy*disenergy;
  13.     % 此操作是为了弥补去趋势时能量的损失
  14.     % 去除位移中的二次项
  15.     p = polyfit(t, disint, 2);
  16.     disint = disint - polyval(p, t);
  17. else
  18.     % 频域积分
  19.     velint =  iomega(acc, ts, 3, 2);
  20.     velint = detrend(velint);
  21.     disint =  iomega(acc, ts, 3, 1);
  22.     % 去除位移中的二次项
  23.     p = polyfit(t, disint, 2);
  24.     disint = disint - polyval(p, t);
  25. end
  26. 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


频域积分的子函数如下（此代码是一个老外编的，在频域内实现积分和微分操作）
  1. function dataout =  iomega(datain, dt, datain_type, dataout_type)
  2. %%%%%%%%%%%%%%%%
  3. %
  4. %   IOMEGA is a MATLAB script for converting displacement, velocity, or
  5. %   acceleration time-series to either displacement, velocity, or
  6. %   acceleration times-series. The script takes an array of waveform data
  7. %   (datain), transforms into the frequency-domain in order to more easily
  8. %   convert into desired output form, and then converts back into the time
  9. %   domain resulting in output (dataout) that is converted into the desired
  10. %   form.
  11. %
  12. %   Variables:
  13. %   ----------
  14. %
  15. %   datain       =   input waveform data of type datain_type
  16. %
  17. %   dataout      =   output waveform data of type dataout_type
  18. %
  19. %   dt           =   time increment (units of seconds per sample)
  20. %
  21. %                    1 - Displacement
  22. %   datain_type  =   2 - Velocity
  23. %                    3 - Acceleration
  24. %
  25. %                    1 - Displacement
  26. %   dataout_type =   2 - Velocity
  27. %                    3 - Acceleration
  28. %
  29. %
  30. %%%%%%%%%%%%%%%%
  31. %   Make sure that datain_type and dataout_type are either 1, 2 or 3
  32. if (datain_type < 1 || datain_type > 3)
  33.     error('Value for datain_type must be a 1, 2 or 3');
  34. elseif (dataout_type < 1 || dataout_type > 3)
  35.     error('Value for dataout_type must be a 1, 2 or 3');
  36. end
  37. %   Determine Number of points (next power of 2), frequency increment
  38. %   and Nyquist frequency
  39. N = 2^nextpow2(max(size(datain)));
  40. df = 1/(N*dt);
  41. Nyq = 1/(2*dt);
  42. %   Save frequency array
  43. iomega_array = 1i*2*pi*(-Nyq : df : Nyq-df);
  44. iomega_exp = dataout_type - datain_type;
  45. %   Pad datain array with zeros (if needed)
  46. size1 = size(datain,1);
  47. size2 = size(datain,2);
  48. if (N-size1 ~= 0 && N-size2 ~= 0)
  49.     if size1 > size2
  50.         datain = vertcat(datain,zeros(N-size1,1));
  51.     else
  52.         datain = horzcat(datain,zeros(1,N-size2));
  53.     end
  54. end
  55. %   Transform datain into frequency domain via FFT and shift output (A)
  56. %   so that zero-frequency amplitude is in the middle of the array
  57. %   (instead of the beginning)
  58. A = fft(datain);
  59. A = fftshift(A);
  60. %   Convert datain of type datain_type to type dataout_type
  61. for j = 1 : N
  62.     if iomega_array(j) ~= 0
  63.         A(j) = A(j) * (iomega_array(j) ^ iomega_exp);
  64.     else
  65.         A(j) = complex(0.0,0.0);
  66.     end
  67. end
  68. %   Shift new frequency-amplitude array back to MATLAB format and
  69. %   transform back into the time domain via the inverse FFT.
  70. A = ifftshift(A);
  71. datain = ifft(A);
  72. %   Remove zeros that were added to datain in order to pad to next
  73. %   biggerst power of 2 and return dataout.
  74. if size1 > size2
  75.     dataout = real(datain(1:size1,size2));
  76. else
  77.     dataout = real(datain(size1,1:size2));
  78. end
  79. return
复制代码

收藏 分享 邀请
上一篇：如何理解卷积运算的物理意义下一篇：新手入门小知识：浅析到底什么是信号与频谱
发表评论

最新评论

引用 mxlzhenzhu 2018-6-11 20:07
这个问题，我研究过。除了模拟积分还没有真正做实验，所得结论如下：
数值积分有积分常数问题，很难处理。数值微分存在高频数值噪声。总体来讲，数值微分优于数值积分。基于傅立叶变换的积分，只适用于周期信号，对非周期信号，不能得到正确结果。
模拟积分也有问题，主要是受到测试系统滤波器常数限制（放大电路回路电容限制），影响低频结果，而滤波器在高频的增益不是无穷大，所以又影响高频积分结果，也是compromised结果。

总之，在加速度，速度和位移之间，通过测试a得到v和d几乎不可能，时间越长，误差越大，通过测试d可以得到v，加速度直接测。
另外一个经验，就是测试d和a，然后kalman滤波，结果那是相当好，适用于非稳态。当然还有一个细节要提醒你，数值微分以后，低通滤波是可以的，只是普通fir滤波有相移，可以想法消除相移。移动平均或者光顺也是一种低通。
引用 mxlzhenzhu 2018-6-11 20:12
想通过测试a，傅立叶积分得到v和d的路径是危险的！难道不能加窗？
Noooo!因为加窗以后，你看时域的v和d都是两端为零的，加窗对v和d是毁灭性的。
引用 impulse 2018-6-12 16:55
mxlzhenzhu 发表于 2018-6-11 20:07
这个问题，我研究过。除了模拟积分还没有真正做实验，所得结论如下：
数值积分有积分常数问题，很难处理。 ...

“在加速度，速度和位移之间，通过测试a得到v和d几乎不可能，时间越长，误差越大”？
怎么是不可能的？我们一直在用，并且也不是你说的时间越长，误差越大，而是随着时间增加，效果越好，前1-2s因为滤波器问题，积分结果有震荡，连续1-2s后就能稳定。
引用 mxlzhenzhu 2018-6-12 18:00
请楼上的，贴出来你积分得到数据吧，我来检查一下精度。
引用 17621618700 2021-3-17 23:46
学习了！
引用 17621618700 2021-3-18 22:44
为什么% 去除位移中的二次项
    p = polyfit(t, disint, 2);
    disint = disint - polyval(p, t);

我们要去除二次项？
引用 TestGuru 2021-3-23 00:50
无论时域或频域积分，都需高通滤波，以消除量化噪声在积分后的累积。

查看全部评论(7)

相关分类

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-26 06:22 , Processed in 0.063983 second(s), 24 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部