非平稳数据平稳化的方法:去直流分量和消除趋势项

a1runner

若观测序列具有增长、下降或其他的变化趋势，或具有周期变化，都属于非平稳序列。
平稳化序列的方法较多。
对于趋势项，可以用计算机拟合相应的函数。事先不假定趋势项是什么函数形式，而给出一个一般多项式的函数形式，然后用逐步回归方法在计算机上加以筛选来确定趋势项的数学表达式。也可以进行差分处理，对于具有线性趋势的序列进行一次差分即可，若属于二次多项式，则需要进行二次差分，等等。对于具有周期变化的序列，可以通过傅立叶谐波分析，提取周期项，或者进行季节差分处理，使之平稳化。
对于要研究的多变量序列，可以分别进行上述处理，将趋势项和周期项处理或提取后，变为平稳序列，进行多维自回归分析。

将时间序列平稳化，而后运用平稳过程的理论和分析手段来建模，分析，预报，从而实现对原时间序列的分析预报。

观察数据序列看是否符合典型分解式的模型：=趋势项m+季节项s+平稳序列Y

趋势项--是数据序列所体现的缓慢变化的趋势

季节项-数据体现的周期性变化

平稳化需要去除其中的趋势项，季节项：

首先对于没有季节项的时间序列：趋势项+平稳序列 可以通过最小二乘法拟合趋势项函数，对原时间序列减去趋势项得到的平稳序列进行分析；或者通过滑动平均来估计趋势项，利用滑动窗口移动/指数平均均来实现对数据序列的平滑，得到趋势项；或者通过差分方法去除趋势项，循环执行差分（1次或2次），直到序列满足一个平稳序列模型。注：原时间序列值增加，相应差分序列波动变化增大，可以通过ln对数化，来减小原时间序列值每步增加的幅度，再作分析。

对于有季节项+趋势项，考虑季节项的删除，设季节项周期为d,趋势项按周期d分段平均-而后原时间序列减去趋势项，将各个周期内对应时间序列值求平均，便得到了相应周期内季节项估计；或者采用滑动窗口平均得到趋势项，原序列减去趋势项，各个周期内对应时间序列wk(k=1,...d)求均值，为了满足周期内季节项之和为0，再减去周期内时刻均值之和均值(w1+..Wd)/d。或者通过延迟d步差分的方法，Ad=xt-x(t-d)=m(t)-m(t-d)+Y(t)-Y(t-d);因为季节项周期为d，所以差为0，然后变成了情况1种的典型分解式趋势项m(t)-m(t-d)，噪声项Y(t)-Y(t-d)，利用前面的方法加以去除趋势项即可分析。

a1runner

在振动测试中采集到的振动信号数据，由于放大器随温度变化产生的零点漂移、传感器频率范围外低频性能的不稳定以及传感器周围的环境干扰等，往往会偏离基线，甚至偏离基线的大小还会随时间变化。
偏离基线随时间变化的整个过程被称为信号的趋势项。
直流分量可理解成一根水平线，而趋势项一般是指信号整体的变换趋势，可以是线性变换，也可是多项式曲线，或者是其他曲线都可以.两者都是慢变成分,一般认为是仪器漂移所导致,所以要去掉.
直流分量就是认为这个漂移不随时间发生变化,因此去零或去直流. 而慢变的漂移就是所谓趋势项,一般通过曲线拟合(短数据),或高通滤波来消除(长数据). 当然,消除趋势项的算法肯定能去直流.

a1runner

基线调整－最小二乘法消除多项式趋势项
王济的《MATLAB在振动信号处理中的应用》一书的原程序：
fni=input('频域积分——输入数据文件名：','s');
fid=fopen(fni,'r');
fno=fscanf(fid,'%s',1);
x=fscanf(fid,'%f',[1,inf]);
status=fclose(fid);

sf=1000;    %采样频率
fmin=0.1;     %最小截止频率
fmax=200;       %最大截止频率
c=9810.5;

n=length(x);
t=0:1/sf:(n-1)/sf;
nfft=2^nextpow2(n);
y=fft(x,nfft);
df=sf/nfft;
ni=round(fmin/df+1);
na=round(fmax/df+1);
dw=2*pi*df;
w1=0:dw:2*pi*(0.5*sf-df);
w2=-2*pi*(0.5*sf-df):dw:-dw;
w=[w1,w2];
a=zeros(1,nfft);
a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1);
real(y)=imag(a);
imag(y)=-real(a);
a=zeros(1,nfft);
a(ni:na)=y(ni:na);
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
y=ifft(a,nfft);
y=real(y(1:n))*c;
subplot(2,1,1);
plot(t,x);
xlabel('时间(s)');
ylabel('加速度(g)');
grid on;
subplot(2,1,2);
plot(t,y);
xlabel('时间(s)');
ylabel('速度(m/s)');
grid on;
fid=fopen('a.txt','w+');
for k=1:n
    fprintf(fid,'%f\n',y(k));
