赛特居士分享 http://blog.sciencenet.cn/u/SciteJushi

博文

[转载]滑动窗口极点模型拟合——混沌序列的低阶逼近与谱估计

已有 723 次阅读 2018-7-9 12:24 |系统分类:科研笔记| 复指数信号, 谱估计, 周期图, 滑动窗口, 混沌 |文章来源:转载

原载http://blog.sina.com.cn/s/blog_729a92140102xkcq.html

               居士在博客中的科技研讨,基于了重要实用背景,根源于几所大学里亲身研发中的需要和问题。由2011年3月和4月贴出的那些自相关序列的曲线图,可看出,可用极点模型去拟合随机过程或混沌序列的单侧自相关序列。这就又到了一个重要研究领域——现代谱估计,也类似于磁共振谱数据分析中拟合FID信号的情形,可以只考虑严格地衰减的分量的模型,在那里,人们可将K空间数据预先与适当的振荡频率为零的指数衰减的实数序列相乘,既强化了衰减性,还能抑制末尾段的噪声的影响(类似于图像重建前在K空间加窗,还抑制成像系统状态的快速切换不完美的窜扰等,美化了图像)。
              不过,本短文考察:经典的非参数化谱估计,即周期图;在获得谱之前,用阶数较低的滑动窗口极点模型,直接逼近混沌吸引子数据的原始序列。
              演示试验程序ExpFitChaos1.sce,如图片1.的右半部所示。程序主体有2层for循环(行号10至31)。内层循环变量,指定数据段xx的序号,对应了滑动窗口的位置。每个数据段的长度,即窗口的宽度,等于窗口的滑动步长的2倍。
             外层循环变量,指定模型“阶数”的取值,约为:窗口宽度值,除以32、64、128、0.1、-1。当阶数小于2时,直接调用Scilab-6.0.0_x64自带的信号处理模块的谱估计函数(居士以其为参照标准,非常感谢作者和免费的Scilab)pspect,估计当前窗口内的数据的谱;否则,利用FFT,简便地计算谱序列,即:
去均值(行号23);范数归一化(行号24);FFT(行号25);取FFT输出序列的模值的平方(行号26),作为当前数据段的谱估计结果。
              这种简便谱序列,表达了功率分布的相对状况,不计绝对的量值。实际上,对于DFT或IDFT,所谓的“时间”、“频率”、“功率”以及它们的“单位”等,也只是虚拟的“名”。可在特定实用技术链中,再赋予切实物理意义。对于均值为零的实数序列的谱,只截取其中半段长。
             这里,在使用pspect时,每次也仅输入了单个窗口内的数据,尽管,形式上仍按其要求,输入了另外两个量,即,移位步长和单个数据段的长度。在通常应用中,它允许一次输入很多段长的序列。Scilab帮助文档,复制于本文后,它描述其输出谱并强调了归一化频率区间[0,1]和函数linspace,反而不像在正常的傅立叶变换。实际上,宜用区间[0,1),视“1”与“0”处相同;联系到模数转换,考虑边界“-0.5”与“0.5”,面对采样频率边界模糊(《小波无罪Shannon采样定理也没有违法》,2015-11-22)问题。
             当指定的阶数不小于2但小于滑动步长的88%时,在去均值之前,先得到低阶极点模型逼近(不过,如果,原均值,相对地太强,影响到了其余分量的数字表达和提取,那么,此前,宜先预处理一次均值)的结果,并以之取代原段数据用于简便地计算谱序列。
             在内、外层循环结束语句之间,统计所有数据段的估计结果,得其:均值序列,标准差序列与均值加上标准差之和的序列的比值。在外层循环结束后,画出统计结果的曲线图,每个阶参数值,对应一条曲线。另外,画出包含前20个数据段的3D吸引子。
              图片1.的左上部的图形窗,显示了用Lorenz吸引子(《Lorenz混沌吸引子以及初始条件敏感性和不敏感性》,2011-03-24)数据时的一组结果。取用了100个数据窗,滑动步长为512点。利用32位平台的旧bin文件组,重新在64位系统中32位平台Scilab-5.3.3上,预先计算生成了吸引子数据文件。32阶逼近时的曲线(黑色),与pspect(红色)以及不做低阶逼近时(青色)的形状,很相似。
              看到中间那幅曲线图时,可以意识到:谱估计的标准差曲线,与均值曲线,太相近了。居士总希望,能减小谱估计的方差。方差越大而平均次数又越少的场合,就可能越难精准地重复和比较,相应估计的平均值曲线。这是个痼疾。
              孤立序列的谱估计,有限地反映了单个长序列的某些统计特性。由低阶极点模型参数,获得了适当的谱估计,这不意味着表达足够好了轨道坐标系中相互间非线性地耦合着的所有分量,不能排除轨线的扭曲变形。
              新浪赛特居士SciteJushi-2018-07-09。

