feat: add Hadamard codes

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

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