From 03935705ccbdc340c18dba6fad79203591b1cad3 Mon Sep 17 00:00:00 2001
From: Mandar Chitre <mandar@nus.edu.sg>
Date: Fri, 6 Sep 2024 17:02:39 +0200
Subject: [PATCH 1/2] fix(dsp): circcorr was mistakenly called circconv (#44)

---
 src/dsp.jl         | 28 ++++++++++++++++------------
 test/tests-core.jl | 11 ++++++++---
 2 files changed, 24 insertions(+), 15 deletions(-)

diff --git a/src/dsp.jl b/src/dsp.jl
index e6aa176..06c9aed 100644
--- a/src/dsp.jl
+++ b/src/dsp.jl
@@ -2,10 +2,11 @@ import DSP: DSP, filt, filtfilt, resample, nextfastfft
 import Statistics: std
 import Peaks: argmaxima, peakproms!
 import Optim: optimize, minimizer, BFGS
+import FFTW: fft, ifft
 
 export fir, removedc, removedc!, demon
 export upconvert, downconvert, rrcosfir, rcosfir
-export mseq, gmseq, circconv, goertzel, pll, hadamard
+export mseq, gmseq, circconv, circcorr, goertzel, pll, hadamard
 export mfilter, findsignal, dzt, idzt
 export istft, whiten, filt, filtfilt, resample, delay!, compose
 
@@ -329,16 +330,19 @@ $(SIGNATURES)
 Computes the circular convolution of `x` and `y`. Both vectors must be the same
 length.
 """
-function circconv(x::AbstractVector, y::AbstractVector=x)
-  if length(x) != length(y)
-    throw(ArgumentError("x and y must be of equal length"))
-  end
-  n = length(x)
-  z = similar(x)
-  for j ∈ 1:n
-    z[j] = circshift(x, j-1)' * y
-  end
-  return z
+function circconv(x::AbstractVector, y::AbstractVector)
+  length(x) == length(y) || throw(ArgumentError("x and y must be of equal length"))
+  ifft(fft(x) .* fft(y))
+end
+
+"""
+$(SIGNATURES)
+Computes the circular correlation of `x` and `y`. Both vectors must be the same
+length.
+"""
+function circcorr(x::AbstractVector, y::AbstractVector=x)
+  length(x) == length(y) || throw(ArgumentError("x and y must be of equal length"))
+  ifft(conj.(fft(x)) .* fft(y))
 end
 
 """
@@ -739,7 +743,7 @@ function findsignal(r, s, n=1; prominence=0.0, finetune=2, mingap=1, mfo=false)
     return (time=Float64[], amplitude=T[], mfo=mfo ? m : empty(m))
   end
   h = m̄[p]
-  ndx = partialsortperm(h, 1:n; rev=true)
+  ndx = partialsortperm(h, 1:min(n,length(h)); rev=true)
   p = p[ndx]
   if finetune == 0
     t = time(Float64.(p), s)
diff --git a/test/tests-core.jl b/test/tests-core.jl
index bfbd7de..8814180 100644
--- a/test/tests-core.jl
+++ b/test/tests-core.jl
@@ -475,20 +475,25 @@ function test_dsp()
   @test y ≈ vcat(zeros(22), 1.0, zeros(22)) atol=0.01
   @test sum(x.^2) ≈ 1.0
 
+  x = rand(rng, 256)
+  @test circcorr(x) ≈ circconv(x, conj.(circshift(reverse(x), 1)))
+
+  x = rand(rng, ComplexF64, 256)
+  @test circcorr(x) ≈ circconv(x, conj.(circshift(reverse(x), 1)))
+
   for j ∈ 2:10
     x = mseq(j)
     @test x isa Array{Float64}
     @test length(x) == 2^j-1
     @test all(abs.(x) .== 1)
-    y = circconv(x)
-    @test all(y .== circconv(x, x))
+    y = circcorr(x)
     @test y[1] ≈ length(x)
     @test all(y[2:end] .≈ -1)
     x = gmseq(j)
     @test x isa Array{ComplexF64}
     @test length(x) == 2^j-1
     @test all(abs.(x) .== 1)
-    y = circconv(x)
+    y = circcorr(x)
     @test y[1] ≈ length(x)
     @test rms(y[2:end]) ≈ 0 atol=1e-10
   end

From 355e682f60d62bb933277af36a28d1159b1368bd Mon Sep 17 00:00:00 2001
From: Mandar Chitre <mandar@nus.edu.sg>
Date: Mon, 9 Sep 2024 14:30:53 +0800
Subject: [PATCH 2/2] test(dsp): add more test cases for circconv and bump
 version

---
 Project.toml       |  2 +-
 test/tests-core.jl | 10 ++++++++--
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2338c60..7a3e64b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SignalAnalysis"
 uuid = "df1fea92-c066-49dd-8b36-eace3378ea47"
 authors = ["Mandar Chitre <mandar@nus.edu.sg>"]
-version = "0.9.0"
+version = "0.10.0"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
diff --git a/test/tests-core.jl b/test/tests-core.jl
index 8814180..ecfb07c 100644
--- a/test/tests-core.jl
+++ b/test/tests-core.jl
@@ -476,10 +476,16 @@ function test_dsp()
   @test sum(x.^2) ≈ 1.0
 
   x = rand(rng, 256)
-  @test circcorr(x) ≈ circconv(x, conj.(circshift(reverse(x), 1)))
+  y = rand(rng, 256)
+  @test circcorr(x) ≈ circconv(conj.(circshift(reverse(x), 1)), x)
+  @test circcorr(x, y) ≈ circconv(conj.(circshift(reverse(x), 1)), y)
+  @test circconv(x, y) ≈ circconv(y, x)
 
   x = rand(rng, ComplexF64, 256)
-  @test circcorr(x) ≈ circconv(x, conj.(circshift(reverse(x), 1)))
+  y = rand(rng, ComplexF64, 256)
+  @test circcorr(x) ≈ circconv(conj.(circshift(reverse(x), 1)), x)
+  @test circcorr(x, y) ≈ circconv(conj.(circshift(reverse(x), 1)), y)
+  @test circconv(x, y) ≈ circconv(y, x)
 
   for j ∈ 2:10
     x = mseq(j)