Akima演算法

wuyu123456發表於2024-05-02

測量資料的內插已有各種方法,如線性內插、多項式內插、樣條函式插值等,但這裡的Akima插值法具有獨特的優點。

線性內插只顧及其附近兩點的影響。

多項式內插時,低階多項式由於引數較少,內插精度很低,而使用高階多項式又會使解不穩定,出現“龍格”現象,即內插函式在插值點與實際資料符合得很好,而在插值點外出現較大的偏差。因此人們又在多項式的基礎上發展了分片多項式,即樣條函式

樣條函式既保持了多項式運算簡單的特點,又避免了多項式階數較高時數值不穩定的缺點,因而得到了廣泛的應用。但在樣條函式插值中,確定任何一個小區間上的多項式,都要考慮所有資料點對它的影響。這不僅擴大了誤差傳播的範圍,還增加了不少工作量。有時只用內插點附近的幾個資料點作為控制點來內插。

Aikma插值法和三次樣條函式一樣考慮了要素導數值的效應,因而得到的整個插值曲線是光滑的。三次樣條函式插值法具有最小模、最佳最優逼近和收斂的特性,而Aikma插值法所得曲線比樣條函式插值曲線更光順,更自然。

https://blog.csdn.net/bluels01/article/details/43561131

Akima演算法是一個插值演算法,即對於一個已知(Xi,yi)的資料集,為了讓曲線看起來平滑、自然,依據現有的資料點,透過插值,多出一些資料點的過程。

C++中的Akima演算法

#include <iostream>
#include <vector>
#include <cmath>

// Akima插值函式
double akimaInterpolation(double x, const std::vector<double>& xData, const std::vector<double>& yData) {
int n = xData.size();

int index = 0;

// Find the interval index
for (int i = 0; i < n - 1; ++i) {
if (x >= xData[i] && x <= xData[i + 1]) {
index = i;
break;
}
}

// 計算斜率
std::vector<double> slopes(n - 1); //初始化n-1個預設值為0的元素
for (int i = 0; i < n - 1; ++i) {
double dx = xData[i + 1] - xData[i];
double dy = yData[i + 1] - yData[i];
slopes[i] = dy / dx; //計算每段之間的斜率
}

// 計算權重
std::vector<double> weights(n - 1); //初始化n-1個預設值為0的元素
for (int i = 2; i < n - 2; ++i) {
weights[i] = std::abs(slopes[i + 1] - slopes[i - 1]); //計算這些權重的目的是確定每個間隔附近的斜坡的“強度”。這些權重隨後用於插值公式中,以確保插值曲線的平滑性和連續性。
}

// 計算插值
double dx = xData[index + 1] - xData[index];
double t = (x - xData[index]) / dx; //引數 t 表示區間內的歸一化位置,取值範圍為 0 到 1

//m0、m1、p0和p1是 A​​kima 插值公式中用於計算插值的係數
double m0 = slopes[index] * dx; //詳見程式碼解析1
double m1 = slopes[index + 1] * dx;
double p0 = (3 * weights[index] - 2 * m0 - m1) / dx; //這裡的3,2係數是怎麼來的詳見程式碼解析2
double p1 = (3 * weights[index + 1] - m0 - 2 * m1) / dx;
//interpolatedValue 這個公式用於計算最終插值結果,詳見程式碼解析3
double interpolatedValue =
yData[index] * (1 - t) * (1 - t) * (1 + 2 * t) +
yData[index + 1] * t * t * (3 - 2 * t) +
p0 * t * (1 - t) * (1 - t) +
p1 * t * t * (t - 1);

return interpolatedValue;
}

int main() {
std::vector<double> xData = {1, 2, 3, 4, 5};
std::vector<double> yData = {1, 3, 2, 5, 4};

// 假設輸入一個x=2.5,y輸出多少?
double interpolatedValue = akimaInterpolation(2.5, xData, yData);
std::cout << "Interpolated value at x = 2.5: " << interpolatedValue << std::endl;

return 0;
}

解釋一下上面的斜率和權重,斜率是透過相鄰點之間 k=dy/dx 來計算。而權重是區間附近斜率對這個區間影響的權重,將點i的左側斜率slopes[i - 1]和右側斜率slopes[i + 1]相減得到,存在weights[i]裡。權重隨後用於插值公式中,以確保插值曲線的平滑性和連續性。

這裡展開講一下:
在 Akima 插值中,插值曲線是透過將分段三次曲線擬合到連續資料點之間的每個區間來構建的。這些三次曲線的斜率在確定插值曲線的形狀和行為方面起著至關重要的作用。目標是確保曲線連續並遵循資料的總體趨勢,同時避免過度振盪。

透過計算權重演算法考慮了相鄰區間之間斜率的變化。權重透過捕獲資料的區域性行為並影響插值過程中每個斜率的“強度”。較大的權重表示區域斜率變化明顯,而較小的權重表示區域較平滑

