gonum.org/v1/gonum@v0.14.0/dsp/window/window_example_test.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_test
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  	"math/cmplx"
    11  
    12  	"gonum.org/v1/gonum/dsp/fourier"
    13  	"gonum.org/v1/gonum/dsp/window"
    14  )
    15  
    16  func Example() {
    17  	// The input sequence is a 2.5 period of the Sin function.
    18  	src := make([]float64, 20)
    19  	k := 5 * math.Pi / float64(len(src)-1)
    20  	for i := range src {
    21  		src[i] = math.Sin(k * float64(i))
    22  	}
    23  
    24  	// Initialize an FFT and perform the analysis.
    25  	fft := fourier.NewFFT(len(src))
    26  	coeff := fft.Coefficients(nil, src)
    27  
    28  	// The result shows that width of the main lobe with center
    29  	// between frequencies 0.1 and 0.15 is small, but that the
    30  	// height of the side lobes is large.
    31  	fmt.Println("Rectangular window (or no window):")
    32  	for i, c := range coeff {
    33  		fmt.Printf("freq=%.4f\tcycles/period, magnitude=%.4f,\tphase=%.4f\n",
    34  			fft.Freq(i), cmplx.Abs(c), cmplx.Phase(c))
    35  	}
    36  
    37  	// Initialize an FFT and perform the analysis on a sequence
    38  	// transformed by the Hamming window function.
    39  	fft = fourier.NewFFT(len(src))
    40  	coeff = fft.Coefficients(nil, window.Hamming(src))
    41  
    42  	// The result shows that width of the main lobe is wider,
    43  	// but height of the side lobes is lower.
    44  	fmt.Println("Hamming window:")
    45  	// The magnitude of all bins has been decreased by β.
    46  	// Generally in an analysis amplification may be omitted, but to
    47  	// make a comparable data, the result should be amplified by -β
    48  	// of the window function — +5.37 dB for the Hamming window.
    49  	//  -β = 20 log_10(amplifier).
    50  	amplifier := math.Pow(10, 5.37/20.0)
    51  	for i, c := range coeff {
    52  		fmt.Printf("freq=%.4f\tcycles/period, magnitude=%.4f,\tphase=%.4f\n",
    53  			fft.Freq(i), amplifier*cmplx.Abs(c), cmplx.Phase(c))
    54  	}
    55  	// Output:
    56  	//
    57  	// Rectangular window (or no window):
    58  	// freq=0.0000	cycles/period, magnitude=2.2798,	phase=0.0000
    59  	// freq=0.0500	cycles/period, magnitude=2.6542,	phase=0.1571
    60  	// freq=0.1000	cycles/period, magnitude=5.3115,	phase=0.3142
    61  	// freq=0.1500	cycles/period, magnitude=7.3247,	phase=-2.6704
    62  	// freq=0.2000	cycles/period, magnitude=1.6163,	phase=-2.5133
    63  	// freq=0.2500	cycles/period, magnitude=0.7681,	phase=-2.3562
    64  	// freq=0.3000	cycles/period, magnitude=0.4385,	phase=-2.1991
    65  	// freq=0.3500	cycles/period, magnitude=0.2640,	phase=-2.0420
    66  	// freq=0.4000	cycles/period, magnitude=0.1530,	phase=-1.8850
    67  	// freq=0.4500	cycles/period, magnitude=0.0707,	phase=-1.7279
    68  	// freq=0.5000	cycles/period, magnitude=0.0000,	phase=0.0000
    69  	// Hamming window:
    70  	// freq=0.0000	cycles/period, magnitude=0.0542,	phase=3.1416
    71  	// freq=0.0500	cycles/period, magnitude=0.8458,	phase=-2.9845
    72  	// freq=0.1000	cycles/period, magnitude=7.1519,	phase=0.3142
    73  	// freq=0.1500	cycles/period, magnitude=8.5907,	phase=-2.6704
    74  	// freq=0.2000	cycles/period, magnitude=2.0804,	phase=0.6283
    75  	// freq=0.2500	cycles/period, magnitude=0.0816,	phase=0.7854
    76  	// freq=0.3000	cycles/period, magnitude=0.0156,	phase=-2.1991
    77  	// freq=0.3500	cycles/period, magnitude=0.0224,	phase=-2.0420
    78  	// freq=0.4000	cycles/period, magnitude=0.0163,	phase=-1.8850
    79  	// freq=0.4500	cycles/period, magnitude=0.0083,	phase=-1.7279
    80  	// freq=0.5000	cycles/period, magnitude=0.0000,	phase=0.0000
    81  }
    82  
    83  func ExampleHamming() {
    84  	src := []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    85  		1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
    86  
    87  	// Window functions change data in place. So, if input data
    88  	// needs to stay unchanged, it must be copied.
    89  	srcCpy := append([]float64(nil), src...)
    90  	// Apply window function to srcCpy.
    91  	dst := window.Hamming(srcCpy)
    92  
    93  	// src is unchanged.
    94  	fmt.Printf("src:    %f\n", src)
    95  	// srcCpy is altered.
    96  	fmt.Printf("srcCpy: %f\n", srcCpy)
    97  	// dst mirrors the srcCpy slice.
    98  	fmt.Printf("dst:    %f\n", dst)
    99  
   100  	// Output:
   101  	//
   102  	// 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]
   103  	// srcCpy: [0.080000 0.104924 0.176995 0.288404 0.427077 0.577986 0.724780 0.851550 0.944558 0.993726 0.993726 0.944558 0.851550 0.724780 0.577986 0.427077 0.288404 0.176995 0.104924 0.080000]
   104  	// dst:    [0.080000 0.104924 0.176995 0.288404 0.427077 0.577986 0.724780 0.851550 0.944558 0.993726 0.993726 0.944558 0.851550 0.724780 0.577986 0.427077 0.288404 0.176995 0.104924 0.080000]
   105  }
   106  
   107  func ExampleValues() {
   108  	src := []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
   109  		1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
   110  
   111  	// Create a Sine Window lookup table.
   112  	sine := window.NewValues(window.Sine, len(src))
   113  
   114  	// Apply the transformation to the src.
   115  	fmt.Printf("dst: %f\n", sine.Transform(src))
   116  
   117  	// Output:
   118  	//
   119  	// dst: [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]
   120  }
   121  
   122  func ExampleValues_TransformTo_gabor() {
   123  	src := []float64{1, 2, 1, 0, -1, -1, -2, -2, -1, -1,
   124  		0, 1, 1, 2, 1, 0, -1, -2, -1, 0}
   125  
   126  	// Create a Gaussian Window lookup table for 4 samples.
   127  	gaussian := window.NewValues(window.Gaussian{0.5}.Transform, 4)
   128  
   129  	// Prepare a destination.
   130  	dst := make([]float64, 8)
   131  
   132  	// Apply the transformation to the src, placing it in dst.
   133  	for i := 0; i < len(src)-len(gaussian); i++ {
   134  		gaussian.TransformTo(dst[0:len(gaussian)], src[i:i+len(gaussian)])
   135  
   136  		// To perform the Gabor transform, we would calculate
   137  		// the FFT on dst for each iteration.
   138  		fmt.Printf("FFT(%f)\n", dst)
   139  	}
   140  
   141  	// Output:
   142  	//
   143  	// FFT([0.135335 1.601475 0.800737 0.000000 0.000000 0.000000 0.000000 0.000000])
   144  	// FFT([0.270671 0.800737 0.000000 -0.135335 0.000000 0.000000 0.000000 0.000000])
   145  	// FFT([0.135335 0.000000 -0.800737 -0.135335 0.000000 0.000000 0.000000 0.000000])
   146  	// FFT([0.000000 -0.800737 -0.800737 -0.270671 0.000000 0.000000 0.000000 0.000000])
   147  	// FFT([-0.135335 -0.800737 -1.601475 -0.270671 0.000000 0.000000 0.000000 0.000000])
   148  	// FFT([-0.135335 -1.601475 -1.601475 -0.135335 0.000000 0.000000 0.000000 0.000000])
   149  	// FFT([-0.270671 -1.601475 -0.800737 -0.135335 0.000000 0.000000 0.000000 0.000000])
   150  	// FFT([-0.270671 -0.800737 -0.800737 0.000000 0.000000 0.000000 0.000000 0.000000])
   151  	// FFT([-0.135335 -0.800737 0.000000 0.135335 0.000000 0.000000 0.000000 0.000000])
   152  	// FFT([-0.135335 0.000000 0.800737 0.135335 0.000000 0.000000 0.000000 0.000000])
   153  	// FFT([0.000000 0.800737 0.800737 0.270671 0.000000 0.000000 0.000000 0.000000])
   154  	// FFT([0.135335 0.800737 1.601475 0.135335 0.000000 0.000000 0.000000 0.000000])
   155  	// FFT([0.135335 1.601475 0.800737 0.000000 0.000000 0.000000 0.000000 0.000000])
   156  	// FFT([0.270671 0.800737 0.000000 -0.135335 0.000000 0.000000 0.000000 0.000000])
   157  	// FFT([0.135335 0.000000 -0.800737 -0.270671 0.000000 0.000000 0.000000 0.000000])
   158  	// FFT([0.000000 -0.800737 -1.601475 -0.135335 0.000000 0.000000 0.000000 0.000000])
   159  }