diff --git a/lib/extractfreq/aweight.go b/lib/extractfreq/aweight.go new file mode 100644 index 0000000..dc424df --- /dev/null +++ b/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)) +} diff --git a/lib/extractfreq/harmonicwindow.go b/lib/extractfreq/harmonicwindow.go index 2dc9100..c211d77 100644 --- a/lib/extractfreq/harmonicwindow.go +++ b/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( - math.Round( - winSizeCents * util.CentDeltaFreq(fFundamental*harmonic) / - deltaF)) + winSize := int(0.5 * fFundamental / 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, - Amps: amps, - Freqs: freqs, - IdxPeak: idxPeak, - AmpPeak: ampPeak, - AmpMean: mean, - AmpStd: std, + Harmonic: harmonic, + Amps: amps, + Freqs: freqs, + IdxPeak: idxPeak, + IdxPeakInterp: float64(idxPeak), + FreqPeakInterp: freqs[idxPeak], + AmpPeak: ampPeak, + AmpMean: mean, + AmpStd: std, } } diff --git a/lib/extractfreq/wip.go b/lib/extractfreq/wip.go index a1c3097..76771a7 100644 --- a/lib/extractfreq/wip.go +++ b/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 { - amps[i] /= ampMax - } + //for i := range amps { + //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) } }