图片1.滑动窗口极点模型逼近混沌序列的试验程序和结果


附1. Scilab-6.0.0_x64中pspect的帮助文档
Scilab Help >> Signal Processing > Spectral estimation > pspect 
pspect
two sided cross-spectral estimate between 2 discrete time signals using the Welch's average periodogram method.
Syntax
[sm [,cwp]]=pspect(sec_step,sec_leng,wtype,x [,y] [,wpar])
[sm [,cwp]]=pspect(sec_step,sec_leng,wtype,nx [,ny] [,wpar])
Arguments
x
vector, the time-domain samples of the first signal.
y
vector, the time-domain samples of the second signal. If y is omitted it is supposed to be equal to x (auto-correlation). If it is present, it must have the same numer of element than x.
nx
a scalar : the number of samples in the x signal. In this case the segments of the x signal are loaded by a user defined function named getx (see below).
ny
a scalar : the number of samples in the y signal. In this case the segments of the y signal are loaded by a user defined function named gety (see below). If present ny must be equal to nx.
sec_step
offset of each data window. The overlap D is given by sec_leng -sec_step. if sec_step==sec_leng/2 50% overlap is made. The overlap
sec_leng
Number of points of the window.
wtype
The window type
're': rectangular
'tr': triangular
'hm': Hamming
'hn': Hann
'kr': Kaiser,in this case the wpar argument must be given
'ch': Chebyshev, in this case the wpar argument must be given
wpar
optional parameters for Kaiser and Chebyshev windows:
'kr': wpar must be a strictly positive number
'ch': wpar must be a 2 element vector [main_lobe_width,side_lobe_height]with 0<main_lobe_width<.5, and side_lobe_height>0
sm
Two sided power spectral estimate in the interval [0,1] of the normalized frequencies. It is a row array with sec_len elements . The array is real in case of auto-correlation and complex in case of cross-correlation.
The associated normalized frequencies array is linspace(0,1,sec_len).
cwp
unspecified Chebyshev window parameter in case of Chebyshev windowing, or an empty matrix.
Description
Computes the cross-spectrum estimate of two signals x and y if both are given and the auto-spectral estimate of x otherwise. Spectral estimate obtained using the modified periodogram method.
The cross spectrum of two signal x and y is defined to be
 
The modified periodogram method of spectral estimation repeatedly calculates the periodogram of windowed sub-sections of the data contained in x and y . These periodograms are then averaged together and normalized by an appropriate constant to obtain the final spectral estimate. It is the averaging process which reduces the variance in the estimate.
For batch processing, the x and y data may be read segment by segment using the getxand gety user defined functions. These functions have the following syntax:
xk=getx(ns,offset) and yk=gety(ns,offset) where ns is the segment size and offset is the index of the first element of the segment in the full signal.
Reference
Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing, Upper Saddle River, NJ: Prentice-Hall, 1999
Examples
rand('normal');rand('seed',0);
x=rand(1:1024-33+1);


//make low-pass filter with eqfir
nf=33;bedge=[0 .1;.125 .5];des=[1 0];wate=[1 1];
h=eqfir(nf,bedge,des,wate);


//filter white data to obtain colored data
h1=[h 0*ones(1:max(size(x))-1)];
x1=[x 0*ones(1:max(size(h))-1)];
hf=fft(h1,-1); xf=fft(x1,-1);y=real(fft(hf.*xf,1));


