package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-06 23:52:47 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

dsp/window: refactor Gaussian and Tukey windows to be each a single type

Browse Source
This commit is contained in:
Dan Kortschak
2020-09-16 06:01:44 +09:30
parent 880b710e3c
commit 8c25d69b14
4 changed files with 181 additions and 222 deletions
Download Patch File Download Diff File Expand all files Collapse all files

108
dsp/window/window.go
View File

@@ -293,111 +293,3 @@ func FlatTop(seq []float64) []float64 {
} }
return seq 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
}

108
dsp/window/window_complex.go
View File

@@ -293,111 +293,3 @@ func FlatTopComplex(seq []complex128) []complex128 {
} }
return seq 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
}

175
dsp/window/window_parametric.go Normal file
View File

@@ -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
}

12
dsp/window/window_test.go
View File

@@ -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{ want: []float64{
0.003866, 0.011708, 0.031348, 0.074214, 0.155344, 0.287499, 0.470444, 0.680632, 0.870660, 0.984728, 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, 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{ want: []float64{
0.135335, 0.201673, 0.287499, 0.392081, 0.511524, 0.638423, 0.762260, 0.870660, 0.951361, 0.994475, 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, 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{ want: []float64{
0.706648, 0.757319, 0.805403, 0.849974, 0.890135, 0.925049, 0.953963, 0.976241, 0.991381, 0.999039, 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, 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. 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.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, 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. 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,
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{ want: []float64{
0.000000, 0.105430, 0.377257, 0.700847, 0.939737, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 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, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.939737, 0.700847, 0.377257, 0.105429, 0.000000,