diff --git a/dsp/window/window.go b/dsp/window/window.go index da6366e1..e6f4bf14 100644 --- a/dsp/window/window.go +++ b/dsp/window/window.go @@ -293,111 +293,3 @@ func FlatTop(seq []float64) []float64 { } return seq } - -// Gaussian can modify a sequence by 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. -// -// The Gaussian window is an adjustable window. -// -// The sequence weights are -// w[k] = exp(-0.5 * ((k-M)/(σ*M))² ), M = (N-1)/2, -// for k=0,1,...,N-1 where N is the length of the window. -// -// The properties of the window depend on the value of σ (sigma). -// It can be used as high or low resolution window, depending of the σ value. -// -// Spectral leakage parameters are summarized in the table: -// | σ=0.3 | σ=0.5 | σ=1.2 | -// -------|---------------------------| -// ΔF_0 | 8 | 3.4 | 2.2 | -// ΔF_0.5 | 1.82 | 1.2 | 0.94 | -// K | 4 | 1.7 | 1.1 | -// ɣ_max | -65 | -31.5 | -15.5 | -// β | -8.52 | -4.48 | -0.96 | -type Gaussian struct { - Sigma float64 -} - -// Transform applies the Gaussian transformation to seq in place, using the value -// of the receiver as the sigma parameter, and returning the result. -func (g Gaussian) Transform(seq []float64) []float64 { - a := float64(len(seq)-1) / 2 - for i := range seq { - x := -0.5 * math.Pow((float64(i)-a)/(g.Sigma*a), 2) - seq[i] *= math.Exp(x) - } - 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 - -// NewValues returns a Values of length n with weights corresponding to the -// provided window function. -func NewValues(window func([]float64) []float64, n int) Values { - v := make(Values, n) - for i := range v { - v[i] = 1 - } - return window(v) -} - -// Transform applies the weights in the receiver to seq in place, returning the -// result. If v is nil, Transform is a no-op, otherwise the length of v must -// match the length of seq. -func (v Values) Transform(seq []float64) []float64 { - if v == nil { - return seq - } - if len(v) != len(seq) { - panic("window: length mismatch") - } - for i, w := range v { - seq[i] *= w - } - return seq -} diff --git a/dsp/window/window_complex.go b/dsp/window/window_complex.go index 29da64c2..19ede438 100644 --- a/dsp/window/window_complex.go +++ b/dsp/window/window_complex.go @@ -293,111 +293,3 @@ func FlatTopComplex(seq []complex128) []complex128 { } return seq } - -// GaussianComplex can modify a sequence by 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. -// -// The Gaussian window is an adjustable window. -// -// The sequence weights are -// w[k] = exp(-0.5 * ((k-M)/(σ*M))² ), M = (N-1)/2, -// for k=0,1,...,N-1 where N is the length of the window. -// -// The properties of the window depend on the value of σ (sigma). -// It can be used as high or low resolution window, depending of the σ value. -// -// Spectral leakage parameters are summarized in the table: -// | σ=0.3 | σ=0.5 | σ=1.2 | -// -------|---------------------------| -// ΔF_0 | 8 | 3.4 | 2.2 | -// ΔF_0.5 | 1.82 | 1.2 | 0.94 | -// K | 4 | 1.7 | 1.1 | -// ɣ_max | -65 | -31.5 | -15.5 | -// β | -8.52 | -4.48 | -0.96 | -type GaussianComplex struct { - Sigma float64 -} - -// Transform applies the Gaussian transformation to seq in place, using the value -// of the receiver as the sigma parameter, and returning the result. -func (g GaussianComplex) Transform(seq []complex128) []complex128 { - a := float64(len(seq)-1) / 2 - for i := range seq { - x := -0.5 * math.Pow((float64(i)-a)/(g.Sigma*a), 2) - seq[i] *= complex(math.Exp(x), 0) - } - 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 - -// NewValuesComplex returns a ValuesComplex of length n with weights corresponding -// to the provided window function. -func NewValuesComplex(window func([]complex128) []complex128, n int) ValuesComplex { - v := make(ValuesComplex, n) - for i := range v { - v[i] = 1 - } - return window(v) -} - -// Transform applies the weights in the receiver to seq in place, returning the -// result. If v is nil, Transform is a no-op, otherwise the length of v must -// match the length of seq. -func (v ValuesComplex) Transform(seq []complex128) []complex128 { - if v == nil { - return seq - } - if len(v) != len(seq) { - panic("window: length mismatch") - } - for i, w := range v { - seq[i] *= w - } - return seq -} diff --git a/dsp/window/window_parametric.go b/dsp/window/window_parametric.go new file mode 100644 index 00000000..5ad1d909 --- /dev/null +++ b/dsp/window/window_parametric.go @@ -0,0 +1,175 @@ +// Copyright ©2020 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 window + +import "math" + +// 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. +// +// The Gaussian window is an adjustable window. +// +// The sequence weights are +// w[k] = exp(-0.5 * ((k-M)/(σ*M))² ), M = (N-1)/2, +// for k=0,1,...,N-1 where N is the length of the window. +// +// The properties of the window depend on the value of σ (sigma). +// It can be used as high or low resolution window, depending of the σ value. +// +// Spectral leakage parameters are summarized in the table: +// | σ=0.3 | σ=0.5 | σ=1.2 | +// -------|---------------------------| +// ΔF_0 | 8 | 3.4 | 2.2 | +// ΔF_0.5 | 1.82 | 1.2 | 0.94 | +// K | 4 | 1.7 | 1.1 | +// ɣ_max | -65 | -31.5 | -15.5 | +// β | -8.52 | -4.48 | -0.96 | +type Gaussian struct { + Sigma float64 +} + +// Transform applies the Gaussian transformation to seq in place, using the value +// of the receiver as the sigma parameter, and returning the result. +func (g Gaussian) Transform(seq []float64) []float64 { + a := float64(len(seq)-1) / 2 + for i := range seq { + x := -0.5 * math.Pow((float64(i)-a)/(g.Sigma*a), 2) + seq[i] *= math.Exp(x) + } + return seq +} + +// TransformComplex applies the Gaussian transformation to seq in place, using the value +// of the receiver as the sigma parameter, and returning the result. +func (g Gaussian) TransformComplex(seq []complex128) []complex128 { + a := float64(len(seq)-1) / 2 + for i := range seq { + x := -0.5 * math.Pow((float64(i)-a)/(g.Sigma*a), 2) + seq[i] *= complex(math.Exp(x), 0) + } + 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 + } +} + +// TransformComplex 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) TransformComplex(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 + } +} + +// Values is an arbitrary real window function. +type Values []float64 + +// NewValues returns a Values of length n with weights corresponding to the +// provided window function. +func NewValues(window func([]float64) []float64, n int) Values { + v := make(Values, n) + for i := range v { + v[i] = 1 + } + return window(v) +} + +// Transform applies the weights in the receiver to seq in place, returning the +// result. If v is nil, Transform is a no-op, otherwise the length of v must +// match the length of seq. +func (v Values) Transform(seq []float64) []float64 { + if v == nil { + return seq + } + if len(v) != len(seq) { + panic("window: length mismatch") + } + for i, w := range v { + seq[i] *= w + } + return seq +} + +// ValuesComplex is an arbitrary complex window function. +type ValuesComplex []complex128 + +// NewValuesComplex returns a ValuesComplex of length n with weights corresponding +// to the provided window function. +func NewValuesComplex(window func([]complex128) []complex128, n int) ValuesComplex { + v := make(ValuesComplex, n) + for i := range v { + v[i] = 1 + } + return window(v) +} + +// Transform applies the weights in the receiver to seq in place, returning the +// result. If v is nil, Transform is a no-op, otherwise the length of v must +// match the length of seq. +func (v ValuesComplex) Transform(seq []complex128) []complex128 { + if v == nil { + return seq + } + if len(v) != len(seq) { + panic("window: length mismatch") + } + for i, w := range v { + seq[i] *= w + } + return seq +} diff --git a/dsp/window/window_test.go b/dsp/window/window_test.go index 5b899b21..1bda1887 100644 --- a/dsp/window/window_test.go +++ b/dsp/window/window_test.go @@ -111,42 +111,42 @@ var windowTests = []struct { }, }, { - name: "Gaussian_0.3", fn: Gaussian{0.3}.Transform, fnCmplx: GaussianComplex{0.3}.Transform, + name: "Gaussian{0.3}.Transform", fn: Gaussian{0.3}.Transform, fnCmplx: Gaussian{0.3}.TransformComplex, 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_0.5", fn: Gaussian{0.5}.Transform, fnCmplx: GaussianComplex{0.5}.Transform, + name: "Gaussian{0.5}.Transform", fn: Gaussian{0.5}.Transform, fnCmplx: Gaussian{0.5}.TransformComplex, 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_1.2", fn: Gaussian{1.2}.Transform, fnCmplx: GaussianComplex{1.2}.Transform, + name: "Gaussian{1.2}.Transform", fn: Gaussian{1.2}.Transform, fnCmplx: Gaussian{1.2}.TransformComplex, 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, + name: "Tukey{1}.Transform", fn: Tukey{1}.Transform, fnCmplx: Tukey{1}.TransformComplex, 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, + name: "Tukey{0}.Transform", fn: Tukey{0}.Transform, fnCmplx: Tukey{0}.TransformComplex, 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, + name: "Tukey{0.5}.Transform", fn: Tukey{0.5}.Transform, fnCmplx: Tukey{0.5}.TransformComplex, 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,