diff --git a/dsp/window/window.go b/dsp/window/window.go index 997c7ca7..6e14082a 100644 --- a/dsp/window/window.go +++ b/dsp/window/window.go @@ -294,7 +294,7 @@ func FlatTop(seq []float64) []float64 { return seq } -// Gaussian can modify a sequence by the Gaussian window and return the result. +// Gaussian can modify a sequence using the Gaussian window and return the result. // See https://en.wikipedia.org/wiki/Window_function#Gaussian_window // and https://www.recordingblogs.com/wiki/gaussian-window for details. // @@ -330,6 +330,53 @@ 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://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2017/174042.pdf page 88 +// +// The Tukey window is an adjustible window. It can be thought of as something +// between a rectangular and a Hann window, with a flat center and tapered edges. +// +// The properties of the window depend on the value of α (alpha). It controls +// the fraction of the window which contains a cosine taper. α = 0.5 gives a +// window whose central 50% is flat and outer quartiles are tapered. α = 1 is +// equivalent to a Hann window. α = 0 is equivalent to a rectangular window. +// 0 <= α <= 1; if α is outside the bounds, it is treated as 0 or 1. +// +// Spectral leakage parameters are summarized in the table: +// | α=0.25 | α=0.5 | α=0.75 | +// -------|---------------------------| +// ΔF_0 | 1.1 | 1.22 | 1.36 | +// ΔF_0.5 | 1.01 | 1.15 | 2.24 | +// K | 1.13 | 1.3 | 2.5 | +// ɣ_max | -14 | -15 | -19 | +// β | -1.11 | -2.5 | -4.01 | +// +// whose values are from A.D. Poularikas, +// "The Handbook of Formulas and Tables for Signal Processing" table 7.1 +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 { + if t.Alpha <= 0 { + return Rectangular(seq) + } else if t.Alpha >= 1 { + return Hann(seq) + } + + 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 9036dc89..7a9a1cc5 100644 --- a/dsp/window/window_complex.go +++ b/dsp/window/window_complex.go @@ -6,7 +6,7 @@ package window import "math" -// Rectangular modifies seq in place by the Rectangular window and returns the result. +// RectangularComplex modifies seq in place by the Rectangular window and returns the result. // See https://en.wikipedia.org/wiki/Window_function#Rectangular_window and // https://www.recordingblogs.com/wiki/rectangular-window for details. // diff --git a/dsp/window/window_test.go b/dsp/window/window_test.go index 2b5e17f0..884d0455 100644 --- a/dsp/window/window_test.go +++ b/dsp/window/window_test.go @@ -12,6 +12,11 @@ import ( "gonum.org/v1/gonum/floats/scalar" ) +const ( + gaussWin = iota + tukeyWin +) + var windowTests = []struct { name string fn func([]float64) []float64 @@ -113,31 +118,53 @@ var windowTests = []struct { }, } -var gausWindowTests = []struct { - name string - sigma float64 - want []float64 +var monoParamWindowTests = []struct { + name string + param float64 + winType int + want []float64 }{ { - name: "Gaussian", sigma: 0.3, + name: "Gaussian", param: 0.3, winType: gaussWin, want: []float64{ 0.006645, 0.018063, 0.043936, 0.095634, 0.186270, 0.324652, 0.506336, 0.706648, 0.882497, 0.986207, 0.986207, 0.882497, 0.706648, 0.506336, 0.324652, 0.186270, 0.095634, 0.043936, 0.018063, 0.006645}, }, { - name: "Gaussian", sigma: 0.5, + name: "Gaussian", param: 0.5, winType: gaussWin, want: []float64{ 0.164474, 0.235746, 0.324652, 0.429557, 0.546074, 0.666977, 0.782705, 0.882497, 0.955997, 0.995012, 0.995012, 0.955997, 0.882497, 0.782705, 0.666977, 0.546074, 0.429557, 0.324652, 0.235746, 0.164474, }, }, { - name: "Gaussian", sigma: 1.2, + name: "Gaussian", param: 1.2, winType: gaussWin, want: []float64{ 0.730981, 0.778125, 0.822578, 0.863552, 0.900293, 0.932102, 0.958357, 0.978532, 0.992218, 0.999132, 0.999132, 0.992218, 0.978532, 0.958357, 0.932102, 0.900293, 0.863552, 0.822578, 0.778125, 0.730981, }, }, + { + name: "Tukey", param: 1, winType: tukeyWin, + want: []float64{ // copied from Hann + 0.006155, 0.054496, 0.146447, 0.273005, 0.421783, 0.578217, 0.726995, 0.853553, 0.945503, 0.993844, + 0.993844, 0.945503, 0.853553, 0.726995, 0.578217, 0.421783, 0.273005, 0.146447, 0.054496, 0.006155, + }, + }, + { + name: "Tukey", param: 0, winType: tukeyWin, + want: []float64{ // copied from rectangular + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + }, + }, + { + name: "Tukey", param: 0.5, winType: tukeyWin, + 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) { @@ -168,20 +195,21 @@ func TestWindows(t *testing.T) { } } -func TestGausWindows(t *testing.T) { +func TestMonoParametricWindows(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) { + for _, test := range monoParamWindowTests { + t.Run(fmt.Sprintf("%s (param=%.1f)", test.name, test.param), func(t *testing.T) { src := make([]float64, 20) for i := range src { src[i] = 1 } - - gaussian := Gaussian{test.sigma} - - dst := gaussian.Transform(src) + type transformer interface { + Transform([]float64) []float64 + } + trans := []transformer{Gaussian{test.param}, Tukey{test.param}}[test.winType] + dst := trans.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) } @@ -190,7 +218,7 @@ func TestGausWindows(t *testing.T) { src[i] = 1 } - dst = NewValues(gaussian.Transform, len(src)).Transform(src) + dst = NewValues(trans.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) } @@ -230,29 +258,31 @@ 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) - } + for _, test := range monoParamWindowTests { + if test.winType == gaussWin { + t.Run(fmt.Sprintf("%sComplex (sigma=%.1f)", test.name, test.param), func(t *testing.T) { + src := make([]complex128, 20) + for i := range src { + src[i] = complex(1, 1) + } - gaussian := GaussianComplex{test.sigma} + gaussian := GaussianComplex{test.param} - 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) - } + 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) - } + 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) - } - }) + 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) + } + }) + } } }