Fix circonv implementation and add circcorr

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"

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)

test/tests-core.jl

@@ -475,20 +475,31 @@ 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)
+  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)
+  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)
     @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