gonum.org/v1/gonum@v0.14.0/dsp/window/window_parametric.go (about)

     1  // Copyright ©2020 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package window
     6  
     7  import "math"
     8  
     9  // Gaussian can modify a sequence using the Gaussian window and return the
    10  // result.
    11  // See https://en.wikipedia.org/wiki/Window_function#Gaussian_window
    12  // and https://www.recordingblogs.com/wiki/gaussian-window for details.
    13  //
    14  // The Gaussian window is an adjustable window.
    15  //
    16  // The sequence weights are
    17  //
    18  //	w[k] = exp(-0.5 * ((k-M)/(σ*M))² ), M = (N-1)/2,
    19  //
    20  // for k=0,1,...,N-1 where N is the length of the window.
    21  //
    22  // The properties of the window depend on the value of σ (sigma).
    23  // It can be used as high or low resolution window, depending of the σ value.
    24  //
    25  // Spectral leakage parameters are summarized in the table:
    26  //
    27  //	       |  σ=0.3  |  σ=0.5 |  σ=1.2 |
    28  //	-------|---------------------------|
    29  //	ΔF_0   |   8     |   3.4  |   2.2  |
    30  //	ΔF_0.5 |   1.82  |   1.2  |   0.94 |
    31  //	K      |   4     |   1.7  |   1.1  |
    32  //	ɣ_max  | -65     | -31.5  | -15.5  |
    33  //	β      |  -8.52  |  -4.48 |  -0.96 |
    34  type Gaussian struct {
    35  	Sigma float64
    36  }
    37  
    38  // Transform applies the Gaussian transformation to seq in place, using the
    39  // value of the receiver as the sigma parameter, and returning the result.
    40  func (g Gaussian) Transform(seq []float64) []float64 {
    41  	a := float64(len(seq)-1) / 2
    42  	for i := range seq {
    43  		x := -0.5 * math.Pow((float64(i)-a)/(g.Sigma*a), 2)
    44  		seq[i] *= math.Exp(x)
    45  	}
    46  	return seq
    47  }
    48  
    49  // TransformComplex applies the Gaussian transformation to seq in place,
    50  // using the value of the receiver as the sigma parameter, and returning
    51  // the result.
    52  func (g Gaussian) TransformComplex(seq []complex128) []complex128 {
    53  	a := float64(len(seq)-1) / 2
    54  	for i, v := range seq {
    55  		x := -0.5 * math.Pow((float64(i)-a)/(g.Sigma*a), 2)
    56  		w := math.Exp(x)
    57  		seq[i] = complex(w*real(v), w*imag(v))
    58  	}
    59  	return seq
    60  }
    61  
    62  // Tukey can modify a sequence using the Tukey window and return the result.
    63  // See https://en.wikipedia.org/wiki/Window_function#Tukey_window
    64  // and https://www.recordingblogs.com/wiki/tukey-window for details.
    65  //
    66  // The Tukey window is an adjustable window.
    67  //
    68  // The sequence weights are
    69  //
    70  //	w[k] = 0.5 * (1 + cos(π*(|k - M| - αM)/((1-α) * M))), |k - M| ≥ αM
    71  //	     = 1, |k - M| < αM
    72  //
    73  // with M = (N - 1)/2 for k=0,1,...,N-1 where N is the length of the window.
    74  //
    75  // Spectral leakage parameters are summarized in the table:
    76  //
    77  //	       |  α=0.3 |  α=0.5 |  α=0.7 |
    78  //	-------|--------------------------|
    79  //	ΔF_0   |   1.33 |   1.22 |   1.13 |
    80  //	ΔF_0.5 |   1.28 |   1.16 |   1.04 |
    81  //	K      |   0.67 |   0.61 |   0.57 |
    82  //	ɣ_max  | -18.2  | -15.1  | -13.8  |
    83  //	β      |  -1.41 |  -2.50 |  -3.74 |
    84  type Tukey struct {
    85  	Alpha float64
    86  }
    87  
    88  // Transform applies the Tukey transformation to seq in place, using the
    89  // value of the receiver as the Alpha parameter, and returning the result.
    90  func (t Tukey) Transform(seq []float64) []float64 {
    91  	switch {
    92  	case t.Alpha <= 0:
    93  		return Rectangular(seq)
    94  	case t.Alpha >= 1:
    95  		return Hann(seq)
    96  	default:
    97  		alphaL := t.Alpha * float64(len(seq)-1)
    98  		width := int(0.5*alphaL) + 1
    99  		for i := range seq[:width] {
   100  			w := 0.5 * (1 - math.Cos(2*math.Pi*float64(i)/alphaL))
   101  			seq[i] *= w
   102  			seq[len(seq)-1-i] *= w
   103  		}
   104  		return seq
   105  	}
   106  }
   107  
   108  // TransformComplex applies the Tukey transformation to seq in place, using
   109  // the value of the receiver as the Alpha parameter, and returning the result.
   110  func (t Tukey) TransformComplex(seq []complex128) []complex128 {
   111  	switch {
   112  	case t.Alpha <= 0:
   113  		return RectangularComplex(seq)
   114  	case t.Alpha >= 1:
   115  		return HannComplex(seq)
   116  	default:
   117  		alphaL := t.Alpha * float64(len(seq)-1)
   118  		width := int(0.5*alphaL) + 1
   119  		for i, v := range seq[:width] {
   120  			w := 0.5 * (1 - math.Cos(2*math.Pi*float64(i)/alphaL))
   121  			v = complex(w*real(v), w*imag(v))
   122  			seq[i] = v
   123  			seq[len(seq)-1-i] = v
   124  		}
   125  		return seq
   126  	}
   127  }
   128  
   129  // Values is an arbitrary real window function.
   130  type Values []float64
   131  
   132  // NewValues returns a Values of length n with weights corresponding to the
   133  // provided window function.
   134  func NewValues(window func([]float64) []float64, n int) Values {
   135  	v := make(Values, n)
   136  	for i := range v {
   137  		v[i] = 1
   138  	}
   139  	return window(v)
   140  }
   141  
   142  // Transform applies the weights in the receiver to seq in place, returning the
   143  // result. If v is nil, Transform is a no-op, otherwise the length of v must
   144  // match the length of seq.
   145  func (v Values) Transform(seq []float64) []float64 {
   146  	if v == nil {
   147  		return seq
   148  	}
   149  	if len(v) != len(seq) {
   150  		panic("window: length mismatch")
   151  	}
   152  	for i, w := range v {
   153  		seq[i] *= w
   154  	}
   155  	return seq
   156  }
   157  
   158  // TransformTo applies the weights in the receiver to src placing the result
   159  // in dst. If v is nil, TransformTo is a no-op, otherwise the length of v must
   160  // match the length of src and dst.
   161  func (v Values) TransformTo(dst, src []float64) {
   162  	if v == nil {
   163  		return
   164  	}
   165  	if len(v) != len(src) {
   166  		panic("window: seq length mismatch")
   167  	}
   168  	if len(v) != len(dst) {
   169  		panic("window: dst length mismatch")
   170  	}
   171  	for i, w := range v {
   172  		dst[i] = w * src[i]
   173  	}
   174  }
   175  
   176  // TransformComplex applies the weights in the receiver to seq in place,
   177  // returning the result. If v is nil, TransformComplex is a no-op, otherwise
   178  // the length of v must match the length of seq.
   179  func (v Values) TransformComplex(seq []complex128) []complex128 {
   180  	if v == nil {
   181  		return seq
   182  	}
   183  	if len(v) != len(seq) {
   184  		panic("window: length mismatch")
   185  	}
   186  	for i, w := range v {
   187  		sv := seq[i]
   188  		seq[i] = complex(w*real(sv), w*imag(sv))
   189  	}
   190  	return seq
   191  }
   192  
   193  // TransformComplexTo applies the weights in the receiver to src placing the
   194  // result in dst. If v is nil, TransformComplexTo is a no-op, otherwise the
   195  // length of v must match the length of src and dst.
   196  func (v Values) TransformComplexTo(dst, src []complex128) {
   197  	if v == nil {
   198  		return
   199  	}
   200  	if len(v) != len(src) {
   201  		panic("window: seq length mismatch")
   202  	}
   203  	if len(v) != len(dst) {
   204  		panic("window: dst length mismatch")
   205  	}
   206  	for i, w := range v {
   207  		sv := src[i]
   208  		dst[i] = complex(w*real(sv), w*imag(sv))
   209  	}
   210  }