你好,
我正在努力將測量數據擬合到三角 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價值。
我不知道這是否是您想要的,但您可以使用非線性平台來解決這個問題。附件是您的數據,以及預測器列(使用參數進行估計)和求解這些參數的腳本。我承認,我花了相當長的時間才對允許解決方案收斂的參數進行了初步猜測(可能是由於我對正弦函數缺乏熟悉 - 我已經很久沒有使用此類函數了)。但是,從概念上講,非線性平台將實現您想要的。