在執行插值時,將權重合併到插值公式中以調整相鄰斜率的貢獻。權重充當控制不同區間斜率之間平衡的係數。此調整有助於平滑插值曲線並減少由異常值或噪聲資料點引起的突然變化。

程式碼解析1
m0、m1、p0和p1是 A​​kima 插值公式中用於計算插值的係數:

m0:此變數表示左相鄰區間的調整斜率。將當前區間 (slopes[index]) 的斜率乘以區間寬度 ( dx)得到。

m1:此變數表示右相鄰區間的調整斜率。將下一個區間 (slopes[index + 1]) 的斜率乘以區間的寬度 (dx)得到。

p0:此變數表示左相鄰區間的調整權重。它是使用當前區間 ( weights[index]) 的權重、左側區間的調整斜率 ( m0) 和右側區間的調整斜率 (m1) 計算得到。由公式(3 * weights[index] - 2 * m0 - m1) / dx確定左相鄰區間對插值的貢獻。

p1:此變數表示右相鄰區間的調整權重。它是使用下一個區間的權重 ( weights[index + 1])、左側區間的調整斜率 (m0) 和右側區間的調整斜率 (m1) 計算得到。由公式(3 * weights[index + 1] - m0 - 2 * m1) / dx確定右相鄰區間對插值的貢獻。

程式碼解析2
公式(3 * weights[index] - 2 * m0 - m1) / dx 和 (3 * weights[index + 1] - m0 - 2 * m1) / dx 是基於Akima插值方案推匯出來的。
為了理解推導,讓我們考慮 Akima 插值方案的一般形式:

y(x) = p0(x) * y0 + p1(x) * y1 + q0(x) * m0 + q1(x) * m1
在此等式中,y(x)表示特定座標處的插值x。y0和y1是x兩側的資料點, m0和m1是與資料點關聯的斜率。項p0(x)和p1(x)是資料點的權重係數,q0(x)和q1(x)是資料點關聯的斜率的權重係數。
為了確定p0(x)和p1(x),Akima 擬合使用三次多項式來確保平滑性和連續性。這些權重係數由斜率的區域性行為決定。

透過考慮Akima插值方案,我們可以推匯出程式碼中使用的具體權重公式:

對於p0(x):權重函式p0(x)決定了左鄰域的貢獻。在程式碼中,(3 * weights[index] - 2 * m0 - m1) / dx代表p0(x).
選擇特定係數3、-2和1是為了平衡斜率的影響並確保間隔邊界處的連續性。這些係數是透過數學分析和最佳化確定的。

對於p1(x):權重函式p1(x)決定了右鄰區間的貢獻。在程式碼中,(3 * weights[index + 1] - m0 - 2 * m1) / dx代表p1(x).同樣,選擇係數3、-1和-2以實現插值曲線的連續性和平滑性。

匯出這些公式中的特定係數是為了最大限度地減少插值誤差並保持曲線的連續性。它們是透過數學分析和最佳化技術確定的,以確保生成的曲線與基礎資料點緊密匹配。

程式碼解析3
  yData[index] * (1 - t) * (1 - t) * (1 + 2 * t):這一項代表左邊資料點(yData [index]) 對插值的貢獻。它乘以三次多項式“(1 - t) * (1 - t) * (1 + 2 * t)”,該多項式取決於引數“t”,範圍從 0 到 1。多項式旨在確保左側資料點的平滑過渡和適當加權。
  yData[index + 1] * t * t * (3 - 2 * t):此項表示右側資料點 (yData[index + 1]) 對插值的貢獻。它乘以三次多項式“t * t * (3 - 2 * t)”。與上面類似,這個多項式確保了右側資料點的平滑過渡和適當加權。
p0 * t * (1 - t) * (1 - t):此項表示左側相鄰區間的調整權重 (p0) 對插值的貢獻。它乘以三次多項式“t * (1 - t) * (1 - t)”。該多項式表示左側相鄰區間對插值的影響。
p1 * t * t * (t - 1):此項表示右相鄰區間的調整權重 (p1) 對插值的貢獻。它乘以三次多項式“t * t * (t - 1)”。該多項式表示右側鄰區間對插值的影響。
該方程結合了相鄰資料點的貢獻及其相應的權重來計算最終的插值。引數 t 表示區間內的歸一化位置,取值範圍為 0 到 1。它決定了相鄰資料點及其對應區間的相對權重。應用於資料點和權重的三次多項式確保插值曲線的平滑性和連續性。
把上面這些影響因素加一起就是插值點的函式值interpolatedValue了。

備註:以上資料主要來源於以下連結

https://blog.csdn.net/silent_dusbin/article/details/131316458?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-1-131316458-blog-125541966.235%5Ev43%5Epc_blog_bottom_relevance_base3&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-1-131316458-blog-125541966.235%5Ev43%5Epc_blog_bottom_relevance_base3&utm_relevant_index=2

相關文章