Bug detected in deconvolution operation deconv()


The DSP package contains some routines very vital and basic to Fourier analysis. However the simple deconvultion operation is not working. Consider the following simple code

using DSP;
a = zeros( Float64, 5); a[3]=10;
b = zeros( Float64, 5); b[3]=0.1;
c = conv(a,b);
d = deconv(c,a); display(d');

Thus a, b represent the Fourier basis of simple Fourier functions and the coefficients are inverses of each other. The vector c is thus all zeros with one 1. Deconvoluting against c however doesnt give b. It gives the error :

ERROR: ArgumentError: filter vector a[1] must be nonzero
 [1] filt!(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\iamsu\.julia\packages\DSP\RUXAo\src\dspbase.jl:34
 [2] deconv(::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\iamsu\.julia\packages\DSP\RUXAo\src\dspbase.jl:18
 [3] top-level scope at none:0

Since the mathematics is simple and correct, this reveals a bug in the routine.

I can’t comment on the bug or the maths, but I’d recommend posting this as an issue on https://github.com/JuliaDSP/DSP.jl/issues directly.


Indeed, the deconv function does not produce polynomial division, as advertised, if the denominator has its first coefficient = 0. The function divrem produces a correct result in dividing c by a:

using Polynomials, DSP
a = [0, 0, 10.0, 0, 0]
b = [0, 0, 0.1, 0, 0]
c = conv(a,b)
A = Polynomial(a)
C = Polynomial(c)
D, r = divrem(C,A);
d = coeffs(D)

Thats a nice fix. I am doing it explicitly while using the Shannon-Nyquist sampling frequency.

Thanks, I will try to.