|
| 1 | +// created by christian privitelli on 20/10/2023 |
| 2 | + |
| 3 | +import Accelerate |
| 4 | +import Foundation |
| 5 | + |
| 6 | +public class LinearPredictiveCoding { |
| 7 | + public func computeLpc(_ y: [Float], order: Int) -> [Float] { |
| 8 | + var arCoeffs = [Float](repeating: 0.0, count: order + 1) |
| 9 | + var arCoeffsPrev = [Float](repeating: 0.0, count: order + 1) |
| 10 | + |
| 11 | + arCoeffs[0] = 1 |
| 12 | + arCoeffsPrev[0] = 1 |
| 13 | + |
| 14 | + // forward and backward prediction error. used to calculate reflection coefficient. |
| 15 | + var fwdError = Array(y.dropFirst()) |
| 16 | + var bwdError = Array(y.dropLast()) |
| 17 | + |
| 18 | + var den = computeInitialDen(fwdError: fwdError, bwdError: bwdError) |
| 19 | + var reflectionCoefficient: Float = 0 |
| 20 | + |
| 21 | + for i in 0 ..< order { |
| 22 | + reflectionCoefficient = computeReflectionCoefficient(fwdError: fwdError, bwdError: bwdError, den: den) |
| 23 | + |
| 24 | + swap(&arCoeffsPrev, &arCoeffs) |
| 25 | + |
| 26 | + for j in 1 ... (i + 1) where j < arCoeffsPrev.count && (i - j + 1) < arCoeffsPrev.count { |
| 27 | + var values: [Float] = [arCoeffsPrev[j], arCoeffsPrev[i - j + 1]] |
| 28 | + var result: Float = 0.0 |
| 29 | + |
| 30 | + vDSP_dotpr(&values, 1, [1.0, reflectionCoefficient], 1, &result, vDSP_Length(2)) |
| 31 | + |
| 32 | + arCoeffs[j] = result |
| 33 | + } |
| 34 | + |
| 35 | + computePredictionErrors(fwdError: &fwdError, bwdError: &bwdError, reflectionCoefficient: &reflectionCoefficient) |
| 36 | + |
| 37 | + den = computeUpdatedDen( |
| 38 | + fwdError: fwdError, |
| 39 | + bwdError: bwdError, |
| 40 | + reflectionCoefficient: reflectionCoefficient, |
| 41 | + den: den |
| 42 | + ) |
| 43 | + |
| 44 | + // shift up forward error |
| 45 | + fwdError = Array(fwdError.dropFirst()) |
| 46 | + bwdError = Array(bwdError.dropLast()) |
| 47 | + } |
| 48 | + |
| 49 | + return arCoeffs |
| 50 | + } |
| 51 | + |
| 52 | + private func computeInitialDen(fwdError: [Float], bwdError: [Float]) -> Float { |
| 53 | + var fwdSq: Float = 0.0 |
| 54 | + var bwdSq: Float = 0.0 |
| 55 | + vDSP_dotpr(fwdError, 1, fwdError, 1, &fwdSq, vDSP_Length(fwdError.count)) |
| 56 | + vDSP_dotpr(bwdError, 1, bwdError, 1, &bwdSq, vDSP_Length(bwdError.count)) |
| 57 | + return fwdSq + bwdSq |
| 58 | + } |
| 59 | + |
| 60 | + private func computeUpdatedDen(fwdError: [Float], bwdError: [Float], reflectionCoefficient: Float, den: Float) -> Float { |
| 61 | + let q = 1.0 - reflectionCoefficient * reflectionCoefficient |
| 62 | + let lastElementBwd = bwdError.last! * bwdError.last! |
| 63 | + let firstElementFwd = fwdError.first! * fwdError.first! |
| 64 | + return q * den - (lastElementBwd + firstElementFwd) |
| 65 | + } |
| 66 | + |
| 67 | + private func computeReflectionCoefficient(fwdError: [Float], bwdError: [Float], den: Float) -> Float { |
| 68 | + let epsilon = Float.leastNonzeroMagnitude |
| 69 | + |
| 70 | + var result: Float = 0.0 |
| 71 | + vDSP_dotpr(bwdError, 1, fwdError, 1, &result, vDSP_Length(bwdError.count)) |
| 72 | + return -2 * result / (den + epsilon) |
| 73 | + } |
| 74 | + |
| 75 | + private func computePredictionErrors(fwdError: inout [Float], bwdError: inout [Float], reflectionCoefficient: inout Float) { |
| 76 | + // We need to pass the state before fwdError is modified in the first vDSP_vsma into the second vDSP_vsma that modifies bwdError. |
| 77 | + let fwdErrorTemp = fwdError |
| 78 | + |
| 79 | + /// vDSP_vsma |
| 80 | + // for (n = 0; n < N; ++n) |
| 81 | + // D[n] = A[n]*B + C[n]; |
| 82 | + vDSP_vsma(bwdError, 1, &reflectionCoefficient, fwdError, 1, &fwdError, 1, vDSP_Length(fwdError.count)) |
| 83 | + vDSP_vsma(fwdErrorTemp, 1, &reflectionCoefficient, bwdError, 1, &bwdError, 1, vDSP_Length(bwdError.count)) |
| 84 | + } |
| 85 | +} |
0 commit comments