From 23fa31767f9c2363520f840e75d77d799f555219 Mon Sep 17 00:00:00 2001
From: Mandar Chitre <mandar@nus.edu.sg>
Date: Thu, 6 Jun 2024 16:42:49 +0530
Subject: [PATCH] feat: add Hadamard codes

---
 src/dsp.jl         | 26 +++++++++++++++++++++++++-
 test/tests-core.jl |  9 +++++++++
 2 files changed, 34 insertions(+), 1 deletion(-)

diff --git a/src/dsp.jl b/src/dsp.jl
index 5d10191..19b9892 100644
--- a/src/dsp.jl
+++ b/src/dsp.jl
@@ -5,7 +5,7 @@ import Optim: optimize, minimizer, BFGS
 
 export fir, removedc, removedc!, demon
 export upconvert, downconvert, rrcosfir, rcosfir
-export mseq, gmseq, circconv, goertzel, pll
+export mseq, gmseq, circconv, goertzel, pll, hadamard
 export sfilt, sfiltfilt, sresample, mfilter, findsignal
 export istft, whiten, filt, filtfilt, resample, delay!, compose
 
@@ -296,6 +296,30 @@ function gmseq(m, θ=atan(√(2^maximum(m)-1)))
   cos(θ) .+ 1im * sin(θ) .* x
 end
 
+"""
+    hadamard(k)
+
+Generate a Walsh-Hadamard matrix of size `2ᵏ × 2ᵏ`. Each row of the matrix
+is orthogonal to all other rows.
+"""
+function hadamard(k)
+  n = 2^k
+  [(-1)^count_ones(x&y) for x in 0:n-1, y in 0:n-1]
+end
+
+"""
+    hadamard(i, k)
+
+Generate a vector with the entries of row `i` of a Walsh-Hadamard matrix of
+size `2ᵏ × 2ᵏ`. Rows are numbered from `0` to `2ᵏ-1`, so that `i = 1` is the
+first non-trivial (not all ones) Hadamard sequence.
+"""
+function hadamard(i, k)
+  n = 2^k
+  0 ≤ i < n || throw(ArgumentError("i must be in the range 0 to 2^k-1"))
+  [(-1)^count_ones(x&i) for x in 0:n-1]
+end
+
 """
 $(SIGNATURES)
 Computes the circular convolution of `x` and `y`. Both vectors must be the same
diff --git a/test/tests-core.jl b/test/tests-core.jl
index 063adf8..9645c4c 100644
--- a/test/tests-core.jl
+++ b/test/tests-core.jl
@@ -498,6 +498,15 @@ function test_dsp()
   @test mseq((1,3)) == mseq(3)
   @test all(gmseq(3, 0) .== 1.0)
 
+  for j ∈ [1,2,4,8]
+    H = hadamard(j)
+    @test H isa Matrix{Int}
+    @test H' * H == 2^j * I(2^j)
+    for i ∈ 0:2^j-1
+      @test hadamard(i, j) == H[i+1,:]
+    end
+  end
+
   x = cw(10kHz, 0.1s, 80kHz)
   @test abs(goertzel(x, 10kHz)) ≈ length(x) atol=1e-6
   @test abs(goertzel(x, 9kHz)) ≈ 0 atol=1e-6