WIP: extractfreq: a-weighting

lib/extractfreq/aweight.go

@@ -0,0 +1,9 @@
+package extractfreq
+
+import "math"
+
+func AWeight(f float64) float64 {
+	return (148693636.0 * f * f * f * f) / ((424.36 + f*f) * math.Sqrt(
+		(f*f+11599.29)*(f*f+544496.41)) *
+		(f*f + 148693636.0))
+}

lib/extractfreq/harmonicwindow.go

@@ -2,9 +2,7 @@ package extractfreq
 
 import (
 	"log"
-	"math"
 
-	"git.crumpington.com/public/jlaudio/lib/util"
 	"gonum.org/v1/gonum/floats"
 	"gonum.org/v1/gonum/stat"
 )
@@ -31,14 +29,11 @@ func NewHarmonicWindow(
 
 	// 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)
+	//winSizeCents := float64(150)
 
 	// To start, we create a window around the approximate harmonic location.
 	idxPeak := int((fFundamental * harmonic) / deltaF)
-	winSize := int(
+	winSize := int(0.5 * fFundamental / deltaF)
-		math.Round(
-			winSizeCents * util.CentDeltaFreq(fFundamental*harmonic) /
-				deltaF))
 	amps := rawAmps[idxPeak-winSize : idxPeak+winSize]
 	freqs := rawFreqs[idxPeak-winSize : idxPeak+winSize]
 
@@ -54,14 +49,17 @@ func NewHarmonicWindow(
 	ampPeak := amps[idxPeak]
 
 	log.Printf("%0.2f: %0.4f", harmonic, (ampPeak-mean)/std)
+	// TODO: Actual interpolation.
 
 	return &HarmonicWindow{
-		Harmonic: harmonic,
+		Harmonic:       harmonic,
-		Amps:     amps,
+		Amps:           amps,
-		Freqs:    freqs,
+		Freqs:          freqs,
-		IdxPeak:  idxPeak,
+		IdxPeak:        idxPeak,
-		AmpPeak:  ampPeak,
+		IdxPeakInterp:  float64(idxPeak),
-		AmpMean:  mean,
+		FreqPeakInterp: freqs[idxPeak],
-		AmpStd:   std,
+		AmpPeak:        ampPeak,
+		AmpMean:        mean,
+		AmpStd:         std,
 	}
 }

lib/extractfreq/wip.go

@@ -3,6 +3,7 @@ package extractfreq
 import (
 	"fmt"
 	"log"
+	"math"
 	"math/cmplx"
 
 	"git.crumpington.com/public/jlaudio/lib/flac"
@@ -13,7 +14,12 @@ import (
 	"gonum.org/v1/plot/vg"
 )
 
-func PlotFFT(path string, nSamples int, targetFreq float64) {
+// Steps:
+//
+// * Create amplitude spectrum (abs of FFT).
+// * Multiply by A-weight curve
+// * Take log10 of result
+func PlotFFT(path string, nSamples int, fFundamental float64) {
 	L, R, err := flac.Load(path)
 	if err != nil {
 		panic(err)
@@ -32,19 +38,21 @@ func PlotFFT(path string, nSamples int, targetFreq float64) {
 	cCoeffs := fft.Coefficients(nil, data)
 	amps := make([]float64, len(cCoeffs))
 	freqs := make([]float64, len(cCoeffs))
-	ampMax := float64(0)
 	for i := range freqs {
-		amps[i] = cmplx.Abs(cCoeffs[i])
-		if amps[i] > ampMax {
-			ampMax = amps[i]
-		}
 		freqs[i] = fft.Freq(i) * 48000
+		if freqs[i] < 20 {
+			amps[i] = 0
+		} else {
+			amps[i] = math.Log10(cmplx.Abs(cCoeffs[i]) * AWeight(freqs[i]))
+		}
 	}
 
 	// Normalize coefficients.
-	for i := range amps {
+	//for i := range amps {
-		amps[i] /= ampMax
+	//amps[i] /= ampMax
-	}
+	//}
+
+	// Find highest N peaks and sort by size.
 
 	/*
 		// Use a moving window to search for local peaks.
@@ -63,11 +71,42 @@ func PlotFFT(path string, nSamples int, targetFreq float64) {
 	*/
 	// Find frequency window around midi note.
 
-	harmonics := []float64{0.25, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}
+	p, err := plot.New()
+	if err != nil {
+		panic(err)
+	}
+
+	log.Printf("Making points...")
+	pts := make(plotter.XYs, 0, len(amps))
+	for i := range amps {
+		if freqs[i] < 20 || freqs[i] > 16000 {
+			continue
+		}
+		pts = append(pts, plotter.XY{
+			X: freqs[i],
+			Y: amps[i],
+		})
+	}
+
+	log.Printf("Plotting...")
+	if err = plotutil.AddLinePoints(p, "fft", pts); err != nil {
+		panic(err)
+	}
+
+	//p.Y.Scale = plot.LogScale{}
+	p.X.Scale = plot.LogScale{}
+
+	log.Printf("Saving...")
+	if err := p.Save(16*vg.Inch, 8*vg.Inch, "fft.png"); err != nil {
+		panic(err)
+	}
+
+	harmonics := []float64{0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}
+	hws := []*HarmonicWindow{}
 	for _, harmonic := range harmonics {
 
 		hw := NewHarmonicWindow(
-			targetFreq,
+			fFundamental,
 			harmonic,
 			amps,
 			freqs)
@@ -99,5 +138,27 @@ func PlotFFT(path string, nSamples int, targetFreq float64) {
 		if err := p.Save(8*vg.Inch, 8*vg.Inch, filename); err != nil {
 			panic(err)
 		}
+
+		hws = append(hws, hw)
+	}
+
+	pts = make(plotter.XYs, len(hws))
+	for i, hw := range hws[:11] {
+		pts[i].X = hw.Harmonic
+		pts[i].Y = hw.FreqPeakInterp
+	}
+
+	p, err = plot.New()
+	if err != nil {
+		panic(err)
+	}
+
+	if err = plotutil.AddScatters(p, "", pts); err != nil {
+		panic(err)
+	}
+
+	log.Printf("Saving...")
+	if err := p.Save(8*vg.Inch, 8*vg.Inch, "harmonics.png"); err != nil {
+		panic(err)
 	}
 }