From 78c90707635d1cc0ae4bc70a5cfd9292445fd82d Mon Sep 17 00:00:00 2001 From: "J. David Lee" Date: Tue, 3 Dec 2019 21:34:29 +0100 Subject: [PATCH] WIP: extractfreq --- lib/extractfreq/harmonicwindow.go | 56 +++++++++++++++++++ lib/extractfreq/wip.go | 89 +++++++++++++++++++++++++++++++ lib/extractfreq/wip_test.go | 7 +++ lib/util/notefreq.go | 14 ++++- lib/util/notefreq_test.go | 47 +++++++++++++++- 5 files changed, 210 insertions(+), 3 deletions(-) create mode 100644 lib/extractfreq/harmonicwindow.go create mode 100644 lib/extractfreq/wip.go create mode 100644 lib/extractfreq/wip_test.go diff --git a/lib/extractfreq/harmonicwindow.go b/lib/extractfreq/harmonicwindow.go new file mode 100644 index 0000000..034350e --- /dev/null +++ b/lib/extractfreq/harmonicwindow.go @@ -0,0 +1,56 @@ +package extractfreq + +import ( + "math" + + "git.crumpington.com/public/jlaudio/lib/util" + "gonum.org/v1/gonum/stat" +) + +type HarmonicWindow struct { + TargetFreq float64 + Amps []float64 + Freqs []float64 + Idx0 int + IdxPeak int + IdxPeakInterp float64 + FreqPeakInterp float64 + AmpPeak float64 + AmpMean float64 + 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, +) *HarmonicWindow { + deltaF := freqs[1] - freqs[0] // Size of frequency bins. + fCenter := freqs[idxPeak] // Center frequency. + + // Create a window approximately 200 cents in width. + deltaIdx := int(math.Round(deltaF / (100 * util.CentDeltaFreq(fCenter)))) + idxMin := idxPeak - deltaIdx + idxMax := idxPeak + deltaIdx + + if idxMin < 0 { + idxMin = 0 + } + + 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 new file mode 100644 index 0000000..62f8268 --- /dev/null +++ b/lib/extractfreq/wip.go @@ -0,0 +1,89 @@ +package extractfreq + +import ( + "log" + "math/cmplx" + + "git.crumpington.com/public/jlaudio/lib/flac" + "gonum.org/v1/gonum/fourier" + "gonum.org/v1/plot" + "gonum.org/v1/plot/plotter" + "gonum.org/v1/plot/plotutil" + "gonum.org/v1/plot/vg" +) + +func PlotFFT(path string, nSamples int, targetFreq float64) { + L, R, err := flac.Load(path) + if err != nil { + panic(err) + } + + if nSamples < len(L) { + nSamples = len(L) + } + + data := make([]float64, nSamples) + for i := range data { + data[i] = float64(L[i]) + float64(R[i]) + } + + fft := fourier.NewFFT(nSamples) + 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 + } + + // Normalize coefficients. + for i := range amps { + amps[i] /= ampMax + } + + /* + // 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. + + p, err := plot.New() + if err != nil { + panic(err) + } + + log.Printf("Making points...") + pts := make(plotter.XYs, len(amps)) + for i := range amps { + pts[i].X = freqs[i] + pts[i].Y = amps[i] + } + + 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...") + if err := p.Save(8*vg.Inch, 8*vg.Inch, "plot.png"); err != nil { + panic(err) + } +} diff --git a/lib/extractfreq/wip_test.go b/lib/extractfreq/wip_test.go new file mode 100644 index 0000000..4cf6370 --- /dev/null +++ b/lib/extractfreq/wip_test.go @@ -0,0 +1,7 @@ +package extractfreq + +import "testing" + +func TestPlotFFT(t *testing.T) { + PlotFFT("test_files/440.flac", 96000, 440) +} diff --git a/lib/util/notefreq.go b/lib/util/notefreq.go index 9f85cd0..316a93b 100644 --- a/lib/util/notefreq.go +++ b/lib/util/notefreq.go @@ -3,6 +3,16 @@ package util import "math" // Return the frequency of a given MIDI note number. -func NoteFreq(n int) float64 { - return (math.Pow(2, (float64(n)-69)/12)) * 440 +func NoteFreq(n float64) float64 { + return (math.Pow(2, (n-69)/12)) * 440 +} + +// Return the midi note number for the given frequency. +func FreqNote(f float64) float64 { + return 12*math.Log2(f/440) + 69 +} + +// Return the one-cent change in frequency at the given frequency. +func CentDeltaFreq(f float64) float64 { + return f * (math.Pow(2, 1/1200.0) - 1) } diff --git a/lib/util/notefreq_test.go b/lib/util/notefreq_test.go index 8e3bac6..c59848b 100644 --- a/lib/util/notefreq_test.go +++ b/lib/util/notefreq_test.go @@ -7,7 +7,7 @@ import ( func TestNoteFreq(t *testing.T) { type TestCase struct { - In int + In float64 Out float64 } @@ -29,3 +29,48 @@ func TestNoteFreq(t *testing.T) { } } } + +func TestFreqNote(t *testing.T) { + type TestCase struct { + Out float64 + In float64 + } + + testCases := []TestCase{ + {108, 4186.01}, + {81, 880}, + {69, 440}, + {58, 233.08}, + {57, 220}, + {56, 207.65}, + {55, 196}, + {36, 65.41}, + } + + for _, tc := range testCases { + out := FreqNote(tc.In) + if fmt.Sprintf("%.2f", out) != fmt.Sprintf("%.2f", tc.Out) { + t.Fatalf("%v != %v", out, tc.Out) + } + } +} + +func TestCentDeltaFreq(t *testing.T) { + type TestCase struct { + In, Out float64 + } + + testCases := []TestCase{ + {27.5, 0.016}, + {55, 0.032}, + {130.8128, 0.076}, + {440, 0.254}, + } + + for _, tc := range testCases { + out := CentDeltaFreq(tc.In) + if fmt.Sprintf("%.3f", out) != fmt.Sprintf("%.3f", tc.Out) { + t.Fatalf("%v: %v != %v", tc.In, out, tc.Out) + } + } +}