swift - 如何使三次样条插值算法更快？

Question

我意识到，当我使用“la_solve”时，由于“la_solve”可以处理的样本数量有限，我会收到 Malloc 错误。

我注意到如果我使用大约 900 个样本，程序运行良好。

因此，我修改了我的代码，以在样本 > 901 的情况下查看 901 个样本的连续 block 中的整个数据。

但是，我使用了相当幼稚的技术，并且确信这可以变得更加高效和稳健:(1) 我需要一个条件，如果样本小于 900，则不要分块； (2) 在附加分块数据时注意 x 和 y 值的重复； (3) 使算法更加高效、快速和稳健。

我最需要帮助的是代码末尾的“//调用重新采样器功能”部分

谢谢!!!

import UIKit
import Accelerate.vecLib.LinearAlgebra

public class Resample {

public class func resample(xi: [Double], xii: [Double], a: [Double]) ->[Double]
{
    // xi - original time stamps (original x)
    // xii - desired time stamps (desired x)
    // a - orignal y values
    // ->[Double] - interpolated y values

    let ni = xii.count // Desired x count

    let n = xi.count // Actual x count

    var h: [Double] = Array(repeating:0, count:n-1)

    for j in 0..<n-1 {
        h[j] = (xi[j+1] - xi[j])
    }

    var A: [Double] = Array(repeating:0, count:(n)*(n))

    A[0] = 1.0

    A[(n*n)-1] = 1.0

    for i in 1..<(n-1) {
        A[(n+1)*i-1] = h[i-1]
        A[(n+1)*i] = 2*(h[i-1] + h[i])
        A[(n+1)*i+1] = h[i]
    }

    var b: [Double] = Array(repeating:0, count:n)

    for i in 1..<n-1 {
        b[i] = (3/h[i])*(a[i+1]-a[i]) - (3/h[i-1])*(a[i]-a[i-1])
    }

    let matA = la_matrix_from_double_buffer(A, la_count_t(n), la_count_t(n), la_count_t(n), la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))

    let vecB = la_matrix_from_double_buffer(b, la_count_t(n), 1, 1, la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))

    let vecCj = la_solve(matA, vecB)

    var cj: [Double] = Array(repeating: 0.0, count: n)

    let status = la_matrix_to_double_buffer(&cj, 1, vecCj)

    if status == la_status_t(LA_SUCCESS) {
        //print(cj.count)
    }
    else {
        print("Failure: \(status)")
    }

    var bj: [Double] = Array(repeating:0, count:n-1)

    for i in 0..<n-1 {
        bj[i] = (1/h[i])*(a[i+1]-a[i]) - (1/3*h[i])*(2*cj[i]+cj[i+1])
    }

    var dj: [Double] = Array(repeating:0, count:n-1)

    for i in 0..<n-1 {
        dj[i] = (1/(3*h[i])) * (cj[i+1]-cj[i])
    }

    var P: [Double] = Array(repeating: 0.0, count: (n-1)*4)

    for i in 0..<n-1 {
        P[(i*4)] = dj[i]
        P[(i*4)+1] = cj[i]
        P[(i*4)+2] = bj[i]
        P[(i*4)+3] = a[i]
    }

    var ai: [Double] = Array(repeating: 0.0, count: ni)


    var jl = Double(1)

    for i in 0..<ni {

        let inter = Double(xii[i])

        while (inter > xi[Int(jl)] && (Int(jl) < n)) {
            jl += 1

        }

        let ind = 4*Int(jl)-4

        var pp = Array(P[ind...ind+3])

        let xx = Double(xii[i]-xi[Int(jl)-1])

        ai[i] = pp[0]*pow(xx, 3) + pp[1]*pow(xx, 2) + pp[2]*xx + pp[3]
    }
    return(ai)
}
}



 // Calling the Re-sampler FUNCTION

let x = Array(stride(from: 0, through: 16649, by: 1.0)) // Orig X (Given, say)

let y = Array(stride(from: 0, through: 16649, by: 1.0)) // Orig Y (Given, say)

let f = 2.0 // Hz (Desired Sampling Freq., say)

let g = Array(stride(from: 0, through: x.count-1, by: 900))
// Helper array for chunking (901 samples each time, say)

var xxx = [Double]() // Results for appended chunks X

var yyy = [Double]() // Results for appended chunks Y

for i in 0..<g.count-1 { // loop through

let xc = Array(x[Int(g[i])...Int(g[i+1])]) // x chunk (count = 901)

let yc = Array(y[Int(g[i])...Int(g[i+1])]) // y chunk (count = 901)

let xci = Array(stride(from: xc[0], through: xc[xc.count-1], by: 1.0/f))
// Desired time stamps for chunk of 901 samples

let yci = Resample.resample(xi: xc, xii: xci, a: yc) // Interpolate chunk

xxx += xci // Append interpolation X

yyy += yci // Append interpolation X

}

if(Int(g[g.count-1])<x.count){
// If helper function last index < original no. of samples

let glb = (Int(g[g.count-1])) // last chunk begin index

let gle = (x.count) // last chunk end index

let xce = Array(x[glb...gle-1]) // last x chunk

let yce = Array(y[glb...gle-1]) // last y chunk

let xcei = Array(stride(from: xce[0], through: xce[xce.count-1], by: 1.0/f))
// Desired time stamps for last chunk of samples

let ycei = Resample.resample(xi: xce, xii: xcei, a: yce) // Interpolate last chunk

xxx += xcei // Append interpolation X for last chunk

yyy += ycei // Append interpolation X for last chunk
}

print(xxx) // would like to identify repeated x values and remove them
print(yyy) // remove y values corresponsding to repated x values (as found above)

// Calling the Re-sampler FUNCTION

Answer 1

请注意，您正在使用二进制搜索寻找每个 ind，而它的值是单调递增的。所以你可以只检查下一个 x 值是否保持在相同的间隔内或移动到下一个。伪代码:

 jl = 1  //before cycle

 //inside cycle
  while (current_x > xi[jl]) and (jl < n)
    jl++
  ind = 4 * jl - 4

la_solve 是求解线性方程组的通用方法吗？如果是，那么您最好对三对角线系统使用专用方法(三次复杂度与线性复杂度)。 Classic text . Python implementation . Wiki .对于长输入数组很重要。