dsp/window: add Tukey

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
 

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.
 //

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 {
+var monoParamWindowTests = []struct {
 	name    string
-	sigma float64
+	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,14 +258,15 @@ 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) {
+	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) {
@@ -255,6 +284,7 @@ func TestGausWindowComplex(t *testing.T) {
 			})
 		}
 	}
+}
 
 func equalApprox(seq1 []complex128, seq2 []float64, tol float64) bool {
 	if len(seq1) != len(seq2) {