diff --git a/dsp/window/doc.go b/dsp/window/doc.go index 97255b1c..5112518b 100644 --- a/dsp/window/doc.go +++ b/dsp/window/doc.go @@ -40,4 +40,9 @@ // // The ɣ_max parameter is the maximum level of the side lobes of the // normalized spectrum, in decibels. +// +// Window interval +// +// The window intervals used in the window functions are offset by a half +// to ensure that the window extends to exclude the flanking zeros. package window // import "gonum.org/v1/gonum/dsp/window" diff --git a/dsp/window/window.go b/dsp/window/window.go index fa791deb..3be482f5 100644 --- a/dsp/window/window.go +++ b/dsp/window/window.go @@ -30,14 +30,14 @@ func Rectangular(seq []float64) []float64 { // Sine window is a high-resolution window. // // The sequence weights are -// w[k] = sin(π*k/(N-1)), +// w[k] = sin(π*(k+1/2)/N), // 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) + k := math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= math.Sin(k * float64(i)) + seq[i] *= math.Sin(k * (float64(i) + 0.5)) } return seq } @@ -49,14 +49,14 @@ func Sine(seq []float64) []float64 { // The Lanczos window is a high-resolution window. // // The sequence weights are -// w[k] = sinc(2*k/(N-1) - 1), +// w[k] = sinc(2*(k+1/2)/N - 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) + k := 2 / float64(len(seq)) for i := range seq { - x := math.Pi * (k*float64(i) - 1) + x := math.Pi * (k*(float64(i)+0.5) - 1) if x == 0 { // Avoid NaN. continue @@ -73,14 +73,14 @@ func Lanczos(seq []float64) []float64 { // The Triangular window is a high-resolution window. // // The sequence weights are -// w[k] = 1 - |k/A -1|, A=(N-1)/2, +// w[k] = 1 - |(k + 1/2 - N/2)/(N/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 + a := float64(len(seq)) / 2 for i := range seq { - seq[i] *= 1 - math.Abs(float64(i)/a-1) + seq[i] *= 1 - math.Abs((float64(i)+0.5-a)/a) } return seq } @@ -92,14 +92,14 @@ func Triangular(seq []float64) []float64 { // The Hann window is a high-resolution window. // // The sequence weights are -// w[k] = 0.5*(1 - cos(2*π*k/(N-1))), +// w[k] = 0.5*(1 - cos(2*π*(k+1/2)/N)), // 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) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= 0.5 * (1 - math.Cos(k*float64(i))) + seq[i] *= 0.5 * (1 - math.Cos(k*(float64(i)+0.5))) } return seq } @@ -111,7 +111,7 @@ func Hann(seq []float64) []float64 { // 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)), +// w[k] = 0.62 - 0.48*|(k+1/2)/N-0.5| - 0.38*cos(2*π*(k+1/2)/N), // 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. @@ -122,9 +122,9 @@ func BartlettHann(seq []float64) []float64 { a2 = 0.38 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= a0 - a1*math.Abs(float64(i)/float64(len(seq)-1)-0.5) - a2*math.Cos(k*float64(i)) + seq[i] *= a0 - a1*math.Abs((float64(i)+0.5)/float64(len(seq))-0.5) - a2*math.Cos(k*(float64(i)+0.5)) } return seq } @@ -137,7 +137,7 @@ func BartlettHann(seq []float64) []float64 { // the highest ɣ_max. // // The sequence weights are -// w[k] = 25/46 - 21/46 * cos(2*π*k/(N-1)), +// w[k] = 25/46 - 21/46 * cos(2*π*(k+1/2)/N), // 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. @@ -147,9 +147,9 @@ func Hamming(seq []float64) []float64 { a1 = 1 - a0 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= a0 - a1*math.Cos(k*float64(i)) + seq[i] *= a0 - a1*math.Cos(k*(float64(i)+0.5)) } return seq } @@ -161,7 +161,7 @@ func Hamming(seq []float64) []float64 { // 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)), +// w[k] = 0.42 - 0.5*cos(2*π*(k+1/2)/N) + 0.08*cos(4*π*(k+1/2)/N), // 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. @@ -172,9 +172,9 @@ func Blackman(seq []float64) []float64 { a2 = 0.08 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) } return seq @@ -187,8 +187,8 @@ func Blackman(seq []float64) []float64 { // 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)), +// w[k] = 0.35875 - 0.48829*cos(2*π*(k+1/2)/N) + +// 0.14128*cos(4*π*(k+1/2)/N) - 0.01168*cos(6*π*(k+1/2)/N), // 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. @@ -200,9 +200,9 @@ func BlackmanHarris(seq []float64) []float64 { a3 = 0.01168 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) } return seq @@ -215,8 +215,8 @@ func BlackmanHarris(seq []float64) []float64 { // 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)), +// w[k] = 0.355768 - 0.487396*cos(2*π*(k+1/2)/N) + 0.144232*cos(4*π*(k+1/2)/N) - +// 0.012604*cos(6*π*(k+1/2)/N), // 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. @@ -228,9 +228,9 @@ func Nuttall(seq []float64) []float64 { a3 = 0.012604 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) } return seq @@ -243,8 +243,8 @@ func Nuttall(seq []float64) []float64 { // 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)), +// w[k] = 0.3635819 - 0.4891775*cos(2*π*(k+1/2)/N) + 0.1365995*cos(4*π*(k+1/2)/N) - +// 0.0106411*cos(6*π*(k+1/2)/N), // 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. @@ -256,9 +256,9 @@ func BlackmanNuttall(seq []float64) []float64 { a3 = 0.0106411 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) } return seq @@ -271,9 +271,9 @@ func BlackmanNuttall(seq []float64) []float64 { // 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)), +// w[k] = 0.21557895 - 0.41663158*cos(2*π*(k+1/2)/(N-1)) + +// 0.277263158*cos(4*π*(k+1/2)/N) - 0.083578947*cos(6*π*(k+1/2)/N) + +// 0.006947368*cos(4*π*(k+1/2)/N), // 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. @@ -286,9 +286,9 @@ func FlatTop(seq []float64) []float64 { a4 = 0.006947368 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= a0 - a1*math.Cos(x) + a2*math.Cos(2*x) - a3*math.Cos(3*x) + a4*math.Cos(4*x) } return seq @@ -301,10 +301,10 @@ func FlatTop(seq []float64) []float64 { // The Gaussian window is an adjustable window. // // The sequence weights are -// w[k] = exp(-0.5 * ((k-M)/(σ*M))² ), M = (N-1)/2, +// w[k] = exp(-0.5 * ((k + 1/2 - M)/(σ*M))² ), M = N/2, // for k=0,1,...,N-1 where N is the length of the window. // -// The properties of window depends on the σ (sigma) argument. +// The properties of the window depend 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: @@ -316,9 +316,9 @@ func FlatTop(seq []float64) []float64 { // ɣ_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 + a := float64(len(seq)) / 2 for i := range seq { - x := -0.5 * math.Pow((float64(i)-a)/(sigma*a), 2) + x := -0.5 * math.Pow(((float64(i)+0.5)-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 index 29e6ceda..39c2d5c8 100644 --- a/dsp/window/window_complex.go +++ b/dsp/window/window_complex.go @@ -30,14 +30,14 @@ func RectangularComplex(seq []complex128) []complex128 { // Sine window is a high-resolution window. // // The sequence weights are -// w[k] = sin(π*k/(N-1)), +// w[k] = sin(π*(k+1/2)/N), // 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) + k := math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= complex(math.Sin(k*float64(i)), 0) + seq[i] *= complex(math.Sin(k*(float64(i)+0.5)), 0) } return seq } @@ -49,14 +49,14 @@ func SineComplex(seq []complex128) []complex128 { // The Lanczos window is a high-resolution window. // // The sequence weights are -// w[k] = sinc(2*k/(N-1) - 1), +// w[k] = sinc(2*(k+1/2)/N - 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) + k := 2 / float64(len(seq)) for i := range seq { - x := math.Pi * (k*float64(i) - 1) + x := math.Pi * (k*(float64(i)+0.5) - 1) if x == 0 { // Avoid NaN. continue @@ -73,14 +73,14 @@ func LanczosComplex(seq []complex128) []complex128 { // The Triangular window is a high-resolution window. // // The sequence weights are -// w[k] = 1 - |k/A -1|, A=(N-1)/2, +// w[k] = 1 - |(k + 1/2 - N/2)/(N/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 + a := float64(len(seq)) / 2 for i := range seq { - seq[i] *= complex(1-math.Abs(float64(i)/a-1), 0) + seq[i] *= complex(1-math.Abs((float64(i)+0.5-a)/a), 0) } return seq } @@ -92,14 +92,14 @@ func TriangularComplex(seq []complex128) []complex128 { // The Hann window is a high-resolution window. // // The sequence weights are -// w[k] = 0.5*(1 - cos(2*π*k/(N-1))), +// w[k] = 0.5*(1 - cos(2*π*(k+1/2)/N)), // 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) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= complex(0.5*(1-math.Cos(k*float64(i))), 0) + seq[i] *= complex(0.5*(1-math.Cos(k*(float64(i)+0.5))), 0) } return seq } @@ -111,7 +111,7 @@ func HannComplex(seq []complex128) []complex128 { // 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)), +// w[k] = 0.62 - 0.48*|(k+1/2)/N-0.5| - 0.38*cos(2*π*(k+1/2)/N), // 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. @@ -122,9 +122,9 @@ func BartlettHannComplex(seq []complex128) []complex128 { a2 = 0.38 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) 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) + seq[i] *= complex(a0-a1*math.Abs((float64(i)+0.5)/float64(len(seq))-0.5)-a2*math.Cos(k*(float64(i)+0.5)), 0) } return seq } @@ -137,7 +137,7 @@ func BartlettHannComplex(seq []complex128) []complex128 { // the highest ɣ_max. // // The sequence weights are -// w[k] = 25/46 - 21/46 * cos(2*π*k/(N-1)), +// w[k] = 25/46 - 21/46 * cos(2*π*(k+1/2)/N), // 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. @@ -147,9 +147,9 @@ func HammingComplex(seq []complex128) []complex128 { a1 = 1 - a0 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - seq[i] *= complex(a0-a1*math.Cos(k*float64(i)), 0) + seq[i] *= complex(a0-a1*math.Cos(k*(float64(i)+0.5)), 0) } return seq } @@ -161,7 +161,7 @@ func HammingComplex(seq []complex128) []complex128 { // 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)), +// w[k] = 0.42 - 0.5*cos(2*π*(k+1/2)/N) + 0.08*cos(4*π*(k+1/2)/N), // 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. @@ -172,9 +172,9 @@ func BlackmanComplex(seq []complex128) []complex128 { a2 = 0.08 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x), 0) } return seq @@ -187,8 +187,8 @@ func BlackmanComplex(seq []complex128) []complex128 { // 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)), +// w[k] = 0.35875 - 0.48829*cos(2*π*(k+1/2)/N) + +// 0.14128*cos(4*π*(k+1/2)/N) - 0.01168*cos(6*π*(k+1/2)/N), // 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. @@ -200,9 +200,9 @@ func BlackmanHarrisComplex(seq []complex128) []complex128 { a3 = 0.01168 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x), 0) } return seq @@ -215,8 +215,8 @@ func BlackmanHarrisComplex(seq []complex128) []complex128 { // 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)), +// w[k] = 0.355768 - 0.487396*cos(2*π*(k+1/2)/N) + 0.144232*cos(4*π*(k+1/2)/N) - +// 0.012604*cos(6*π*(k+1/2)/N), // 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. @@ -228,9 +228,9 @@ func NuttallComplex(seq []complex128) []complex128 { a3 = 0.012604 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x), 0) } return seq @@ -243,8 +243,8 @@ func NuttallComplex(seq []complex128) []complex128 { // 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)), +// w[k] = 0.3635819 - 0.4891775*cos(2*π*(k+1/2)/N) + 0.1365995*cos(4*π*(k+1/2)/N) - +// 0.0106411*cos(6*π*(k+1/2)/N), // 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. @@ -256,9 +256,9 @@ func BlackmanNuttallComplex(seq []complex128) []complex128 { a3 = 0.0106411 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) seq[i] *= complex(a0-a1*math.Cos(x)+a2*math.Cos(2*x)-a3*math.Cos(3*x), 0) } return seq @@ -271,9 +271,9 @@ func BlackmanNuttallComplex(seq []complex128) []complex128 { // 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)), +// w[k] = 0.21557895 - 0.41663158*cos(2*π*(k+1/2)/(N-1)) + +// 0.277263158*cos(4*π*(k+1/2)/N) - 0.083578947*cos(6*π*(k+1/2)/N) + +// 0.006947368*cos(4*π*(k+1/2)/N), // 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. @@ -286,9 +286,9 @@ func FlatTopComplex(seq []complex128) []complex128 { a4 = 0.006947368 ) - k := 2 * math.Pi / float64(len(seq)-1) + k := 2 * math.Pi / float64(len(seq)) for i := range seq { - x := k * float64(i) + x := k * (float64(i) + 0.5) 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 @@ -301,10 +301,10 @@ func FlatTopComplex(seq []complex128) []complex128 { // The Gaussian window is an adjustable window. // // The sequence weights are -// w[k] = exp(-0.5 * ((k-M)/(σ*M))² ), M = (N-1)/2, +// w[k] = exp(-0.5 * ((k + 1/2 - M)/(σ*M))² ), M = N/2, // for k=0,1,...,N-1 where N is the length of the window. // -// The properties of window depends on the σ (sigma) argument. +// The properties of the window depend 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: @@ -316,9 +316,9 @@ func FlatTopComplex(seq []complex128) []complex128 { // ɣ_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 + a := float64(len(seq)) / 2 for i := range seq { - x := -0.5 * math.Pow((float64(i)-a)/(sigma*a), 2) + x := -0.5 * math.Pow(((float64(i)+0.5)-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 deleted file mode 100644 index 91d887e6..00000000 --- a/dsp/window/window_complex_test.go +++ /dev/null @@ -1,167 +0,0 @@ -// 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 index 9918ae2c..842b9f61 100644 --- a/dsp/window/window_example_test.go +++ b/dsp/window/window_example_test.go @@ -67,16 +67,16 @@ func Example() { // 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.0000 cycles/period, magnitude=0.0506, phase=0.0000 + // freq=0.0500 cycles/period, magnitude=0.5386, phase=-2.9845 + // freq=0.1000 cycles/period, magnitude=7.3350, phase=0.3142 + // freq=0.1500 cycles/period, magnitude=8.9523, phase=-2.6704 + // freq=0.2000 cycles/period, magnitude=1.7979, phase=0.6283 + // freq=0.2500 cycles/period, magnitude=0.0957, phase=0.7854 + // freq=0.3000 cycles/period, magnitude=0.0050, phase=-2.1991 + // freq=0.3500 cycles/period, magnitude=0.0158, phase=-2.0420 + // freq=0.4000 cycles/period, magnitude=0.0125, phase=-1.8850 + // freq=0.4500 cycles/period, magnitude=0.0065, phase=-1.7279 // freq=0.5000 cycles/period, magnitude=0.0000, phase=0.0000 } @@ -100,6 +100,6 @@ func ExampleHamming() { // 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] + // srcCpy: [0.092577 0.136714 0.220669 0.336222 0.472063 0.614894 0.750735 0.866288 0.950242 0.994379 0.994379 0.950242 0.866288 0.750735 0.614894 0.472063 0.336222 0.220669 0.136714 0.092577] + // dst: [0.092577 0.136714 0.220669 0.336222 0.472063 0.614894 0.750735 0.866288 0.950242 0.994379 0.994379 0.950242 0.866288 0.750735 0.614894 0.472063 0.336222 0.220669 0.136714 0.092577] } diff --git a/dsp/window/window_test.go b/dsp/window/window_test.go index b9bc5c94..4c5de149 100644 --- a/dsp/window/window_test.go +++ b/dsp/window/window_test.go @@ -5,82 +5,110 @@ package window import ( + "fmt" "testing" "gonum.org/v1/gonum/floats" ) var windowTests = []struct { - name string - fn func([]float64) []float64 - want []float64 - winLen int + name string + fn func([]float64) []float64 + fnCmplx func([]complex128) []complex128 + want []float64 }{ { - 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: "Rectangular", fn: Rectangular, fnCmplx: RectangularComplex, + 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: "Sine", fn: Sine, fnCmplx: SineComplex, + want: []float64{ + 0.078459, 0.233445, 0.382683, 0.522499, 0.649448, 0.760406, 0.852640, 0.923880, 0.972370, 0.996917, + 0.996917, 0.972370, 0.923880, 0.852640, 0.760406, 0.649448, 0.522499, 0.382683, 0.233445, 0.078459, + }, }, { - 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}, + name: "Lanczos", fn: Lanczos, fnCmplx: LanczosComplex, + want: []float64{ + 0.052415, 0.170011, 0.300105, 0.436333, 0.57162, 0.698647, 0.810332, 0.900316, 0.963398, 0.995893, + 0.995893, 0.963398, 0.900316, 0.810332, 0.698647, 0.57162, 0.436333, 0.300105, 0.170011, 0.052415, + }, }, // 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: "LanczosOdd", fn: Lanczos, fnCmplx: LanczosComplex, + want: []float64{ + 0.049813, 0.161128, 0.284164, 0.413497, 0.543076, 0.666582, 0.777804, 0.871026, 0.941379, 0.985147, + 1, + 0.985147, 0.941379, 0.871026, 0.777804, 0.666582, 0.543076, 0.413497, 0.284164, 0.161128, 0.049813, + }, }, { - 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: "Triangular", fn: Triangular, fnCmplx: TriangularComplex, + want: []float64{ + 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, + 0.95, 0.85, 0.75, 0.65, 0.55, 0.45, 0.35, 0.25, 0.15, 0.05, + }, }, { - 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: "Hann", fn: Hann, fnCmplx: HannComplex, + want: []float64{ + 0.006155, 0.054496, 0.146447, 0.273005, 0.421783, 0.578217, 0.726995, 0.853553, 0.945503, 0.993844, + 0.993844, 0.945503, 0.853553, 0.726995, 0.578217, 0.421783, 0.273005, 0.146447, 0.054496, 0.006155, + }, }, { - name: "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: "BartlettHann", fn: BartlettHann, fnCmplx: BartlettHannComplex, + want: []float64{ + 0.016678, 0.077417, 0.171299, 0.291484, 0.428555, 0.571445, 0.708516, 0.828701, 0.922582, 0.983322, + 0.983322, 0.922582, 0.828701, 0.708516, 0.571445, 0.428555, 0.291484, 0.171299, 0.077417, 0.016678, + }, }, { - 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: "Hamming", fn: Hamming, fnCmplx: HammingComplex, + want: []float64{ + 0.092577, 0.136714, 0.220669, 0.336222, 0.472063, 0.614894, 0.750735, 0.866288, 0.950242, 0.994379, + 0.994379, 0.950242, 0.866288, 0.750735, 0.614894, 0.472063, 0.336222, 0.220669, 0.136714, 0.092577, + }, }, { - 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: "Blackman", fn: Blackman, fnCmplx: BlackmanComplex, + want: []float64{ + 0.002240, 0.021519, 0.066446, 0.145982, 0.265698, 0.422133, 0.599972, 0.773553, 0.912526, 0.989929, + 0.989929, 0.912526, 0.773553, 0.599972, 0.422133, 0.265698, 0.145982, 0.066446, 0.021519, 0.002240, + }, }, { - 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: "BlackmanHarris", fn: BlackmanHarris, fnCmplx: BlackmanHarrisComplex, + want: []float64{ + 0.000429, 0.004895, 0.021735, 0.065564, 0.153302, 0.295468, 0.485851, 0.695764, 0.878689, 0.985801, + 0.985801, 0.878689, 0.695764, 0.485851, 0.295468, 0.153302, 0.065564, 0.021735, 0.004895, 0.000429, + }, }, { - 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: "Nuttall", fn: Nuttall, fnCmplx: NuttallComplex, + want: []float64{ + 0.000315, 0.004300, 0.020039, 0.062166, 0.148072, 0.289119, 0.479815, 0.691497, 0.876790, 0.985566, + 0.985566, 0.876790, 0.691497, 0.479815, 0.289119, 0.148072, 0.062166, 0.020039, 0.004300, 0.000315, + }, }, { - 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: "BlackmanNuttall", fn: BlackmanNuttall, fnCmplx: BlackmanNuttallComplex, + want: []float64{ + 0.000859, 0.006348, 0.025205, 0.071718, 0.161975, 0.305361, 0.494863, 0.701958, 0.881398, 0.986132, + 0.986132, 0.881398, 0.701958, 0.494863, 0.305361, 0.161975, 0.071718, 0.025205, 0.006348, 0.000859, + }, }, { - 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}, + name: "FlatTop", fn: FlatTop, fnCmplx: FlatTopComplex, + want: []float64{ + -0.001079, -0.007892, -0.026872, -0.056135, -0.069724, -0.015262, 0.157058, 0.444135, 0.760699, 0.970864, + 0.970864, 0.760699, 0.444135, 0.157058, -0.015262, -0.069724, -0.056135, -0.026872, -0.007892, -0.001079, + }, }, } @@ -90,19 +118,24 @@ var gausWindowTests = []struct { 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.3, + want: []float64{ + 0.006645, 0.018063, 0.043936, 0.095634, 0.186270, 0.324652, 0.506336, 0.706648, 0.882497, 0.986207, + 0.986207, 0.882497, 0.706648, 0.506336, 0.324652, 0.186270, 0.095634, 0.043936, 0.018063, 0.006645}, }, { - name: "Gaussian (sigma=0.5)", 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: 0.5, + want: []float64{ + 0.164474, 0.235746, 0.324652, 0.429557, 0.546074, 0.666977, 0.782705, 0.882497, 0.955997, 0.995012, + 0.995012, 0.955997, 0.882497, 0.782705, 0.666977, 0.546074, 0.429557, 0.324652, 0.235746, 0.164474, + }, }, { - name: "Gaussian (sigma=1.2)", 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}, + name: "Gaussian", sigma: 1.2, + want: []float64{ + 0.730981, 0.778125, 0.822578, 0.863552, 0.900293, 0.932102, 0.958357, 0.978532, 0.992218, 0.999132, + 0.999132, 0.992218, 0.978532, 0.958357, 0.932102, 0.900293, 0.863552, 0.822578, 0.778125, 0.730981, + }, }, } @@ -111,7 +144,7 @@ func TestWindows(t *testing.T) { for _, test := range windowTests { t.Run(test.name, func(t *testing.T) { - src := make([]float64, test.winLen) + src := make([]float64, len(test.want)) for i := range src { src[i] = 1 } @@ -119,7 +152,7 @@ func TestWindows(t *testing.T) { 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) + t.Errorf("unexpected result for window function %q:\ngot:%#.6v\nwant:%#v", test.name, dst, test.want) } }) } @@ -134,15 +167,69 @@ func TestGausWindows(t *testing.T) { } for _, test := range gausWindowTests { - t.Run(test.name, func(t *testing.T) { - // Copy the input since we are mutating it. + t.Run(fmt.Sprintf("%s (sigma=%.1f)", test.name, test.sigma), func(t *testing.T) { 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) + t.Errorf("unexpected result for window function %q:\ngot:%#.6v\nwant:%#v", test.name, dst, test.want) } }) } } + +func TestWindowsComplex(t *testing.T) { + const tol = 1e-6 + + for _, test := range windowTests { + t.Run(test.name+"Complex", func(t *testing.T) { + src := make([]complex128, len(test.want)) + for i := range src { + src[i] = complex(1, 1) + } + + dst := test.fnCmplx(src) + + if !equalApprox(dst, test.want, tol) { + t.Errorf("unexpected result for window function %q:\ngot:%#.6v\nwant:%#.6v", 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 gausWindowTests { + t.Run(fmt.Sprintf("%sComplex (sigma=%.1f)", test.name, test.sigma), func(t *testing.T) { + 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:%#.6v\nwant:%#.6v", 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 +}