Skip to content

Commit

Permalink
feat:add modes to conv function (JuliaDSP#399)
Browse files Browse the repository at this point in the history
  • Loading branch information
@Longhao-Chen
Longhao-Chen committed Feb 22, 2021
1 parent 99c01c3 commit 946464b
Showing 1 changed file with 45 additions and 13 deletions.
58 changes: 45 additions & 13 deletions src/dspbase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -679,13 +679,31 @@ end

# May switch argument order
"""
conv(u,v)
conv(u,v; mode = :full)
Convolution of two arrays. Uses either FFT convolution or overlap-save,
depending on the size of the input. `u` and `v` can be N-dimensional arrays,
with arbitrary indexing offsets, but their axes must be a `UnitRange`.
:full — Return the full 2-D convolution.
:same — Return the central part of the convolution, which is the same size as u.
:valid — Return only parts of the convolution that are computed without zero-padded edges.
"""
function conv(u::AbstractArray{T, N},
function conv(u, v; mode::Symbol = :full)
if mode == :full
conv_full(u, v)
elseif mode == :same
conv_same(u, v)
elseif mode == :valid
conv_valid(u, v)
else
throw(ArgumentError("mode keyword argument must be either :full or :same or :valid"))
end
end

function conv_full(u::AbstractArray{T, N},
v::AbstractArray{T, N}) where {T<:BLAS.BlasFloat, N}
su = size(u)
sv = size(v)
Expand All @@ -696,38 +714,52 @@ function conv(u::AbstractArray{T, N},
end
end

function conv(u::AbstractArray{<:BLAS.BlasFloat, N},
function conv_full(u::AbstractArray{<:BLAS.BlasFloat, N},
v::AbstractArray{<:BLAS.BlasFloat, N}) where N
fu, fv = promote(u, v)
conv(fu, fv)
conv_full(fu, fv)
end

conv(u::AbstractArray{<:Integer, N}, v::AbstractArray{<:Integer, N}) where {N} =
conv_full(u::AbstractArray{<:Integer, N}, v::AbstractArray{<:Integer, N}) where {N} =
round.(Int, conv(float(u), float(v)))

conv(u::AbstractArray{<:Number, N}, v::AbstractArray{<:Number, N}) where {N} =
conv_full(u::AbstractArray{<:Number, N}, v::AbstractArray{<:Number, N}) where {N} =
conv(float(u), float(v))

function conv(u::AbstractArray{<:Number, N},
function conv_full(u::AbstractArray{<:Number, N},
v::AbstractArray{<:BLAS.BlasFloat, N}) where N
conv(float(u), v)
conv_full(float(u), v)
end

function conv(u::AbstractArray{<:BLAS.BlasFloat, N},
function conv_full(u::AbstractArray{<:BLAS.BlasFloat, N},
v::AbstractArray{<:Number, N}) where N
conv(u, float(v))
conv_full(u, float(v))
end

function conv(A::AbstractArray{<:Number, M},
function conv_full(A::AbstractArray{<:Number, M},
B::AbstractArray{<:Number, N}) where {M, N}
if (M < N)
conv(cat(A, dims=N)::AbstractArray{eltype(A), N}, B)
conv_full(cat(A, dims=N)::AbstractArray{eltype(A), N}, B)
else
@assert M > N
conv(A, cat(B, dims=M)::AbstractArray{eltype(B), M})
conv_full(A, cat(B, dims=M)::AbstractArray{eltype(B), M})
end
end

function conv_same(u::AbstractArray{T, N}, v::AbstractArray{T, N}) where {T, N}
su = size(u)
sv = size(v)
conv_res = conv_full(u, v)
conv_res[Int(floor.(sv[1]/2 + 1)):Int(floor.(sv[1]/2) + su[1]), Int(floor.(sv[2]/2 + 1)):Int(floor.(sv[2]/2) + su[2]), axes(conv_res)[3:end]...]
end

function conv_valid(u::AbstractArray{T}, v::AbstractArray{T})::AbstractArray{T} where {T <: Number}
su = size(u)
sv = size(v)
conv_res = conv_full(u, v)
conv_res[sv[1]:su[1], sv[2]:su[2], axes(conv_res)[3:end]...]

This comment has been minimized.

Sign in to view

Copy link
@danilo-bc

danilo-bc Sep 16, 2021 •
edited
Loading

I'm interested in the discussion of JuliaDSP#399, so I came to check this commit.

I've checked the Matlab definition for valid mode. They state

Only those parts of the convolution that are computed without the zero-padded edges. Using this option, length(w) is max(length(u)-length(v)+1,0), except when length(v) is zero. If length(v) = 0, then length(w) = length(u).

Here (https://www.mathworks.com/help/matlab/ref/conv.html)

This seems weird, since conv(u, v) == conv(v, u). I think this line has the same pitfall. We can see that if su = 2 and sv = 1, su - sv != sv - su, which means indexing will go wrong and return an empty result.

Please confirm this behavior. The solution would probably be to have su and sl representing upper and lower bounds according to which one of su or sv is greater than the other.

This comment has been minimized.

Sign in to view

Copy link
@Longhao-Chen

Longhao-Chen via email Sep 21, 2021

Author Owner

This comment has been minimized.

Sign in to view

Copy link
@danilo-bc

danilo-bc Sep 21, 2021

I suggest checking once again. This is a simple test I made using python and Numpy

from numpy import convolve as conv
import numpy as np
a = np.array([1, 3, 5, -10, 2])
b = np.array([-30, 2, -3])
conv(a, b, mode='valid')
# output from REPL: array([-147,  301,  -95])
conv(b, a, mode='valid')
# output from REPL: array([-147,  301,  -95])

What was your test vector on Octave? Can you reproduce my results on Octave? Please warn me if I did something wrong.

Best regards!

This comment has been minimized.

Sign in to view

Copy link
@Longhao-Chen

Longhao-Chen via email Sep 21, 2021

Author Owner

This comment has been minimized.

Sign in to view

Copy link
@danilo-bc

danilo-bc Sep 22, 2021

Great point!

Somehow, Numpy's implementation is the odd one out.

I asked a friend to test the same behavior on Matlab and I tested locally on Scilab (another FOSS alternative, like Octave). Both yield the same behavior as you've shown.

In this case, I think your implementation is replicating this desired behavior. If someone wants, replicate Numpy's behavior by picking the 1st input as the largest of the two (in their desired dimension).
end

"""
conv(u,v,A)
Expand Down

0 comments on commit 946464b

Please sign in to comment.