From d7a8d4ec742a01ee04c27e2fd1dc0a558c8ace64 Mon Sep 17 00:00:00 2001
From: "J. David Lee" <johnl@crumpington.com>
Date: Thu, 5 Dec 2019 21:37:16 +0100
Subject: [PATCH] WIP

---
 lib/extractfreq/freq.go           | 19 +++++++++
 lib/extractfreq/harmonicwindow.go | 71 ++++++++++++++++++-------------
 lib/extractfreq/wip.go            | 52 +++++++++++++---------
 3 files changed, 93 insertions(+), 49 deletions(-)
 create mode 100644 lib/extractfreq/freq.go

diff --git a/lib/extractfreq/freq.go b/lib/extractfreq/freq.go
new file mode 100644
index 0000000..b4427bd
--- /dev/null
+++ b/lib/extractfreq/freq.go
@@ -0,0 +1,19 @@
+package extractfreq
+
+type ExtractFreqArgs struct {
+	Data       []float64
+	ApproxFreq float64
+	Log        bool
+	LogPrefix  string
+	LogDir     string
+}
+
+/*
+func ExtractFreq(args ExtractFreqArgs) float64 {
+	fft := fourier.NewFFT(len(args.Data))
+	cCoeffs := fft.Coefficients(nil, args.Data)
+	amps := make([]float64, len(cCoeffs))
+	freqs := make([]float64, len(cCoeffs))
+	ampMax := float64(0)
+}
+*/
diff --git a/lib/extractfreq/harmonicwindow.go b/lib/extractfreq/harmonicwindow.go
index 034350e..2dc9100 100644
--- a/lib/extractfreq/harmonicwindow.go
+++ b/lib/extractfreq/harmonicwindow.go
@@ -1,17 +1,18 @@
 package extractfreq
 
 import (
+	"log"
 	"math"
 
 	"git.crumpington.com/public/jlaudio/lib/util"
+	"gonum.org/v1/gonum/floats"
 	"gonum.org/v1/gonum/stat"
 )
 
 type HarmonicWindow struct {
-	TargetFreq     float64
+	Harmonic       float64 // 0.25, 0.5, 1, 2, 3, ...
 	Amps           []float64
 	Freqs          []float64
-	Idx0           int
 	IdxPeak        int
 	IdxPeakInterp  float64
 	FreqPeakInterp float64
@@ -20,37 +21,47 @@ type HarmonicWindow struct {
 	AmpStd         float64
 }
 
-// Construct a window around the harmonic peak with a width of 100 cents.
 func NewHarmonicWindow(
-	amps []float64,
-	freqs []float64,
-	targetFreq float64,
-	idxPeak int,
+	fFundamental float64,
+	harmonic float64,
+	rawAmps []float64,
+	rawFreqs []float64,
 ) *HarmonicWindow {
-	deltaF := freqs[1] - freqs[0] // Size of frequency bins.
-	fCenter := freqs[idxPeak]     // Center frequency.
+	deltaF := rawFreqs[1] - rawFreqs[0]
 
-	// Create a window approximately 200 cents in width.
-	deltaIdx := int(math.Round(deltaF / (100 * util.CentDeltaFreq(fCenter))))
-	idxMin := idxPeak - deltaIdx
-	idxMax := idxPeak + deltaIdx
+	// TODO: At higher harmonic numbers, the harmonics are going to be spaced
+	// more cloesly together. That means that the window size must shrink.
+	winSizeCents := float64(150)
 
-	if idxMin < 0 {
-		idxMin = 0
+	// To start, we create a window around the approximate harmonic location.
+	idxPeak := int((fFundamental * harmonic) / deltaF)
+	winSize := int(
+		math.Round(
+			winSizeCents * util.CentDeltaFreq(fFundamental*harmonic) /
+				deltaF))
+	amps := rawAmps[idxPeak-winSize : idxPeak+winSize]
+	freqs := rawFreqs[idxPeak-winSize : idxPeak+winSize]
+
+	// Next, we refine the window be centering it on the actual peak.
+	idxPeak += (floats.MaxIdx(amps) - winSize)
+	amps = rawAmps[idxPeak-winSize : idxPeak+winSize]
+	freqs = rawFreqs[idxPeak-winSize : idxPeak+winSize]
+	idxPeak = winSize
+
+	// Next we compute statistical properties to determine if this peak is worth
+	// using.
+	mean, std := stat.MeanStdDev(amps, nil)
+	ampPeak := amps[idxPeak]
+
+	log.Printf("%0.2f: %0.4f", harmonic, (ampPeak-mean)/std)
+
+	return &HarmonicWindow{
+		Harmonic: harmonic,
+		Amps:     amps,
+		Freqs:    freqs,
+		IdxPeak:  idxPeak,
+		AmpPeak:  ampPeak,
+		AmpMean:  mean,
+		AmpStd:   std,
 	}
-
-	if idxMax > len(freqs) {
-		idxMax = len(freqs)
-	}
-
-	hw := &HarmonicWindow{
-		TargetFreq: targetFreq,
-		Amps:       amps[idxMin:idxMax],
-		Freqs:      freqs[idxMin:idxMax],
-		IdxPeak:    idxPeak,
-		AmpPeak:    amps[idxPeak],
-	}
-
-	hw.AmpMean, hw.AmpStd = stat.MeanStdDev(hw.Amps, nil)
-	return hw
 }
diff --git a/lib/extractfreq/wip.go b/lib/extractfreq/wip.go
index 62f8268..a1c3097 100644
--- a/lib/extractfreq/wip.go
+++ b/lib/extractfreq/wip.go
@@ -1,6 +1,7 @@
 package extractfreq
 
 import (
+	"fmt"
 	"log"
 	"math/cmplx"
 
@@ -62,28 +63,41 @@ func PlotFFT(path string, nSamples int, targetFreq float64) {
 	*/
 	// Find frequency window around midi note.
 
-	p, err := plot.New()
-	if err != nil {
-		panic(err)
-	}
+	harmonics := []float64{0.25, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}
+	for _, harmonic := range harmonics {
 
-	log.Printf("Making points...")
-	pts := make(plotter.XYs, len(amps))
-	for i := range amps {
-		pts[i].X = freqs[i]
-		pts[i].Y = amps[i]
-	}
+		hw := NewHarmonicWindow(
+			targetFreq,
+			harmonic,
+			amps,
+			freqs)
 
-	log.Printf("Plotting...")
-	if err = plotutil.AddLines(p, "fft", pts); err != nil {
-		panic(err)
-	}
+		p, err := plot.New()
+		if err != nil {
+			panic(err)
+		}
 
-	p.Y.Scale = plot.LogScale{}
-	//p.X.Scale = plot.LogScale{}
+		log.Printf("Making points...")
+		pts := make(plotter.XYs, 0, len(hw.Amps))
+		for i := range hw.Amps {
+			pts = append(pts, plotter.XY{
+				X: hw.Freqs[i],
+				Y: hw.Amps[i],
+			})
+		}
 
-	log.Printf("Saving...")
-	if err := p.Save(8*vg.Inch, 8*vg.Inch, "plot.png"); err != nil {
-		panic(err)
+		log.Printf("Plotting...")
+		if err = plotutil.AddLines(p, "fft", pts); err != nil {
+			panic(err)
+		}
+
+		//p.Y.Scale = plot.LogScale{}
+		//p.X.Scale = plot.LogScale{}
+
+		log.Printf("Saving...")
+		filename := fmt.Sprintf("plot-%03.2f.png", harmonic)
+		if err := p.Save(8*vg.Inch, 8*vg.Inch, filename); err != nil {
+			panic(err)
+		}
 	}
 }