diff --git a/dsp/transform/doc.go b/dsp/transform/doc.go new file mode 100644 index 00000000..bb4f0b02 --- /dev/null +++ b/dsp/transform/doc.go @@ -0,0 +1,6 @@ +// Copyright ©2024 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +// Package transform provides important transforms on signals used in digital signal processing. +package transform // import "gonum.org/v1/gonum/dsp/transform" diff --git a/dsp/transform/hilbert.go b/dsp/transform/hilbert.go new file mode 100644 index 00000000..bf444ee1 --- /dev/null +++ b/dsp/transform/hilbert.go @@ -0,0 +1,76 @@ +// Copyright ©2024 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package transform + +import ( + "gonum.org/v1/gonum/dsp/fourier" +) + +// Hilbert implements an approximate Hilbert transform that allows calculation +// of an approximate analytical signal of a real signal, and determine the +// real envelope of a signal. +type Hilbert struct { + fft *fourier.CmplxFFT + work []complex128 +} + +// NewHilbert returns a new Hilbert transformer for signals of size n. +// The transform is most efficient when n is a product of small primes. +// n must not be less than one. +func NewHilbert(n int) *Hilbert { + return &Hilbert{ + fft: fourier.NewCmplxFFT(n), + work: make([]complex128, n), + } +} + +// Len returns the length of signals that are valid input for this Hilbert transform. +func (h *Hilbert) Len() int { + return len(h.work) +} + +// AnalyticSignal computes the analytical signal of a real signal, and stores +// the result in the dst slice, returning it. +// +// If the dst slice is nil, a new slice will be created and returned. The dst slice +// must be the same length as the input signal, otherwise the method will panic. +func (h *Hilbert) AnalyticSignal(dst []complex128, signal []float64) []complex128 { + if len(signal) != h.Len() { + panic("transform: input signal length mismatch") + } + + if dst == nil { + dst = make([]complex128, len(signal)) + } else if len(dst) != h.Len() { + panic("transform: destination length mismatch") + } + + for i, v := range signal { + h.work[i] = complex(v, 0) + } + + // Forward FFT of the signal. + coeff := h.fft.Coefficients(dst, h.work) + for i := range h.work { + h.work[i] = 0 + } + + // Multiply positive frequencies by 2, zero out negative frequencies. + // However, leave dc unchanged (and nyquist when n%2 == 0). + h.work[0] = coeff[0] + for i, d := range coeff[1 : len(coeff)/2+1] { + h.work[i+1] = d * 2 + } + if len(coeff)%2 == 0 { + h.work[len(coeff)/2] = coeff[len(coeff)/2] + } + + // Normalize the results so they have a similar amplitude to the input + unnorm := h.fft.Sequence(dst, h.work) + for i, u := range unnorm { + unnorm[i] = u / complex(float64(len(unnorm)), 0) + } + return unnorm +} diff --git a/dsp/transform/hilbert_example_test.go b/dsp/transform/hilbert_example_test.go new file mode 100644 index 00000000..2f571dd4 --- /dev/null +++ b/dsp/transform/hilbert_example_test.go @@ -0,0 +1,56 @@ +// Copyright ©2024 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package transform + +import ( + "fmt" + "math/cmplx" +) + +func ExampleHilbert_AnalyticSignal() { + // Samples is a set of real amplitudes that make up a signal. + samples := []float64{1, 0, 2, 0, 4, 0, 2, 0} + + // Initialize a Hilbert transform and 'demodulate' to get the + // analytic signal. + // The result is the complex I/Q (In-Phase / Quadrature) demodulation + // of the input signal. + h := NewHilbert(len(samples)) + iqSamples := h.AnalyticSignal(nil, samples) + + // We can compute the instantaneous amplitude of the signal + // (or 'envelope') using absolute value. Analyzing the envelope + // is an easy way to measure changes in amplitude over time in a + // signal. + envelope := make([]float64, len(samples)) + for ind, iq := range iqSamples { + envelope[ind] = cmplx.Abs(iq) + } + + // We can also compute the instantaneous phase of each part of the + // signal using the 4-quadrant arc-tangent. With multiple samples, + // the instantaneous phase can be used to estimate instantaneous + // frequency of a signal. + phase := make([]float64, len(samples)) + for ind, iq := range iqSamples { + phase[ind] = cmplx.Phase(iq) + } + + for i, iq := range iqSamples { + fmt.Printf("ind=%d -> I=%.4f, Q=%.4f, envelope=%.4f, phase=%.4f\n", + i, real(iq), imag(iq), envelope[i], phase[i]) + } + + // Output: + // + // ind=0 -> I=1.0000, Q=0.0000, envelope=1.0000, phase=0.0000 + // ind=1 -> I=-0.0000, Q=-0.8107, envelope=0.8107, phase=-1.5708 + // ind=2 -> I=2.0000, Q=0.0000, envelope=2.0000, phase=0.0000 + // ind=3 -> I=-0.0000, Q=-1.3107, envelope=1.3107, phase=-1.5708 + // ind=4 -> I=4.0000, Q=0.0000, envelope=4.0000, phase=0.0000 + // ind=5 -> I=0.0000, Q=1.3107, envelope=1.3107, phase=1.5708 + // ind=6 -> I=2.0000, Q=0.0000, envelope=2.0000, phase=0.0000 + // ind=7 -> I=0.0000, Q=0.8107, envelope=0.8107, phase=1.5708 +} diff --git a/dsp/transform/hilbert_test.go b/dsp/transform/hilbert_test.go new file mode 100644 index 00000000..3a671cdb --- /dev/null +++ b/dsp/transform/hilbert_test.go @@ -0,0 +1,57 @@ +// Copyright ©2024 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package transform + +import ( + "testing" + + "gonum.org/v1/gonum/cmplxs" +) + +var hilbertAnalyticSignalTests = []struct { + name string + in []float64 + want []complex128 +}{ + { + name: "zeros", + in: []float64{0, 0, 0, 0}, + want: []complex128{0, 0, 0, 0}}, + { + name: "whole_components", + in: []float64{1, 2, 3, 4}, + want: []complex128{1 + 1i, 2 - 1i, 3 - 1i, 4 + 1i}, + }, + { + name: "irrational_imaginary_components", + in: []float64{1, 2, 3, 4, 5}, + want: []complex128{ + 1 + 1.7013016167i, + 2 - 1.3763819204i, + 3 - 0.6498393924i, + 4 - 1.3763819204i, + 5 + 1.7013016167i, + }, + }, +} + +func TestHilbertAnalytic(t *testing.T) { + const tol = 1e-10 + + for _, test := range hilbertAnalyticSignalTests { + t.Run(test.name, func(t *testing.T) { + h := NewHilbert(len(test.in)) + if h.Len() != len(test.in) { + t.Errorf("unexpected Hilbert transform length: got:%d, want:%d", h.Len(), len(test.in)) + } + + dst := make([]complex128, len(test.in)) + got := h.AnalyticSignal(dst, test.in) + if !cmplxs.EqualApprox(got, test.want, tol) { + t.Errorf("unexpected Hilbert transform result:\ngot: %v\nwant:%v", got, test.want) + } + }) + } +}