Bug detected in deconvolution operation deconv()

Specific Domains Signal and Image Processing
,
1

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
Stacktrace:
 [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.

2

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.

2 Likes
3

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)
3 Likes
4

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

5

Thanks, I will try to.