LPC burg: ar_burg() interface - returns reflection coeffs - same as … #171
Conversation
…matlab call as well, lpc() - call - returns 2 params - same as Levinson version
|
Hi, I wrote the function @staticfloat integrated here before (see: #156 ) - this is great, but to be able to use this version in my own code - I also need the Reflection-Coefficients as a return value (they can be re-evaluated - but there's a numerical stability cost) |
|
note: the non passing CI test seems to be due to some non-related incompatibility issue with Julia nightly (0.7 line) - something about some functions to have moved out of base into DSP |
src/lpc.jl
Outdated
|
|
||
| # Dispatch types for lpc() | ||
| mutable struct LPCBurg; end | ||
| mutable struct LPCLevinson; end | ||
|
|
||
| """ | ||
| `lpc_burg(x::AbstractVector, p::Int)` | ||
| `ar_burg(x::AbstractVector, p::Int)` |
There was a problem hiding this comment.
this should be arburg() not ar_burg()
| @@ -59,7 +66,7 @@ function implements the mathematics described in [1]. | |||
| """ | |||
| function lpc{T <: Number}(x::AbstractVector{T}, p::Int, ::LPCLevinson) | |||
| R_xx = xcorr(x,x)[length(x):end] | |||
| a = zeros(p,p) | |||
| a = zeros(T,p,p) | |||
There was a problem hiding this comment.
I've tested the versions for type stability - the Burg versions works with Float(16/32/64)/BigFloat - where as the Levinson version works only with Float(32/64) - because of conv() doing an FFT (which is alright for now I guess).
Anyhow - the Levinson version can also have the reflection coeffs returns eventually.
Test is cool, I'll try to add a Ladder(Lattice) filter - where the reflection coeffs are used, and that can also be added to the test
There was a problem hiding this comment.
Is this accurate for integer T?
There was a problem hiding this comment.
I would try to add a test with a Fixed-Point implementation ( like https://github.com/JuliaMath/FixedPointNumbers.jl ) the accepted base time might have to be "Real" , "Integer" as is - should probably be excluded, but stability of Rational + FixedPoint might be nice to have. A 'trait' of "invertible" (or 'semi-field' or anything similar ) might be required.
There was a problem hiding this comment.
Ah I wasn't even thinking about fixed point arithmetic. That would be cool.
My point with these comments is that the current PR will error if x is an abstract vector with Integer elements, because of the type assertions and container types used in the code. I think you should either change the type assertions and container types to reflect the actual types, or restrict the type bounds in the method definitions to types that will work as you expect (e.g. AbstractFloat).
Because restricting the type bounds will change the API for this function, I think it would be best to be more careful about the container types. If T = eltype(x) and T <: AbstractFloat or T <: FixedPoint then prediciton_error, reflection_coeffs, and a will be type T (or have eltype T). If T <: Integer then prediction_error and reflection_coeffs will be type Float64, and the Levinson algorithm will error because xcorr errors with integer inputs.
There was a problem hiding this comment.
Fixed-Point is quite common for LPC implementations, it's an import use-case as LPC is used in some analysis and compression algorithms in some embedded systems with Integer-Only hardware (still I believe); the few integer LPC's I've seen had the Fixed-Point arithmetics inlined by hand (for example in Speex), but this could actually be an interesting use case to see whether we can produce efficient machine code with an abstract type here (might change the abstraction for that). Anyhow, I think I've even tested this with Rational before, it's quite cool that it works.
src/lpc.jl
Outdated
| # Initialize prediction error wtih the variance of the signal | ||
| prediction_err = (sum(abs2, x) ./ length(x))[1] | ||
| prediction_err :: T = sum(abs2, x) ./ length(x) |
There was a problem hiding this comment.
/ instead of ./? Or is there some reason for the broadcasting I'm not aware of?
There was a problem hiding this comment.
no, change it to / if you see fit.
This was not changed in the current commit, I probably wrote it as to be explicit about doing a scalar operation (rather than a broadcast)
| # Initialize prediction error wtih the variance of the signal | ||
| prediction_err = (sum(abs2, x) ./ length(x))[1] | ||
| prediction_err :: T = sum(abs2, x) / length(x) |
There was a problem hiding this comment.
Is prediction_err always type T? What if T is an integer? It may be good to change the tests to include integer data.
| # Initialize prediction error wtih the variance of the signal | ||
| prediction_err = (sum(abs2, x) ./ length(x))[1] | ||
| prediction_err :: T = sum(abs2, x) / length(x) | ||
| reflection_coeffs = zeros(T, p) |
There was a problem hiding this comment.
If x is a vector of integers, then it seems like reflections_coeffs would not also be integers, would they? I'm not sure if T is the right type of reflection_coeffs.
| end | ||
|
|
||
| function lpc{T <: Number}(x::AbstractVector{T}, p::Int, ::LPCBurg) | ||
| a, prediction_err, _ = arburg(x, p) |
There was a problem hiding this comment.
I have added _ to take care of the reflection_coeffs return value here. Is that ok?
| @@ -59,7 +66,7 @@ function implements the mathematics described in [1]. | |||
| """ | |||
| function lpc{T <: Number}(x::AbstractVector{T}, p::Int, ::LPCLevinson) | |||
| R_xx = xcorr(x,x)[length(x):end] | |||
| a = zeros(p,p) | |||
| a = zeros(T,p,p) | |||
There was a problem hiding this comment.
Is this accurate for integer T?
|
Superseded by #517, although more tests e.g. for |
…matlab call as well, lpc() - call - returns 2 params - same as Levinson version