你好,
我正在努力将测量数据拟合到三角 Ksin(ax+b)+c 函数。这是我的流程的本质,我正在寻找一种通过 jsl 代码自动适应它的方法。
我编写了一段代码来标准化数据,以便振荡在 -1 和 1 之间。
最后我想得到 a 和 b 常数(K 和 c 也是如此,但它更容易)
附上我的数据类型的示例。
BR,
约阿夫
你可以使用过零检测器,这将是一小段 JSL,用于遍历您的数据并在每次数据从负值变为正值时进行注意。 这是正弦波周期的开始。 到下一个过零的距离就是周期。
您可能会遇到多次交叉的噪声数据问题。
过零比峰值更容易检测。
我用了这个想法毛刺。
您好,感谢大家花费时间。 对于我的应用程序来说,这看起来是一个很好的方法。但是-你是如何“自动”获得经期长度的? 我没听懂..
约阿夫
太好了,这正是我所缺乏的直觉。相反,我操纵非线性平台中的参数值,直到图形看起来接近数据,但了解参数代表什么(并从数据中获得良好的初始猜测)是一种更好的方法。
我认为使用我所说的正弦函数的标准形式会有所帮助:Y = K*sin(A*(x -
) + C。
在您的数据中:
K = 值范围的 1/2
A = 2*pi / 周期
B = 从 0 向右移动到正弦波的起点
C = Y 值的平均值
可以使用分布平台和图形生成器中的这些定义找到非线性回归的起始估计:
K = 1/2 *(最大值 - 最小值)= 1/2*(47 - 17.6) = 14.7
A = 2 * 3.14 / 880 = 0.007
乙 = 662
C = 34.1


非线性r 回归使用这种方法离子收敛得非常快。该模型使用这些值在 5 次迭代中收敛。 见下文

以下是如何获得估计值A使用 FFT 函数的参数。 关于秒和小时的评论是假设大约一个小时内每 4 秒收集一次数据; FFT 不关心单位。 FFT 对数据中的频率进行分档,使得 DC 偏移分档之后的第一个分档代表可以在输入数据中形成一个完整周期的最低频率,下一个分档是 2 个完整周期,...中间的分档是 ndata/ 2个周期。 之后的垃圾箱是第一个垃圾箱的镜像。
dt = Open( "z:\example data.jmp" );
// by hand, I made a formula column to look at x-lag(x) and
// decided the first 79 and last 8 were not representative...
rowstart = 80;
rowstop = 937;
timeStart = dt:x[rowstart];
timeStop = dt:x[rowstop];
ydata = dt[rowstart :: rowstop, {y}]; // throw away a few with irregular time
{r, i} = FFT( {ydata} );
nBins = N Items( r );
sampleRate = (timeStop - timeStart) / (nBins - 1); // imagine there were only 2 bins...thus -1.
// the samplerate is ~4 seconds/sample...
// the highest frequency this can resolve is 1 cycle/ 8 seconds. (.125 Hz)
// frequencies higher than that will cause aliasing unless they are filtered out before this.
// the bin [nBins/2] represents the maximum frequency that can be resolved
lastBinFrequency = 1 / (2 * sampleRate); // ~ 0.125 Hz
// the first bin [1] is the DC, or offset from 0, component. We'll ignore it.
// the next bin [2] represents the frequency that has one cycle over the length of the data.
// the lowest non-DC frequency is one cycle per hour (the length of time the data was collected).
binWidth = 1 / (timeStop - timeStart); // ~0.0003Hz, or ~1 cycle per hour
// note: lastBinFrequency == (nbins-1)*binWidth/2 -- each bin is this wide, in frequency
// the remaining bins are a mirror of the first and ignored.
// combine real and imaginary bin components to get an amplitude
p = Sqrt( r ^ 2 + i ^ 2 );
dtGraph2 = As Table( ((rowstart :: rowstop) - rowstart)` || p );
dtGraph2 << Graph Builder(
Size( 1300, 500 ),
Show Control Panel( 0 ),
Show Legend( 0 ),
Variables( X( :Col1 ), Y( :Col2 ) ),
Elements( Points( X, Y, Legend( 5 ) ) ),
SendToReport(
Dispatch( {}, "Col1", ScaleBox, {Format( "Best", 12 ), Min( 0 ), Max( nBins / 2 ), Inc( 100 ), Minor Ticks( 4 )} ),
Dispatch( {}, "Col2", ScaleBox, {Format( "Best", 12 ), Min( 0 ), Max( 30000 ), Inc( 5000 ), Minor Ticks( 1 )} ),
Dispatch( {}, "graph title", TextEditBox, {Set Text( "FFT Bins" )} ),
Dispatch( {}, "X title", TextEditBox, {Set Text( "bin number" )} ),
Dispatch( {}, "Y title", TextEditBox, {Set Text( "size" )} )
)
);
// there is a stong DC component in this data (all values > 0) so bin 1 is the max.
If( Loc Max( p ) != 1,
Throw( "something odd" )
);
p[1] = 0; // clear the DC value
// next max after DC:
maxBinIndex = Loc Max( p ); // 5. the data repeats about (5-1)==4 times in this frame.
Write(
"\!nBiggest frequency component at ",
binWidth * (maxBinIndex - 1),
" Hz or ",
Round( 1 / (binWidth * (maxBinIndex - 1)) ),
" seconds per cycle"
);
// where does a==.007 come from?
secondsPerCycle = 1 / (binWidth * (maxBinIndex - 1));
totalSeconds = timeStop-timeStart;
cyclesInSample = totalSeconds/secondsPerCycle;//about 4
timerange = timeStop-timeStart;
a_from_fft = cyclesInSample*2*pi()/timerange; // .0073...
show(a_from_fft);
One strong response in 4th bin after DC binOne strong response in 4th bin after DC binDC bin 之后的第 4 个 bin 中出现强烈响应
最大频率分量为 0.00116448326055313 Hz 或每个周期 859 秒
a_from_fft = 0.007316664113164;
DC 分量与数据的平均值有关, C价值。
我不知道这是否是您想要的,但您可以使用非线性平台来解决这个问题。附件是您的数据,以及预测器列(使用参数进行估计)和求解这些参数的脚本。我承认,我花了相当长的时间才对允许解决方案收敛的参数进行了初步猜测(可能是由于我对正弦函数缺乏熟悉 - 我已经很久没有使用此类函数了)。但是,从概念上讲,非线性平台将实现您想要的。