From 880b710e3cb6a2fe5a8d51a9211cae3cc54ee93a Mon Sep 17 00:00:00 2001 From: Dan Kortschak Date: Tue, 15 Sep 2020 14:21:53 +0930 Subject: [PATCH] dsp/window: add Tukey window function and simplify tests --- dsp/window/window.go | 43 ++++++++++++++++ dsp/window/window_complex.go | 43 ++++++++++++++++ dsp/window/window_test.go | 95 +++++++++--------------------------- 3 files changed, 110 insertions(+), 71 deletions(-) diff --git a/dsp/window/window.go b/dsp/window/window.go index 44f50b01..da6366e1 100644 --- a/dsp/window/window.go +++ b/dsp/window/window.go @@ -330,6 +330,49 @@ func (g Gaussian) Transform(seq []float64) []float64 { return seq } +// Tukey can modify a sequence using the Tukey window and return the result. +// See https://en.wikipedia.org/wiki/Window_function#Tukey_window +// and https://www.recordingblogs.com/wiki/tukey-window for details. +// +// The Tukey window is an adjustable window. +// +// The sequence weights are +// w[k] = 0.5 * (1 + cos(π*(|k - M| - αM)/((1-α) * M))), |k - M| ≥ αM +// = 1, |k - M| < αM +// with M = (N - 1)/2 for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters are summarized in the table: +// | α=0.3 | α=0.5 | α=0.7 | +// -------|--------------------------| +// ΔF_0 | 1.33 | 1.22 | 1.13 | +// ΔF_0.5 | 1.28 | 1.16 | 1.04 | +// K | 0.67 | 0.61 | 0.57 | +// ɣ_max | -18.2 | -15.1 | -13.8 | +// β | -1.41 | -2.50 | -3.74 | +type Tukey struct { + Alpha float64 +} + +// Transform applies the Tukey transformation to seq in place, using the value +// of the receiver as the Alpha parameter, and returning the result +func (t Tukey) Transform(seq []float64) []float64 { + switch { + case t.Alpha <= 0: + return Rectangular(seq) + case t.Alpha >= 1: + return Hann(seq) + default: + alphaL := t.Alpha * float64(len(seq)-1) + width := int(0.5*alphaL) + 1 + for i := range seq[:width] { + w := 0.5 * (1 - math.Cos(2*math.Pi*float64(i)/alphaL)) + seq[i] *= w + seq[len(seq)-1-i] *= w + } + return seq + } +} + // Values is an arbitrary real window function. type Values []float64 diff --git a/dsp/window/window_complex.go b/dsp/window/window_complex.go index e705f11b..29da64c2 100644 --- a/dsp/window/window_complex.go +++ b/dsp/window/window_complex.go @@ -330,6 +330,49 @@ func (g GaussianComplex) Transform(seq []complex128) []complex128 { return seq } +// TukeyComplex can modify a sequence using the Tukey window and return the result. +// See https://en.wikipedia.org/wiki/Window_function#Tukey_window +// and https://www.recordingblogs.com/wiki/tukey-window for details. +// +// The Tukey window is an adjustable window. +// +// The sequence weights are +// w[k] = 0.5 * (1 + cos(π*(|k - M| - αM)/((1-α) * M))), |k - M| ≥ αM +// = 1, |k - M| < αM +// with M = (N - 1)/2 for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters are summarized in the table: +// | α=0.3 | α=0.5 | α=0.7 | +// -------|--------------------------| +// ΔF_0 | 1.33 | 1.22 | 1.13 | +// ΔF_0.5 | 1.28 | 1.16 | 1.04 | +// K | 0.67 | 0.61 | 0.57 | +// ɣ_max | -18.2 | -15.1 | -13.8 | +// β | -1.41 | -2.50 | -3.74 | +type TukeyComplex struct { + Alpha float64 +} + +// Transform applies the Tukey transformation to seq in place, using the value +// of the receiver as the Alpha parameter, and returning the result +func (t TukeyComplex) Transform(seq []complex128) []complex128 { + switch { + case t.Alpha <= 0: + return RectangularComplex(seq) + case t.Alpha >= 1: + return HannComplex(seq) + default: + alphaL := t.Alpha * float64(len(seq)-1) + width := int(0.5*alphaL) + 1 + for i := range seq[:width] { + w := complex(0.5*(1-math.Cos(2*math.Pi*float64(i)/alphaL)), 0) + seq[i] *= w + seq[len(seq)-1-i] *= w + } + return seq + } +} + // ValuesComplex is an arbitrary complex window function. type ValuesComplex []complex128 diff --git a/dsp/window/window_test.go b/dsp/window/window_test.go index 10d4c56a..5b899b21 100644 --- a/dsp/window/window_test.go +++ b/dsp/window/window_test.go @@ -5,7 +5,6 @@ package window import ( - "fmt" "testing" "gonum.org/v1/gonum/floats" @@ -111,34 +110,48 @@ var windowTests = []struct { 0.967756, 0.737755, 0.402051, 0.115529, -0.037444, -0.070137, -0.045939, -0.017675, -0.003687, -0.000421, }, }, -} - -var gausWindowTests = []struct { - name string - sigma float64 - want []float64 -}{ { - name: "Gaussian", sigma: 0.3, + name: "Gaussian_0.3", fn: Gaussian{0.3}.Transform, fnCmplx: GaussianComplex{0.3}.Transform, want: []float64{ 0.003866, 0.011708, 0.031348, 0.074214, 0.155344, 0.287499, 0.470444, 0.680632, 0.870660, 0.984728, 0.984728, 0.870660, 0.680632, 0.470444, 0.287499, 0.155344, 0.074214, 0.031348, 0.011708, 0.003866, }, }, { - name: "Gaussian", sigma: 0.5, + name: "Gaussian_0.5", fn: Gaussian{0.5}.Transform, fnCmplx: GaussianComplex{0.5}.Transform, want: []float64{ 0.135335, 0.201673, 0.287499, 0.392081, 0.511524, 0.638423, 0.762260, 0.870660, 0.951361, 0.994475, 0.994475, 0.951361, 0.870660, 0.762260, 0.638423, 0.511524, 0.392081, 0.287499, 0.201673, 0.135335, }, }, { - name: "Gaussian", sigma: 1.2, + name: "Gaussian_1.2", fn: Gaussian{1.2}.Transform, fnCmplx: GaussianComplex{1.2}.Transform, want: []float64{ 0.706648, 0.757319, 0.805403, 0.849974, 0.890135, 0.925049, 0.953963, 0.976241, 0.991381, 0.999039, 0.999039, 0.991381, 0.976241, 0.953963, 0.925049, 0.890135, 0.849974, 0.805403, 0.757319, 0.706648, }, }, + { + name: "Tukey_1", fn: Tukey{1}.Transform, fnCmplx: TukeyComplex{1}.Transform, + want: []float64{ // Hann window case. + 0.000000, 0.027091, 0.105430, 0.226526, 0.377257, 0.541290, 0.700848, 0.838641, 0.939737, 0.993181, + 0.993181, 0.939737, 0.838641, 0.700848, 0.541290, 0.377257, 0.226526, 0.105430, 0.027091, 0.000000, + }, + }, + { + name: "Tukey_0", fn: Tukey{0}.Transform, fnCmplx: TukeyComplex{0}.Transform, + want: []float64{ // Rectangular window case. + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + }, + }, + { + name: "Tukey_0.5", fn: Tukey{0.5}.Transform, fnCmplx: TukeyComplex{0.5}.Transform, + want: []float64{ + 0.000000, 0.105430, 0.377257, 0.700847, 0.939737, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, + 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.939737, 0.700847, 0.377257, 0.105429, 0.000000, + }, + }, } func TestWindows(t *testing.T) { @@ -169,36 +182,6 @@ func TestWindows(t *testing.T) { } } -func TestGausWindows(t *testing.T) { - t.Parallel() - const tol = 1e-6 - - for _, test := range gausWindowTests { - t.Run(fmt.Sprintf("%s (sigma=%.1f)", test.name, test.sigma), func(t *testing.T) { - src := make([]float64, 20) - for i := range src { - src[i] = 1 - } - - gaussian := Gaussian{test.sigma} - - dst := gaussian.Transform(src) - if !floats.EqualApprox(dst, test.want, tol) { - t.Errorf("unexpected result for window function %q:\ngot:%#.6v\nwant:%#v", test.name, dst, test.want) - } - - for i := range src { - src[i] = 1 - } - - dst = NewValues(gaussian.Transform, len(src)).Transform(src) - if !floats.EqualApprox(dst, test.want, tol) { - t.Errorf("unexpected result for lookup window function %q:\ngot:%#.6v\nwant:%#.6v", test.name, dst, test.want) - } - }) - } -} - func TestWindowsComplex(t *testing.T) { t.Parallel() const tol = 1e-6 @@ -227,36 +210,6 @@ func TestWindowsComplex(t *testing.T) { } } -func TestGausWindowComplex(t *testing.T) { - t.Parallel() - const tol = 1e-6 - - for _, test := range gausWindowTests { - t.Run(fmt.Sprintf("%sComplex (sigma=%.1f)", test.name, test.sigma), func(t *testing.T) { - src := make([]complex128, 20) - for i := range src { - src[i] = complex(1, 1) - } - - gaussian := GaussianComplex{test.sigma} - - dst := gaussian.Transform(src) - if !equalApprox(dst, test.want, tol) { - t.Errorf("unexpected result for window function %q:\ngot:%#.6v\nwant:%#.6v", test.name, dst, test.want) - } - - for i := range src { - src[i] = complex(1, 1) - } - - dst = NewValuesComplex(gaussian.Transform, len(src)).Transform(src) - if !equalApprox(dst, test.want, tol) { - t.Errorf("unexpected result for lookup window function %q:\ngot:%#.6v\nwant:%#.6v", test.name, dst, test.want) - } - }) - } -} - func equalApprox(seq1 []complex128, seq2 []float64, tol float64) bool { if len(seq1) != len(seq2) { return false