diff --git a/lib/extractfreq/aweight.go b/lib/extractfreq/aweight.go index dc424df..86a156c 100644 --- a/lib/extractfreq/aweight.go +++ b/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)) } diff --git a/lib/extractfreq/wip.go b/lib/extractfreq/wip.go index 76771a7..0fa9ff4 100644 --- a/lib/extractfreq/wip.go +++ b/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) } }