From dc40d47f3c7720ab4f1481bd1b032bf702499295 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Mon, 18 Mar 2024 23:28:29 +0800 Subject: [PATCH 1/5] clean 2d conv CartesianIndices(A) --- src/dspbase.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/dspbase.jl b/src/dspbase.jl index afb7c326..8987ea3f 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -794,13 +794,13 @@ function conv(u::AbstractVector{T}, v::Transpose{T,<:AbstractVector}, A::Abstrac if any(conv_axis_with_offset, (axes(u)..., axes(v)..., axes(A)...)) throw(ArgumentError("offset axes not supported")) end - m = length(u)+size(A,1)-1 - n = length(v)+size(A,2)-1 + m = length(u) + size(A, 1) - 1 + n = length(v) + size(A, 2) - 1 B = zeros(T, m, n) - B[1:size(A,1),1:size(A,2)] = A - u = fft([u;zeros(T,m-length(u))]) - v = fft([v transpose(zeros(T,n-length(v)))]) - C = ifft(fft(B) .* (u * v)) + B[CartesianIndices(A)] = A + u, v = fft.(_zeropad.((u, transpose(v)), (m, n))) + p = plan_fft(B) + C = inv(p) * ((p * B) .*= u .* transpose(v)) if T <: Real return real(C) end From 9bf0d9583922e5eba12780bd65353ea366cdc49e Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sun, 17 Mar 2024 22:54:43 +0800 Subject: [PATCH 2/5] remove inbounds from `unsafe_conv_kern_os!` slightly faster on v1.10, about the same on v1.6 --- src/dspbase.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dspbase.jl b/src/dspbase.jl index 8987ea3f..e7db6d3f 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -402,7 +402,7 @@ function unsafe_conv_kern_os_edge!( # on an edge. # # First make all entries equal to the center blocks: - @inbounds copyto!(perimeter_range, 1, center_block_ranges, 1, N) + copyto!(perimeter_range, 1, center_block_ranges, 1, N) # For the dimensions chosen to be on an edge (edge_dims), get the # ranges of the blocks that would need to be padded (lie on an edge) @@ -421,7 +421,7 @@ function unsafe_conv_kern_os_edge!( # The center region for non-edge dimensions has been specified above, # so finish specifying the region of the perimeter for this edge # block - @inbounds for (i, dim) in enumerate(edge_dims) + for (i, dim) in enumerate(edge_dims) perimeter_range[dim] = perimeter_edge_ranges[i] end @@ -429,7 +429,7 @@ function unsafe_conv_kern_os_edge!( block_region = CartesianIndices( NTuple{N, UnitRange{Int}}(perimeter_range) ) - @inbounds for block_pos in block_region + for block_pos in block_region # Figure out which portion of the input data should be transformed block_idx = convert(NTuple{N, Int}, block_pos) @@ -571,7 +571,7 @@ function unsafe_conv_kern_os!(out, # Portion of buffer with valid result of convolution valid_buff_region = CartesianIndices(UnitRange.(sv, nffts)) # Iterate over block indices (not data indices) that do not need to be padded - @inbounds for block_pos in CartesianIndices(center_block_ranges) + for block_pos in CartesianIndices(center_block_ranges) # Calculate portion of data to transform block_idx = convert(NTuple{N, Int}, block_pos) From 99f9247fc96221a252307f4287fb718808f5971c Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 11:55:33 +0800 Subject: [PATCH 3/5] vt for clarity --- src/dspbase.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dspbase.jl b/src/dspbase.jl index e7db6d3f..32f82fa1 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -798,9 +798,9 @@ function conv(u::AbstractVector{T}, v::Transpose{T,<:AbstractVector}, A::Abstrac n = length(v) + size(A, 2) - 1 B = zeros(T, m, n) B[CartesianIndices(A)] = A - u, v = fft.(_zeropad.((u, transpose(v)), (m, n))) + u, vt = fft.(_zeropad.((u, transpose(v)), (m, n))) p = plan_fft(B) - C = inv(p) * ((p * B) .*= u .* transpose(v)) + C = inv(p) * ((p * B) .*= u .* transpose(vt)) if T <: Real return real(C) end From 027b8ba50a0ee1452729b3b3fc9d747cc99964c2 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 14:38:19 +0800 Subject: [PATCH 4/5] Apply suggestion to avoid broadcasting Co-authored-by: Martin Holters --- src/dspbase.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dspbase.jl b/src/dspbase.jl index 32f82fa1..f53a49d2 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -798,7 +798,8 @@ function conv(u::AbstractVector{T}, v::Transpose{T,<:AbstractVector}, A::Abstrac n = length(v) + size(A, 2) - 1 B = zeros(T, m, n) B[CartesianIndices(A)] = A - u, vt = fft.(_zeropad.((u, transpose(v)), (m, n))) + u = fft(_zeropad(u, m)) + vt = fft(_zeropad(transpose(v), n)) p = plan_fft(B) C = inv(p) * ((p * B) .*= u .* transpose(vt)) if T <: Real From e8f87c90c06220744afda7a659dc6dc11addad57 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 16:52:08 +0800 Subject: [PATCH 5/5] don't plan fft Co-authored-by: Martin Holters --- src/dspbase.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/dspbase.jl b/src/dspbase.jl index f53a49d2..ddb75730 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -800,8 +800,7 @@ function conv(u::AbstractVector{T}, v::Transpose{T,<:AbstractVector}, A::Abstrac B[CartesianIndices(A)] = A u = fft(_zeropad(u, m)) vt = fft(_zeropad(transpose(v), n)) - p = plan_fft(B) - C = inv(p) * ((p * B) .*= u .* transpose(vt)) + C = ifft(fft(B) .*= u .* transpose(vt)) if T <: Real return real(C) end