end
status=fclose(fid);

a1runner

数值积分结果误差大，究其原因，主要有两个方面。
第一，积分过程中受到“平移项”（直流偏置）的影响。可以通过简单的减去积分结果的均值来消除平移项的影响。国内外绝大部分测试软件中的数值积分都有此项功能。
第二，积分过程中受到“趋势项”的影响。
所谓趋势项是在测试信号中存在线性项或缓变的非线性项成分。工程实际测量的信号大部分是复杂周期信号与随机信号的混合，而且周期信号往往是研究对象。趋势项的存在会使数值积分的结果产生很大的误差，严重地背离真实情况。因此，在测试信号（积分求速度、位移时）中常要先消除趋势项，这是积分中的一个重要的中间步骤。
产生趋势项的原因主要有以下几种：
1）采样时未对原始信号进行适当处理，如高通滤波，使得信号中含有周期大于采样时间的极低频信号。
2）由于外界原因，包括传感器或仪器的零点漂移、基础运动等引起的信号波形偏移。
3）在截取记录时，样本长度选择不当。
4）由于操作不当，信号经过积分放大后产生趋势项，如零点未调准产生的常数，经积分后成为一条斜直线，低频噪声经积分放大后成为缓慢变化的趋势项。
可见产生趋势项的原因是很复杂的，要在实际采集的数据中完全避免趋势项十分困难。
去趋势项:
1）通过高通滤波去趋势项
测量信号中的趋势项和实际信号的周期比较要长很多，此时可以先把测到的信号进行高通滤波处理，滤除低频趋势项。此种消除趋势项的方法要求实验人员必须对所测结构的基频有一个大致的了解（可以通过有限元计算或先对实测信号进行频谱分析确定），在对结构基频了解的基础上，设置一个合适的高通滤波器的截止频率（大量测试数据表明截止频率设为基频的60%左右时消除趋势项的影响效果较好）。让测得的数据先经去直流和高通滤波后再进行数值积分，可以取得较为满意的积分效果。
2）EMD分解法去趋势项
根据定义，EMD分解方法假定任何待分解信号都由一组固有振动模式构成，并据此将信号分解为若干固有模式函数（intrinsic mode function，IMF）的和与趋势项的叠加。
具体实施时，可在EMD分解中指定每个IMF分量的时延或频率（时延的倒数），以保证在分解得到的IMF中不包含有低于此设定频率值的信号（积分要消除的趋势项）。例如，在第i个IMF分量Ci的筛选过程中，设定“间歇检测准则”对应的频率fc，则完成筛选过程后，分量Ci中不含有频率低于fc的成分。可见，这里的频率fc类似于高通滤波器中设定的截止频率值。具体做法是是，第一步进行数值积分消除直流的影响，第二步将信号进行EMD分解，提取出合适的趋势项，用各个IMF分量重构得到消除了趋势项的数值积分结果。
——张永强、宋建江、屠良尧、薛姗，软件数值积分误差原因分析及改进办法，机械强度，第28卷第3期，p419~423
本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/ap0108220/archive/2010/03/16/5385211.aspx

a1runner

如果趋势项的波动频率与信号频率偏离比较远，则可以用：
1）直接对数据差分
2）高通滤波
如果数据比较短，则最好用最小二乘法拟合趋势项，然后从原始数据中减去这个趋势项

afnuaa

请教，关于去直流去趋势项，有没有比较好的书对这方面有系统介绍，谢谢

太阳vib

回复 5 # a1runner 的帖子

楼主能不能上传个高通滤波去趋势项的程序啊，这样有个实例对概念理解也是有帮助的

liuxianming456

此贴不回天诛地灭！
太深刻了，以前只知道去除直流分量

倔强的笨蛋

a1runner 发表于 2010-12-20 15:02

                               
登录/注册后可看大图

如果趋势项的波动频率与信号频率偏离比较远，则可以用：
1）直接对数据差分
2）高通滤波

想问下楼主，楼主说的数据较短，是大约多少个点啊？

lyw234

xuexi,学习，谢谢楼主。

sumier

太阳vib 发表于 2012-7-24 17:23

                               
登录/注册后可看大图

回复 5 # a1runner 的帖子

楼主能不能上传个高通滤波去趋势项的程序啊，这样有个实例对概念理解也是有帮助 ...

请问你没有找到高通滤波去趋势项的程序呀，我最近正在搞这个，不会做啊，求赐教