//plot magnitude of filter
h2=[h 0*ones(1:968)];hf2=fft(h2,-1);hf2=real(hf2.*conj(hf2));
hsize=max(size(hf2));fr=(1:hsize)/hsize;plot(fr,log(hf2));


//pspect example
sm=pspect(100,200,'tr',y);smsize=max(size(sm));fr=(1:smsize)/smsize;
plot(fr,log(sm));
rand('unif');
 
 


See also
cspect — two sided cross-spectral estimate between 2 discrete time signals using the correlation method
pspect — two sided cross-spectral estimate between 2 discrete time signals using the Welch's average periodogram method.
mese — maximum entropy spectral estimation
window — compute symmetric window of various type



附2. Scilab-6.0.0_x64中linspace的帮助文档


Scilab Help >> Elementary Functions > Matrix generation > linspace 
linspace
generates linearly spaced numbers between 2 reached bounds
Syntax
row = linspace(x1, x2)
row = linspace(x1, x2, n)
Matrix = linspace(Col1, Col2)
Matrix = linspace(Col1, Col2, n)
Arguments
x1, x2
Real or complex scalars: Bounds between which values must be generated.
Col1, Col2
Column vectors of same heights of real or complex numbers.
n
integer number of requested values or columns. Default value: 100
row
row vector of n numbers.
Matrix
Matrix with n columns of numbers.
Description
linspace(x1, x2) generates a row vector of n equally spaced values ranging exactly from x1 to x2.
 
The syntax y1:y2 or y1:step:y2 like 1:0.1:%pi does the same but fixes the starting bound y1 and the step. The y2 is used as stopping bound to not be overstepped. The last value actually generated may not reach it. y2 is then not included in the result.
Instead of fixing the step to a given value, linspace fixes the final bound x2 to be exactly reached, and computes the step accordingly.
If x1 or x2 are complex numbers, then linspace(x1,x2) interpolates separately the real and the imaginary parts of x1 and x2.
If some column vectors Col1 and Col2 are provided, linspace works in a row-wise way: the resulting Matrix has the same number of rows, and n columns. We get Matrix(i,:) = linspace(Col1(i), Col2(i), n).
Examples
linspace(1, %pi, 0)         // n = 0
linspace(1, 2, 10)          // x2 > x1 : increasing values
linspace(2, 1, 10)          // x2 < x1 : decreasing values
linspace(1+%i, 2-2*%i, 5)      // with complex numbers
linspace([1:4]', [5:8]', 10)   // with input columns
 
 


--> linspace(1, %pi, 0)    // n = 0
 ans  =
    []


--> linspace(1, 2, 10)    // x2 > x1 : increasing values
 ans  =
   1.   1.111   1.222   1.333   1.444   1.556   1.667   1.778   1.889   2.



--> linspace(2, 1, 10)    // x2 < x1 : decreasing values
 ans  =
   2.   1.889   1.778   1.667   1.556   1.444   1.333   1.222   1.111   1.



--> linspace(1+%i, 2-2*%i, 5)      // with complex numbers
 ans  =
   1. +i     1.25 +0.25i   1.5 -0.5i   1.75 -1.25i   2. -2.i



--> linspace([1:4]', [5:8]', 10)   // with input columns
 ans  =
   1.   1.444   1.889   2.333   2.778   3.222   3.667   4.111   4.556   5.
   2.   2.444   2.889   3.333   3.778   4.222   4.667   5.111   5.556   6.
   3.   3.444   3.889   4.333   4.778   5.222   5.667   6.111   6.556   7.
   4.   4.444   4.889   5.333   5.778   6.222   6.667   7.111   7.556   8.
See also
colon — (:) colon operator
logspace — logarithmically spaced vector
grand — Random numbers
History
Version
Description
5.4.0
Column vectors can now be provided.
The third input argument (n) must be an integer value.
6.0
linspace(a, b, n<=0) now returns [] instead of b.
bounds are now checked against %inf or %nan values.
 
Report an issue
<< eye
Matrix generation
logspace >>




http://blog.sciencenet.cn/blog-530150-1123054.html

上一篇:[转载]滑动窗口极点模型拟合——用几种宽度的窗口试降噪
下一篇:[转载]滑动窗口极点模型拟合——混沌序列的临界重建

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-11-20 04:48

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部