你好,
我正在努力將測量數據擬合到三角 Ksin(ax+b)+c 函數。這是我的過程的性質,我正在尋找一種自動方式來適應它的 jsl 代碼。
我寫了一個代碼來規範化數據,所以振盪在 -1 和 1 之間。
最後我想得到 a 和 b 常量(K 和 c 也是如此,但它更容易)
我的數據類型的附加示例。
BR,
約夫
我不知道這是否是您想要的,但您可以使用非線性平台解決這個問題。附件是您的數據,以及預測列(使用參數進行估計)和求解這些參數的腳本。我承認我花了很長時間才對允許解收斂的參數進行初步猜測(可能是由於我對正弦函數不熟悉——我已經很長時間沒有使用這些函數了)。但是,從概念上講,非線性平台將完成您想要的。
以下是獲得估計值的方法一種使用 FFT 函數的參數。 關於秒和小時的評論假設數據大約每 4 秒收集一次,持續大約一個小時; FFT 不關心單位。 FFT 對數據中的頻率進行分箱,使得 DC 偏移分箱之後的第一個分箱表示可以在輸入數據中完成一個完整週期的最低頻率,下一個分箱是 2 個完整週期,......中間分箱是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 之後的第 4 個 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 來遍歷您的數據並註意每次數據從負數變為正數。 這是正弦波週期的開始。 到下一個零交叉點的距離就是周期。
您可能會遇到似乎多次交叉的嘈雜數據的問題。
零交叉比峰值更容易檢測。
我用這個想法毛刺.