你好,
我正在努力将测量数据拟合到三角 Ksin(ax+b)+c 函数。这是我的过程的性质,我正在寻找一种自动方式来适应它的 jsl 代码。
我写了一个代码来规范化数据,所以振荡在 -1 和 1 之间。
最后我想得到 a 和 b 常量(K 和 c 也是如此,但它更容易)
我的数据类型的附加示例。
BR,
约夫
我不知道这是否是您想要的,但您可以使用非线性平台解决这个问题。附件是您的数据,以及预测列(使用参数进行估计)和求解这些参数的脚本。我承认我花了很长时间才对允许解收敛的参数进行初步猜测(可能是由于我对正弦函数不熟悉——我已经很长时间没有使用这些函数了)。但是,从概念上讲,非线性平台将完成您想要的。
以下是获得估计值的方法一种使用 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);
DC bin 之后的第四个 bin 中的一个强烈响应
最大频率分量为 0.00116448326055313 Hz 或每周期 859 秒
a_from_fft = 0.0073166664113164;
DC 分量与数据的平均值有关, C价值。
我认为使用我所说的正弦函数的标准形式会有所帮助: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


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

太好了,这正是我所缺乏的直觉。相反,我在非线性平台中操纵参数值,直到图形看起来接近数据,但知道参数代表什么(并从数据中获得良好的初始猜测)是一种更好的方法。
您好,感谢您花费的时间。 对于我的应用程序来说,这看起来是一个很好的方法。但是 - 你是如何“自动”获得周期长度的? 我没有。。
约夫
你可以使用一个过零检测器,这将是一小段 JSL 来遍历您的数据并注意每次数据从负数变为正数。 这是正弦波周期的开始。 到下一个零交叉点的距离就是周期。
您可能会遇到似乎多次交叉的嘈杂数据的问题。
零交叉比峰值更容易检测。
我用这个想法毛刺.