diff --git a/dsp/window/doc.go b/dsp/window/doc.go new file mode 100644 index 00000000..97255b1c --- /dev/null +++ b/dsp/window/doc.go @@ -0,0 +1,43 @@ +// 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 provides a set of functions to perform the transformation +// of sequence by different window functions. +// +// Window functions can be used to control spectral leakage parameters +// when performing a Fourier transform on a signal of limited length. +// See https://en.wikipedia.org/wiki/Window_function for more details. +// +// Spectral leakage parameters +// +// Application of window functions to an input will result in changes +// to the frequency content of the signal in an effect called spectral +// leakage. See https://en.wikipedia.org/wiki/Spectral_leakage. +// +// The characteristic changes associated with each window function may +// be described using a set of spectral leakage parameters; β, ΔF_0, ΔF_0.5, +// K and ɣ_max. +// +// The β, attenuation, coefficient of a window is the ratio of the +// constant component of the spectrum resulting from use of the window +// compared to that produced using the rectangular window, expressed in +// a logarithmic scale. +// β_w = 20 log10(A_w / A_rect) dB +// +// The ΔF_0 parameter describes the normalized width of the main lobe of +// the frequency spectrum at zero amplitude. +// +// The ΔF_0.5 parameter describes the normalized width of the main lobe of +// the frequency spectrum at -3 dB (half maximum amplitude). +// +// The K parameter describes the relative width of the main lobe of the +// frequency spectrum produced by the window compared with the rectangular +// window. The rectangular window has the lowest ΔF_0 at a value of 2. +// K_w = ΔF_0_w/ΔF_0_rect. +// The value of K divides windows into high resolution windows (K≤3) and +// low resolution windows (K>3). +// +// The ɣ_max parameter is the maximum level of the side lobes of the +// normalized spectrum, in decibels. +package window // import "gonum.org/v1/gonum/dsp/window" diff --git a/dsp/window/window.go b/dsp/window/window.go new file mode 100644 index 00000000..fa791deb --- /dev/null +++ b/dsp/window/window.go @@ -0,0 +1,325 @@ +// 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" + +// Rectangular 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. +// +// The rectangular window has the lowest width of the main lobe and largest +// level of the side lobes. The result corresponds to a selection of +// limited length sequence of values without any modification. +// +// The sequence weights are +// w[k] = 1, +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 2, ΔF_0.5 = 0.89, K = 1, ɣ_max = -13, β = 0. +func Rectangular(seq []float64) []float64 { + return seq +} + +// Sine modifies seq in place by the Sine window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Sine_window and +// https://www.recordingblogs.com/wiki/sine-window for details. +// +// Sine window is a high-resolution window. +// +// The sequence weights are +// w[k] = sin(π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 3, ΔF_0.5 = 1.23, K = 1.5, ɣ_max = -23, β = -3.93. +func Sine(seq []float64) []float64 { + k := math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= math.Sin(k * float64(i)) + } + return seq +} + +// Lanczos modifies seq in place by the Lanczos window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Lanczos_window and +// https://www.recordingblogs.com/wiki/lanczos-window for details. +// +// The Lanczos window is a high-resolution window. +// +// The sequence weights are +// w[k] = sinc(2*k/(N-1) - 1), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 3.24, ΔF_0.5 = 1.3, K = 1.62, ɣ_max = -26.4, β = -4.6. +func Lanczos(seq []float64) []float64 { + k := 2 / float64(len(seq)-1) + for i := range seq { + x := math.Pi * (k*float64(i) - 1) + if x == 0 { + // Avoid NaN. + continue + } + seq[i] *= math.Sin(x) / x + } + return seq +} + +// Triangular modifies seq in place by the Triangular window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Triangular_window and +// https://www.recordingblogs.com/wiki/triangular-window for details. +// +// The Triangular window is a high-resolution window. +// +// The sequence weights are +// w[k] = 1 - |k/A -1|, A=(N-1)/2, +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.33, K = 2, ɣ_max = -26.5, β = -6. +func Triangular(seq []float64) []float64 { + a := float64(len(seq)-1) / 2 + for i := range seq { + seq[i] *= 1 - math.Abs(float64(i)/a-1) + } + return seq +} + +// Hann modifies seq in place by the Hann window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Hann_and_Hamming_windows +// and https://www.recordingblogs.com/wiki/hann-window for details. +// +// The Hann window is a high-resolution window. +// +// The sequence weights are +// w[k] = 0.5*(1 - cos(2*π*k/(N-1))), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.5, K = 2, ɣ_max = -31.5, β = -6. +func Hann(seq []float64) []float64 { + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= 0.5 * (1 - math.Cos(k*float64(i))) + } + return seq +} + +// BartlettHann modifies seq in place by the Bartlett-Hann window and returns result. +// See https://en.wikipedia.org/wiki/Window_function#Bartlett%E2%80%93Hann_window +// and https://www.recordingblogs.com/wiki/bartlett-hann-window for details. +// +// The Bartlett-Hann window is a high-resolution window. +// +// The sequence weights are +// w[k] = 0.62 - 0.48*|k/(N-1)-0.5| - 0.38*cos(2*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.45, K = 2, ɣ_max = -35.9, β = -6. +func BartlettHann(seq []float64) []float64 { + const ( + a0 = 0.62 + a1 = 0.48 + a2 = 0.38 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= a0 - a1*math.Abs(float64(i)/float64(len(seq)-1)-0.5) - a2*math.Cos(k*float64(i)) + } + return seq +} + +// Hamming modifies seq in place by the Hamming window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Hann_and_Hamming_windows +// and https://www.recordingblogs.com/wiki/hamming-window for details. +// +// The Hamming window is a high-resolution window. Among K=2 windows it has +// the highest ɣ_max. +// +// The sequence weights are +// w[k] = 25/46 - 21/46 * cos(2*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.33, K = 2, ɣ_max = -42, β = -5.37. +func Hamming(seq []float64) []float64 { + const ( + a0 = 25.0 / 46.0 + a1 = 1 - a0 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= a0 - a1*math.Cos(k*float64(i)) + } + return seq +} + +// Blackman modifies seq in place by the Blackman window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Blackman_window and +// https://www.recordingblogs.com/wiki/blackman-window for details. +// +// The Blackman window is a high-resolution window. +// +// The sequence weights are +// w[k] = 0.42 - 0.5*cos(2*π*k/(N-1)) + 0.08*cos(4*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 6, ΔF_0.5 = 1.7, K = 3, ɣ_max = -58, β = -7.54. +func Blackman(seq []float64) []float64 { + const ( + a0 = 0.42 + a1 = 0.5 + a2 = 0.08 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) + } + return seq +} + +// BlackmanHarris modifies seq in place by the Blackman-Harris window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Harris_window +// and https://www.recordingblogs.com/wiki/blackman-harris-window for details. +// +// The Blackman-Harris window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.35875 - 0.48829*cos(2*π*k/(N-1)) + +// 0.14128*cos(4*π*k/(N-1)) - 0.01168*cos(6*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 8, ΔF_0.5 = 1.97, K = 4, ɣ_max = -92, β = -8.91. +func BlackmanHarris(seq []float64) []float64 { + const ( + a0 = 0.35875 + a1 = 0.48829 + a2 = 0.14128 + a3 = 0.01168 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) + } + return seq +} + +// Nuttall modifies seq in place by the Nuttall window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Nuttall_window,_continuous_first_derivative +// and https://www.recordingblogs.com/wiki/nuttall-window for details. +// +// The Nuttall window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.355768 - 0.487396*cos(2*π*k/(N-1)) + 0.144232*cos(4*π*k/(N-1)) - +// 0.012604*cos(6*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 8, ΔF_0.5 = 1.98, K = 4, ɣ_max = -93, β = -9. +func Nuttall(seq []float64) []float64 { + const ( + a0 = 0.355768 + a1 = 0.487396 + a2 = 0.144232 + a3 = 0.012604 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) + } + return seq +} + +// BlackmanNuttall modifies seq in place by the Blackman-Nuttall window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Nuttall_window +// and https://www.recordingblogs.com/wiki/blackman-nuttall-window for details. +// +// The Blackman-Nuttall window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.3635819 - 0.4891775*cos(2*π*k/(N-1)) + 0.1365995*cos(4*π*k/(N-1)) - +// 0.0106411*cos(6*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 8, ΔF_0.5 = 1.94, K = 4, ɣ_max = -98, β = -8.8. +func BlackmanNuttall(seq []float64) []float64 { + const ( + a0 = 0.3635819 + a1 = 0.4891775 + a2 = 0.1365995 + a3 = 0.0106411 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) + } + return seq +} + +// FlatTop modifies seq in place by the Flat Top window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Flat_top_window and +// https://www.recordingblogs.com/wiki/flat-top-window for details. +// +// The Flat Top window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.21557895 - 0.41663158*cos(2*π*k/(N-1)) + +// 0.277263158*cos(4*π*k/(N-1)) - 0.083578947*cos(6*π*k/(N-1)) + +// 0.006947368*cos(4*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 10, ΔF_0.5 = 3.72, K = 5, ɣ_max = -93.0, β = -13.34. +func FlatTop(seq []float64) []float64 { + const ( + a0 = 0.21557895 + a1 = 0.41663158 + a2 = 0.277263158 + a3 = 0.083578947 + a4 = 0.006947368 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) + a4*math.Cos(4*x) + } + return seq +} + +// Gaussian modifies seq in place by the Gaussian window and returns 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 window depends on the σ (sigma) argument. +// 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 | +func Gaussian(seq []float64, sigma float64) []float64 { + a := float64(len(seq)-1) / 2 + for i := range seq { + x := -0.5 * math.Pow((float64(i)-a)/(sigma*a), 2) + seq[i] *= math.Exp(x) + } + return seq +} diff --git a/dsp/window/window_complex.go b/dsp/window/window_complex.go new file mode 100644 index 00000000..29e6ceda --- /dev/null +++ b/dsp/window/window_complex.go @@ -0,0 +1,325 @@ +// 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" + +// Rectangular 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. +// +// The rectangular window has the lowest width of the main lobe and largest +// level of the side lobes. The result corresponds to a selection of +// limited length sequence of values without any modification. +// +// The sequence weights are +// w[k] = 1, +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 2, ΔF_0.5 = 0.89, K = 1, ɣ_max = -13, β = 0. +func RectangularComplex(seq []complex128) []complex128 { + return seq +} + +// SineComplex modifies seq in place by the Sine window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Sine_window and +// https://www.recordingblogs.com/wiki/sine-window for details. +// +// Sine window is a high-resolution window. +// +// The sequence weights are +// w[k] = sin(π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 3, ΔF_0.5 = 1.23, K = 1.5, ɣ_max = -23, β = -3.93. +func SineComplex(seq []complex128) []complex128 { + k := math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= complex(math.Sin(k*float64(i)), 0) + } + return seq +} + +// LanczosComplex modifies seq in place by the Lanczos window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Lanczos_window and +// https://www.recordingblogs.com/wiki/lanczos-window for details. +// +// The Lanczos window is a high-resolution window. +// +// The sequence weights are +// w[k] = sinc(2*k/(N-1) - 1), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 3.24, ΔF_0.5 = 1.3, K = 1.62, ɣ_max = -26.4, β = -4.6. +func LanczosComplex(seq []complex128) []complex128 { + k := 2 / float64(len(seq)-1) + for i := range seq { + x := math.Pi * (k*float64(i) - 1) + if x == 0 { + // Avoid NaN. + continue + } + seq[i] *= complex(math.Sin(x)/x, 0) + } + return seq +} + +// TriangularComplex modifies seq in place by the Triangular window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Triangular_window and +// https://www.recordingblogs.com/wiki/triangular-window for details. +// +// The Triangular window is a high-resolution window. +// +// The sequence weights are +// w[k] = 1 - |k/A -1|, A=(N-1)/2, +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.33, K = 2, ɣ_max = -26.5, β = -6. +func TriangularComplex(seq []complex128) []complex128 { + a := float64(len(seq)-1) / 2 + for i := range seq { + seq[i] *= complex(1-math.Abs(float64(i)/a-1), 0) + } + return seq +} + +// HannComplex modifies seq in place by the Hann window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Hann_and_Hamming_windows +// and https://www.recordingblogs.com/wiki/hann-window for details. +// +// The Hann window is a high-resolution window. +// +// The sequence weights are +// w[k] = 0.5*(1 - cos(2*π*k/(N-1))), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.5, K = 2, ɣ_max = -31.5, β = -6. +func HannComplex(seq []complex128) []complex128 { + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= complex(0.5*(1-math.Cos(k*float64(i))), 0) + } + return seq +} + +// BartlettHannComplex modifies seq in place by the Bartlett-Hann window and returns result. +// See https://en.wikipedia.org/wiki/Window_function#Bartlett%E2%80%93Hann_window +// and https://www.recordingblogs.com/wiki/bartlett-hann-window for details. +// +// The Bartlett-Hann window is a high-resolution window. +// +// The sequence weights are +// w[k] = 0.62 - 0.48*|k/(N-1)-0.5| - 0.38*cos(2*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.45, K = 2, ɣ_max = -35.9, β = -6. +func BartlettHannComplex(seq []complex128) []complex128 { + const ( + a0 = 0.62 + a1 = 0.48 + a2 = 0.38 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= complex(a0-a1*math.Abs(float64(i)/float64(len(seq)-1)-0.5)-a2*math.Cos(k*float64(i)), 0) + } + return seq +} + +// HammingComplex modifies seq in place by the Hamming window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Hann_and_Hamming_windows +// and https://www.recordingblogs.com/wiki/hamming-window for details. +// +// The Hamming window is a high-resolution window. Among K=2 windows it has +// the highest ɣ_max. +// +// The sequence weights are +// w[k] = 25/46 - 21/46 * cos(2*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 4, ΔF_0.5 = 1.33, K = 2, ɣ_max = -42, β = -5.37. +func HammingComplex(seq []complex128) []complex128 { + const ( + a0 = 25.0 / 46.0 + a1 = 1 - a0 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + seq[i] *= complex(a0-a1*math.Cos(k*float64(i)), 0) + } + return seq +} + +// BlackmanComplex modifies seq in place by the Blackman window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Blackman_window and +// https://www.recordingblogs.com/wiki/blackman-window for details. +// +// The Blackman window is a high-resolution window. +// +// The sequence weights are +// w[k] = 0.42 - 0.5*cos(2*π*k/(N-1)) + 0.08*cos(4*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 6, ΔF_0.5 = 1.7, K = 3, ɣ_max = -58, β = -7.54. +func BlackmanComplex(seq []complex128) []complex128 { + const ( + a0 = 0.42 + a1 = 0.5 + a2 = 0.08 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x), 0) + } + return seq +} + +// BlackmanHarrisComplex modifies seq in place by the Blackman-Harris window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Harris_window +// and https://www.recordingblogs.com/wiki/blackman-harris-window for details. +// +// The Blackman-Harris window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.35875 - 0.48829*cos(2*π*k/(N-1)) + +// 0.14128*cos(4*π*k/(N-1)) - 0.01168*cos(6*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 8, ΔF_0.5 = 1.97, K = 4, ɣ_max = -92, β = -8.91. +func BlackmanHarrisComplex(seq []complex128) []complex128 { + const ( + a0 = 0.35875 + a1 = 0.48829 + a2 = 0.14128 + a3 = 0.01168 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x), 0) + } + return seq +} + +// NuttallComplex modifies seq in place by the Nuttall window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Nuttall_window,_continuous_first_derivative +// and https://www.recordingblogs.com/wiki/nuttall-window for details. +// +// The Nuttall window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.355768 - 0.487396*cos(2*π*k/(N-1)) + 0.144232*cos(4*π*k/(N-1)) - +// 0.012604*cos(6*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 8, ΔF_0.5 = 1.98, K = 4, ɣ_max = -93, β = -9. +func NuttallComplex(seq []complex128) []complex128 { + const ( + a0 = 0.355768 + a1 = 0.487396 + a2 = 0.144232 + a3 = 0.012604 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x), 0) + } + return seq +} + +// BlackmanNuttallComplex modifies seq in place by the Blackman-Nuttall window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Nuttall_window +// and https://www.recordingblogs.com/wiki/blackman-nuttall-window for details. +// +// The Blackman-Nuttall window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.3635819 - 0.4891775*cos(2*π*k/(N-1)) + 0.1365995*cos(4*π*k/(N-1)) - +// 0.0106411*cos(6*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 8, ΔF_0.5 = 1.94, K = 4, ɣ_max = -98, β = -8.8. +func BlackmanNuttallComplex(seq []complex128) []complex128 { + const ( + a0 = 0.3635819 + a1 = 0.4891775 + a2 = 0.1365995 + a3 = 0.0106411 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x), 0) + } + return seq +} + +// FlatTopComplex modifies seq in place by the Flat Top window and returns the result. +// See https://en.wikipedia.org/wiki/Window_function#Flat_top_window and +// https://www.recordingblogs.com/wiki/flat-top-window for details. +// +// The Flat Top window is a low-resolution window. +// +// The sequence weights are +// w[k] = 0.21557895 - 0.41663158*cos(2*π*k/(N-1)) + +// 0.277263158*cos(4*π*k/(N-1)) - 0.083578947*cos(6*π*k/(N-1)) + +// 0.006947368*cos(4*π*k/(N-1)), +// for k=0,1,...,N-1 where N is the length of the window. +// +// Spectral leakage parameters: ΔF_0 = 10, ΔF_0.5 = 3.72, K = 5, ɣ_max = -93.0, β = -13.34. +func FlatTopComplex(seq []complex128) []complex128 { + const ( + a0 = 0.21557895 + a1 = 0.41663158 + a2 = 0.277263158 + a3 = 0.083578947 + a4 = 0.006947368 + ) + + k := 2 * math.Pi / float64(len(seq)-1) + for i := range seq { + x := k * float64(i) + seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x)+a4*math.Cos(4*x), 0) + } + return seq +} + +// GaussianComplex modifies seq in place by the Gaussian window and returns 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 window depends on the σ (sigma) argument. +// 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 | +func GaussianComplex(seq []complex128, sigma float64) []complex128 { + a := float64(len(seq)-1) / 2 + for i := range seq { + x := -0.5 * math.Pow((float64(i)-a)/(sigma*a), 2) + seq[i] *= complex(math.Exp(x), 0) + } + return seq +} diff --git a/dsp/window/window_complex_test.go b/dsp/window/window_complex_test.go new file mode 100644 index 00000000..91d887e6 --- /dev/null +++ b/dsp/window/window_complex_test.go @@ -0,0 +1,167 @@ +// 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 ( + "testing" + + "gonum.org/v1/gonum/floats" +) + +// want the same value in imag part as in real part, +// so use one float64 for both +var complexWindowTests = []struct { + name string + fn func([]complex128) []complex128 + want []float64 + winLen int +}{ + { + name: "RectangularComplex", fn: RectangularComplex, winLen: 20, + want: []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, + }, + { + name: "SineComplex", fn: SineComplex, winLen: 20, + want: []float64{0.000000, 0.164595, 0.324699, 0.475947, 0.614213, 0.735724, 0.837166, 0.915773, 0.969400, 0.996584, + 0.996584, 0.969400, 0.915773, 0.837166, 0.735724, 0.614213, 0.475947, 0.324699, 0.164595, 0.000000}, + }, + { + name: "LanczosComplex", fn: LanczosComplex, winLen: 20, + want: []float64{0.000000, 0.115514, 0.247646, 0.389468, 0.532984, 0.669692, 0.791213, 0.889915, 0.959492, 0.995450, + 0.995450, 0.959492, 0.889915, 0.791213, 0.669692, 0.532984, 0.389468, 0.247646, 0.115514, 0.000000}, + }, + // This case tests Lanczos for a NaN condition. The Lanczos NaN condition is k=(N-1)/2, that is when N is odd. + { + name: "LanczosComplexOdd", fn: LanczosComplex, winLen: 21, + want: []float64{0.000000, 0.109292, 0.233872, 0.367883, 0.504551, 0.636620, 0.756827, 0.858394, 0.935489, 0.983632, + 1.000000, 0.983632, 0.935489, 0.858394, 0.756827, 0.636620, 0.504551, 0.367883, 0.233872, 0.109292, 0.000000}, + }, + { + name: "TriangularComplex", fn: TriangularComplex, winLen: 20, + want: []float64{0.000000, 0.105263, 0.210526, 0.315789, 0.421053, 0.526316, 0.631579, 0.736842, 0.842105, 0.947368, + 0.947368, 0.842105, 0.736842, 0.631579, 0.526316, 0.421053, 0.315789, 0.210526, 0.105263, 0.000000}, + }, + { + name: "HannComplex", fn: HannComplex, winLen: 20, + want: []float64{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: "BartlettHannComplex", fn: BartlettHannComplex, winLen: 20, + want: []float64{0.000000, 0.045853, 0.130653, 0.247949, 0.387768, 0.537696, 0.684223, 0.814209, 0.916305, 0.982186, + 0.982186, 0.916305, 0.814209, 0.684223, 0.537696, 0.387768, 0.247949, 0.130653, 0.045853, 0.000000}, + }, + { + name: "HammingComplex", fn: HammingComplex, winLen: 20, + want: []float64{0.086957, 0.111692, 0.183218, 0.293785, 0.431408, 0.581178, 0.726861, 0.852672, 0.944977, 0.993774, + 0.993774, 0.944977, 0.852672, 0.726861, 0.581178, 0.431409, 0.293785, 0.183218, 0.111692, 0.086957}, + }, + { + name: "BlackmanComplex", fn: BlackmanComplex, winLen: 20, + want: []float64{0.000000, 0.010223, 0.045069, 0.114390, 0.226899, 0.382381, 0.566665, 0.752034, 0.903493, 0.988846, + 0.988846, 0.903493, 0.752034, 0.566665, 0.382381, 0.226899, 0.114390, 0.045069, 0.010223, 0.000000}, + }, + { + name: "BlackmanHarrisComplex", fn: BlackmanHarrisComplex, winLen: 20, + want: []float64{0.000060, 0.002018, 0.012795, 0.046450, 0.122540, 0.256852, 0.448160, 0.668576, 0.866426, 0.984278, + 0.984278, 0.866426, 0.668576, 0.448160, 0.256852, 0.122540, 0.046450, 0.012795, 0.002018, 0.000060}, + }, + { + name: "NuttallComplex", fn: NuttallComplex, winLen: 20, + want: []float64{0.000000, 0.001706, 0.011614, 0.043682, 0.117808, 0.250658, 0.441946, 0.664015, 0.864348, 0.984019, + 0.984019, 0.864348, 0.664015, 0.441946, 0.250658, 0.117808, 0.043682, 0.011614, 0.001706, 0.000000}, + }, + { + name: "BlackmanNuttallComplex", fn: BlackmanNuttallComplex, winLen: 20, + want: []float64{0.000363, 0.002885, 0.015360, 0.051652, 0.130567, 0.266629, 0.457501, 0.675215, 0.869392, 0.984644, + 0.984644, 0.869392, 0.675215, 0.457501, 0.266629, 0.130567, 0.051652, 0.015360, 0.002885, 0.000363}, + }, + { + name: "FlatTopComplex", fn: FlatTopComplex, winLen: 20, + want: []float64{-0.000421, -0.003687, -0.017675, -0.045939, -0.070137, -0.037444, 0.115529, 0.402051, 0.737755, 0.967756, + 0.967756, 0.737755, 0.402051, 0.115529, -0.037444, -0.070137, -0.045939, -0.017675, -0.003687, -0.000421}, + }, +} + +// want the same value in imag part as in real part, +// so use one float64 for both +var complexGausWindowTests = []struct { + name string + sigma float64 + want []float64 +}{ + { + name: "GaussianComplex (sigma=0.3)", sigma: 0.3, + 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: "GaussianComplex (sigma=0.5)", sigma: 0.5, + 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: "GaussianComplex (sigma=1.2)", sigma: 1.2, + 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}, + }, +} + +func TestWindowsComplex(t *testing.T) { + const tol = 1e-6 + + for _, test := range complexWindowTests { + t.Run(test.name, func(t *testing.T) { + src := make([]complex128, test.winLen) + for i := range src { + src[i] = complex(1, 1) + } + + dst := test.fn(src) + + if !equalApprox(dst, test.want, tol) { + t.Errorf("unexpected result for window function %q:\ngot:%v\nwant:%v", test.name, dst, test.want) + } + }) + } +} + +func TestGausWindowComplex(t *testing.T) { + const tol = 1e-6 + + src := make([]complex128, 20) + for i := range src { + src[i] = complex(1, 1) + } + + for _, test := range complexGausWindowTests { + t.Run(test.name, func(t *testing.T) { + // Copy the input since we are mutating it. + srcCpy := make([]complex128, len(src)) + copy(srcCpy, src) + dst := GaussianComplex(srcCpy, test.sigma) + + if !equalApprox(dst, test.want, tol) { + t.Errorf("unexpected result for window function %q:\ngot:%v\nwant:%v", test.name, dst, test.want) + } + }) + } +} + +func equalApprox(seq1 []complex128, seq2 []float64, tol float64) bool { + if len(seq1) != len(seq2) { + return false + } + for i := range seq1 { + if !floats.EqualWithinAbsOrRel(real(seq1[i]), seq2[i], tol, tol) { + return false + } + if !floats.EqualWithinAbsOrRel(imag(seq1[i]), seq2[i], tol, tol) { + return false + } + } + return true +} diff --git a/dsp/window/window_example_test.go b/dsp/window/window_example_test.go new file mode 100644 index 00000000..1cff32ab --- /dev/null +++ b/dsp/window/window_example_test.go @@ -0,0 +1,105 @@ +// 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_test + +import ( + "fmt" + "math" + "math/cmplx" + + "gonum.org/v1/gonum/dsp/window" + "gonum.org/v1/gonum/fourier" +) + +func Example() { + // The input sequence is a 2.5 period of the Sin function. + src := make([]float64, 20) + k := 5 * math.Pi / float64(len(src)-1) + for i := range src { + src[i] = math.Sin(k * float64(i)) + } + + // Initialize an FFT and perform the analysis. + fft := fourier.NewFFT(len(src)) + coeff := fft.Coefficients(nil, src) + + // The result shows that width of the main lobe with center + // between frequencies 0.1 and 0.15 is small, but that the + // height of the side lobes is large. + fmt.Println("Rectangular window (or no window):") + for i, c := range coeff { + fmt.Printf("freq=%.4f\tcycles/period, magnitude=%.4f,\tphase=%.4f\n", + fft.Freq(i), cmplx.Abs(c), cmplx.Phase(c)) + } + + // Initialize an FFT and perform the analysis on a sequence + // transformed by the Hamming window function. + fft = fourier.NewFFT(len(src)) + coeff = fft.Coefficients(nil, window.Hamming(src)) + + // The result shows that width of the main lobe is wider, + // but height of the side lobes is lower. + fmt.Println("Hamming window:") + // The magnitude of all bins has been decreased by β. + // Generally in an analysis amplification may be omitted, but to + // make a comparable data, the result should be amplified by -β + // of the window function — +5.37 dB for the Hamming window. + // -β = 20 log_10(amplifier). + amplifier := math.Pow(10, 5.37/20.0) + for i, c := range coeff { + fmt.Printf("freq=%.4f\tcycles/period, magnitude=%.4f,\tphase=%.4f\n", + fft.Freq(i), amplifier*cmplx.Abs(c), cmplx.Phase(c)) + } + // Output: + // + // Rectangular window (or no window): + // freq=0.0000 cycles/period, magnitude=2.2798, phase=0.0000 + // freq=0.0500 cycles/period, magnitude=2.6542, phase=0.1571 + // freq=0.1000 cycles/period, magnitude=5.3115, phase=0.3142 + // freq=0.1500 cycles/period, magnitude=7.3247, phase=-2.6704 + // freq=0.2000 cycles/period, magnitude=1.6163, phase=-2.5133 + // freq=0.2500 cycles/period, magnitude=0.7681, phase=-2.3562 + // freq=0.3000 cycles/period, magnitude=0.4385, phase=-2.1991 + // freq=0.3500 cycles/period, magnitude=0.2640, phase=-2.0420 + // freq=0.4000 cycles/period, magnitude=0.1530, phase=-1.8850 + // freq=0.4500 cycles/period, magnitude=0.0707, phase=-1.7279 + // freq=0.5000 cycles/period, magnitude=0.0000, phase=0.0000 + // Hamming window: + // freq=0.0000 cycles/period, magnitude=0.0218, phase=3.1416 + // freq=0.0500 cycles/period, magnitude=0.8022, phase=-2.9845 + // freq=0.1000 cycles/period, magnitude=7.1723, phase=0.3142 + // freq=0.1500 cycles/period, magnitude=8.6285, phase=-2.6704 + // freq=0.2000 cycles/period, magnitude=2.0420, phase=0.6283 + // freq=0.2500 cycles/period, magnitude=0.0702, phase=0.7854 + // freq=0.3000 cycles/period, magnitude=0.0217, phase=-2.1991 + // freq=0.3500 cycles/period, magnitude=0.0259, phase=-2.0420 + // freq=0.4000 cycles/period, magnitude=0.0184, phase=-1.8850 + // freq=0.4500 cycles/period, magnitude=0.0092, phase=-1.7279 + // freq=0.5000 cycles/period, magnitude=0.0000, phase=0.0000 +} + +func ExampleHamming() { + src := []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} + + // Window functions change data in place. So, if input data + // needs to stay unchanged, it must be copied. + srcCpy := append([]float64(nil), src...) + // Apply window function to srcCpy. + dst := window.Hamming(srcCpy) + + // src is unchanged. + fmt.Printf("src: %f\n", src) + // srcCpy is altered. + fmt.Printf("srcCpy: %f\n", srcCpy) + // dst mirrors the srcCpy slice. + fmt.Printf("dst: %f\n", dst) + + // Output: + // + // src: [1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000] + // srcCpy: [0.086957 0.111692 0.183218 0.293785 0.431409 0.581178 0.726861 0.852672 0.944977 0.993774 0.993774 0.944977 0.852672 0.726861 0.581178 0.431409 0.293785 0.183218 0.111692 0.086957] + // dst: [0.086957 0.111692 0.183218 0.293785 0.431409 0.581178 0.726861 0.852672 0.944977 0.993774 0.993774 0.944977 0.852672 0.726861 0.581178 0.431409 0.293785 0.183218 0.111692 0.086957] +} diff --git a/dsp/window/window_test.go b/dsp/window/window_test.go new file mode 100644 index 00000000..b9bc5c94 --- /dev/null +++ b/dsp/window/window_test.go @@ -0,0 +1,148 @@ +// 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 ( + "testing" + + "gonum.org/v1/gonum/floats" +) + +var windowTests = []struct { + name string + fn func([]float64) []float64 + want []float64 + winLen int +}{ + { + name: "Rectangular", fn: Rectangular, winLen: 20, + want: []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, + }, + { + name: "Sine", fn: Sine, winLen: 20, + want: []float64{0.000000, 0.164595, 0.324699, 0.475947, 0.614213, 0.735724, 0.837166, 0.915773, 0.969400, 0.996584, + 0.996584, 0.969400, 0.915773, 0.837166, 0.735724, 0.614213, 0.475947, 0.324699, 0.164595, 0.000000}, + }, + { + name: "Lanczos", fn: Lanczos, winLen: 20, + want: []float64{0.000000, 0.115514, 0.247646, 0.389468, 0.532984, 0.669692, 0.791213, 0.889915, 0.959492, 0.995450, + 0.995450, 0.959492, 0.889915, 0.791213, 0.669692, 0.532984, 0.389468, 0.247646, 0.115514, 0.000000}, + }, + // This case tests Lanczos for a NaN condition. The Lanczos NaN condition is k=(N-1)/2, that is when N is odd. + { + name: "LanczosOdd", fn: Lanczos, winLen: 21, + want: []float64{0.000000, 0.109292, 0.233872, 0.367883, 0.504551, 0.636620, 0.756827, 0.858394, 0.935489, 0.983632, + 1.000000, 0.983632, 0.935489, 0.858394, 0.756827, 0.636620, 0.504551, 0.367883, 0.233872, 0.109292, 0.000000}, + }, + { + name: "Triangular", fn: Triangular, winLen: 20, + want: []float64{0.000000, 0.105263, 0.210526, 0.315789, 0.421053, 0.526316, 0.631579, 0.736842, 0.842105, 0.947368, + 0.947368, 0.842105, 0.736842, 0.631579, 0.526316, 0.421053, 0.315789, 0.210526, 0.105263, 0.000000}, + }, + { + name: "Hann", fn: Hann, winLen: 20, + want: []float64{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: "BartlettHann", fn: BartlettHann, winLen: 20, + want: []float64{0.000000, 0.045853, 0.130653, 0.247949, 0.387768, 0.537696, 0.684223, 0.814209, 0.916305, 0.982186, + 0.982186, 0.916305, 0.814209, 0.684223, 0.537696, 0.387768, 0.247949, 0.130653, 0.045853, 0.000000}, + }, + { + name: "Hamming", fn: Hamming, winLen: 20, + want: []float64{0.086957, 0.111692, 0.183218, 0.293785, 0.431408, 0.581178, 0.726861, 0.852672, 0.944977, 0.993774, + 0.993774, 0.944977, 0.852672, 0.726861, 0.581178, 0.431409, 0.293785, 0.183218, 0.111692, 0.086957}, + }, + { + name: "Blackman", fn: Blackman, winLen: 20, + want: []float64{0.000000, 0.010223, 0.045069, 0.114390, 0.226899, 0.382381, 0.566665, 0.752034, 0.903493, 0.988846, + 0.988846, 0.903493, 0.752034, 0.566665, 0.382381, 0.226899, 0.114390, 0.045069, 0.010223, 0.000000}, + }, + { + name: "BlackmanHarris", fn: BlackmanHarris, winLen: 20, + want: []float64{0.000060, 0.002018, 0.012795, 0.046450, 0.122540, 0.256852, 0.448160, 0.668576, 0.866426, 0.984278, + 0.984278, 0.866426, 0.668576, 0.448160, 0.256852, 0.122540, 0.046450, 0.012795, 0.002018, 0.000060}, + }, + { + name: "Nuttall", fn: Nuttall, winLen: 20, + want: []float64{0.000000, 0.001706, 0.011614, 0.043682, 0.117808, 0.250658, 0.441946, 0.664015, 0.864348, 0.984019, + 0.984019, 0.864348, 0.664015, 0.441946, 0.250658, 0.117808, 0.043682, 0.011614, 0.001706, 0.000000}, + }, + { + name: "BlackmanNuttall", fn: BlackmanNuttall, winLen: 20, + want: []float64{0.000363, 0.002885, 0.015360, 0.051652, 0.130567, 0.266629, 0.457501, 0.675215, 0.869392, 0.984644, + 0.984644, 0.869392, 0.675215, 0.457501, 0.266629, 0.130567, 0.051652, 0.015360, 0.002885, 0.000363}, + }, + { + name: "FlatTop", fn: FlatTop, winLen: 20, + want: []float64{-0.000421, -0.003687, -0.017675, -0.045939, -0.070137, -0.037444, 0.115529, 0.402051, 0.737755, 0.967756, + 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)", sigma: 0.3, + 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)", sigma: 0.5, + 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)", sigma: 1.2, + 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}, + }, +} + +func TestWindows(t *testing.T) { + const tol = 1e-6 + + for _, test := range windowTests { + t.Run(test.name, func(t *testing.T) { + src := make([]float64, test.winLen) + for i := range src { + src[i] = 1 + } + + dst := test.fn(src) + + if !floats.EqualApprox(dst, test.want, tol) { + t.Errorf("unexpected result for window function %q:\ngot:%v\nwant:%v", test.name, dst, test.want) + } + }) + } +} + +func TestGausWindows(t *testing.T) { + const tol = 1e-6 + + src := make([]float64, 20) + for i := range src { + src[i] = 1 + } + + for _, test := range gausWindowTests { + t.Run(test.name, func(t *testing.T) { + // Copy the input since we are mutating it. + srcCpy := make([]float64, len(src)) + copy(srcCpy, src) + dst := Gaussian(srcCpy, test.sigma) + + if !floats.EqualApprox(dst, test.want, tol) { + t.Errorf("unexpected result for window function %q:\ngot:%v\nwant:%v", test.name, dst, test.want) + } + }) + } +}