WIP: peak finding

lib/extractfreq/aweight.go

@@ -2,8 +2,11 @@ package extractfreq
 
 import "math"
 
+// A-weight filter is supposed to correspond to how a human would perceive
+// loudness.
 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))
+	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/wip.go

@@ -16,9 +16,12 @@ import (
 
 // Steps:
 //
-// * Create amplitude spectrum (abs of FFT).
+// * Create amplitude spectrum (abs of real FFT).
 // * Multiply by A-weight curve
-// * Take log10 of result
+// * Find peaks
+// * Sort peaks by amplitude
+// * Use strongest peak to estimate fFundamental?
+// *
 func PlotFFT(path string, nSamples int, fFundamental float64) {
 	L, R, err := flac.Load(path)
 	if err != nil {
@@ -34,98 +37,96 @@ func PlotFFT(path string, nSamples int, fFundamental float64) {
 		data[i] = float64(L[i]) + float64(R[i])
 	}
 
+	// Make initial fft.
 	fft := fourier.NewFFT(nSamples)
 	cCoeffs := fft.Coefficients(nil, data)
 	amps := make([]float64, len(cCoeffs))
 	freqs := make([]float64, len(cCoeffs))
 	for i := range freqs {
 		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]))
-		}
+		amps[i] = cmplx.Abs(cCoeffs[i])
 	}
 
-	// Normalize coefficients.
-	//for i := range amps {
-	//amps[i] /= ampMax
-	//}
+	// Trim to frequencies of interest.
+	deltaF := freqs[1] - freqs[0]
+	idx0 := int(20 / deltaF)
+	idx1 := int(16000 / deltaF)
+	freqs = freqs[idx0:idx1]
+	amps = amps[idx0:idx1]
 
-	// Find highest N peaks and sort by size.
+	// Multiply by A-weight function.
+	for i := range amps {
+		amps[i] *= AWeight(freqs[i])
+	}
 
-	/*
-		// Use a moving window to search for local peaks.
-		hws := []*HarmonicWindow{}
-
-		// Bypass low-frequency samples.
-		i := 0
-		for i < len(freqs) && freqs[i] < 24 {
-			i++
-		}
-
-		iCenter := i
-		for iC{
-
-		}
-	*/
-	// Find frequency window around midi note.
+	// Find peaks.
+	dFreq := freqs[1] - freqs[0]
+	halfWinSize := int(math.Round((fFundamental / dFreq) / 3))
+	peakInds, peakAmps := getPeaks(halfWinSize, amps)
+	peakInds, peakAmps = sortPeaks(peakInds, peakAmps)
 
 	p, err := plot.New()
 	if err != nil {
 		panic(err)
 	}
 
+	fMin := float64(20)
+	fMax := float64(16000)
+
 	log.Printf("Making points...")
 	pts := make(plotter.XYs, 0, len(amps))
 	for i := range amps {
-		if freqs[i] < 20 || freqs[i] > 16000 {
-			continue
+		if freqs[i] > fMin && freqs[i] < fMax {
+			pts = append(pts, plotter.XY{
+				X: freqs[i],
+				Y: amps[i],
+			})
+		}
+	}
+
+	pts2 := make(plotter.XYs, 0, len(amps))
+	for i := range peakInds {
+		if freqs[peakInds[i]] > fMin && freqs[peakInds[i]] < fMax {
+			pts2 = append(pts2, plotter.XY{
+				X: freqs[peakInds[i]],
+				Y: peakAmps[i],
+			})
 		}
-		pts = append(pts, plotter.XY{
-			X: freqs[i],
-			Y: amps[i],
-		})
 	}
 
 	log.Printf("Plotting...")
-	if err = plotutil.AddLinePoints(p, "fft", pts); err != nil {
+	if err = plotutil.AddScatters(p, "fft", pts, "max", pts2); err != nil {
 		panic(err)
 	}
-
-	//p.Y.Scale = plot.LogScale{}
-	p.X.Scale = plot.LogScale{}
+	p.Add(plotter.NewGrid())
 
 	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 {
+	for i, idx := range peakInds {
+		//amp := peakAmps[i]
 
-		hw := NewHarmonicWindow(
-			fFundamental,
-			harmonic,
-			amps,
-			freqs)
+		// Make window.
+		ampsWin := amps[idx-halfWinSize : idx+halfWinSize]
+		freqWin := freqs[idx-halfWinSize : idx+halfWinSize]
 
+		log.Printf("Making points...")
+		pts := make(plotter.XYs, 0, 2*halfWinSize)
+		for i := range ampsWin {
+			pts = append(pts, plotter.XY{
+				X: freqWin[i],
+				Y: ampsWin[i],
+			})
+		}
+
+		log.Printf("Plotting...")
 		p, err := plot.New()
 		if err != nil {
 			panic(err)
 		}
 
-		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("Plotting...")
 		if err = plotutil.AddLines(p, "fft", pts); err != nil {
 			panic(err)
 		}
@@ -134,31 +135,9 @@ func PlotFFT(path string, nSamples int, fFundamental float64) {
 		//p.X.Scale = plot.LogScale{}
 
 		log.Printf("Saving...")
-		filename := fmt.Sprintf("plot-%03.2f.png", harmonic)
+		filename := fmt.Sprintf("plot-%03d.png", i)
 		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)
 	}
 }