mirror of
https://github.com/gonum/gonum.git
synced 2025-10-06 15:47:01 +08:00
dsp/window: new package for functions to control spectral leakage of FFT
This commit is contained in:
Andrei Blinnikov
43
dsp/window/doc.go
Normal file
43
dsp/window/doc.go
Normal file
@@ -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"
|
||||
325
dsp/window/window.go
Normal file
325
dsp/window/window.go
Normal file
@@ -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
|
||||
}
|
||||
325
dsp/window/window_complex.go
Normal file
325
dsp/window/window_complex.go
Normal file
@@ -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
|
||||
}
|
||||
167
dsp/window/window_complex_test.go
Normal file
167
dsp/window/window_complex_test.go
Normal file
@@ -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
|
||||
}
|
||||
105
dsp/window/window_example_test.go
Normal file
105
dsp/window/window_example_test.go
Normal file
@@ -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]
|
||||
}
|
||||
148
dsp/window/window_test.go
Normal file
148
dsp/window/window_test.go
Normal file
@@ -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)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user