WIP
parent
78c9070763
commit
d7a8d4ec74
|
@ -0,0 +1,19 @@
|
||||||
|
package extractfreq
|
||||||
|
|
||||||
|
type ExtractFreqArgs struct {
|
||||||
|
Data []float64
|
||||||
|
ApproxFreq float64
|
||||||
|
Log bool
|
||||||
|
LogPrefix string
|
||||||
|
LogDir string
|
||||||
|
}
|
||||||
|
|
||||||
|
/*
|
||||||
|
func ExtractFreq(args ExtractFreqArgs) float64 {
|
||||||
|
fft := fourier.NewFFT(len(args.Data))
|
||||||
|
cCoeffs := fft.Coefficients(nil, args.Data)
|
||||||
|
amps := make([]float64, len(cCoeffs))
|
||||||
|
freqs := make([]float64, len(cCoeffs))
|
||||||
|
ampMax := float64(0)
|
||||||
|
}
|
||||||
|
*/
|
|
@ -1,17 +1,18 @@
|
||||||
package extractfreq
|
package extractfreq
|
||||||
|
|
||||||
import (
|
import (
|
||||||
|
"log"
|
||||||
"math"
|
"math"
|
||||||
|
|
||||||
"git.crumpington.com/public/jlaudio/lib/util"
|
"git.crumpington.com/public/jlaudio/lib/util"
|
||||||
|
"gonum.org/v1/gonum/floats"
|
||||||
"gonum.org/v1/gonum/stat"
|
"gonum.org/v1/gonum/stat"
|
||||||
)
|
)
|
||||||
|
|
||||||
type HarmonicWindow struct {
|
type HarmonicWindow struct {
|
||||||
TargetFreq float64
|
Harmonic float64 // 0.25, 0.5, 1, 2, 3, ...
|
||||||
Amps []float64
|
Amps []float64
|
||||||
Freqs []float64
|
Freqs []float64
|
||||||
Idx0 int
|
|
||||||
IdxPeak int
|
IdxPeak int
|
||||||
IdxPeakInterp float64
|
IdxPeakInterp float64
|
||||||
FreqPeakInterp float64
|
FreqPeakInterp float64
|
||||||
|
@ -20,37 +21,47 @@ type HarmonicWindow struct {
|
||||||
AmpStd float64
|
AmpStd float64
|
||||||
}
|
}
|
||||||
|
|
||||||
// Construct a window around the harmonic peak with a width of 100 cents.
|
|
||||||
func NewHarmonicWindow(
|
func NewHarmonicWindow(
|
||||||
amps []float64,
|
fFundamental float64,
|
||||||
freqs []float64,
|
harmonic float64,
|
||||||
targetFreq float64,
|
rawAmps []float64,
|
||||||
idxPeak int,
|
rawFreqs []float64,
|
||||||
) *HarmonicWindow {
|
) *HarmonicWindow {
|
||||||
deltaF := freqs[1] - freqs[0] // Size of frequency bins.
|
deltaF := rawFreqs[1] - rawFreqs[0]
|
||||||
fCenter := freqs[idxPeak] // Center frequency.
|
|
||||||
|
|
||||||
// Create a window approximately 200 cents in width.
|
// TODO: At higher harmonic numbers, the harmonics are going to be spaced
|
||||||
deltaIdx := int(math.Round(deltaF / (100 * util.CentDeltaFreq(fCenter))))
|
// more cloesly together. That means that the window size must shrink.
|
||||||
idxMin := idxPeak - deltaIdx
|
winSizeCents := float64(150)
|
||||||
idxMax := idxPeak + deltaIdx
|
|
||||||
|
|
||||||
if idxMin < 0 {
|
// To start, we create a window around the approximate harmonic location.
|
||||||
idxMin = 0
|
idxPeak := int((fFundamental * harmonic) / deltaF)
|
||||||
|
winSize := int(
|
||||||
|
math.Round(
|
||||||
|
winSizeCents * util.CentDeltaFreq(fFundamental*harmonic) /
|
||||||
|
deltaF))
|
||||||
|
amps := rawAmps[idxPeak-winSize : idxPeak+winSize]
|
||||||
|
freqs := rawFreqs[idxPeak-winSize : idxPeak+winSize]
|
||||||
|
|
||||||
|
// Next, we refine the window be centering it on the actual peak.
|
||||||
|
idxPeak += (floats.MaxIdx(amps) - winSize)
|
||||||
|
amps = rawAmps[idxPeak-winSize : idxPeak+winSize]
|
||||||
|
freqs = rawFreqs[idxPeak-winSize : idxPeak+winSize]
|
||||||
|
idxPeak = winSize
|
||||||
|
|
||||||
|
// Next we compute statistical properties to determine if this peak is worth
|
||||||
|
// using.
|
||||||
|
mean, std := stat.MeanStdDev(amps, nil)
|
||||||
|
ampPeak := amps[idxPeak]
|
||||||
|
|
||||||
|
log.Printf("%0.2f: %0.4f", harmonic, (ampPeak-mean)/std)
|
||||||
|
|
||||||
|
return &HarmonicWindow{
|
||||||
|
Harmonic: harmonic,
|
||||||
|
Amps: amps,
|
||||||
|
Freqs: freqs,
|
||||||
|
IdxPeak: idxPeak,
|
||||||
|
AmpPeak: ampPeak,
|
||||||
|
AmpMean: mean,
|
||||||
|
AmpStd: std,
|
||||||
}
|
}
|
||||||
|
|
||||||
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
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -1,6 +1,7 @@
|
||||||
package extractfreq
|
package extractfreq
|
||||||
|
|
||||||
import (
|
import (
|
||||||
|
"fmt"
|
||||||
"log"
|
"log"
|
||||||
"math/cmplx"
|
"math/cmplx"
|
||||||
|
|
||||||
|
@ -62,28 +63,41 @@ func PlotFFT(path string, nSamples int, targetFreq float64) {
|
||||||
*/
|
*/
|
||||||
// Find frequency window around midi note.
|
// Find frequency window around midi note.
|
||||||
|
|
||||||
p, err := plot.New()
|
harmonics := []float64{0.25, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}
|
||||||
if err != nil {
|
for _, harmonic := range harmonics {
|
||||||
panic(err)
|
|
||||||
}
|
|
||||||
|
|
||||||
log.Printf("Making points...")
|
hw := NewHarmonicWindow(
|
||||||
pts := make(plotter.XYs, len(amps))
|
targetFreq,
|
||||||
for i := range amps {
|
harmonic,
|
||||||
pts[i].X = freqs[i]
|
amps,
|
||||||
pts[i].Y = amps[i]
|
freqs)
|
||||||
}
|
|
||||||
|
|
||||||
log.Printf("Plotting...")
|
p, err := plot.New()
|
||||||
if err = plotutil.AddLines(p, "fft", pts); err != nil {
|
if err != nil {
|
||||||
panic(err)
|
panic(err)
|
||||||
}
|
}
|
||||||
|
|
||||||
p.Y.Scale = plot.LogScale{}
|
log.Printf("Making points...")
|
||||||
//p.X.Scale = plot.LogScale{}
|
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("Saving...")
|
log.Printf("Plotting...")
|
||||||
if err := p.Save(8*vg.Inch, 8*vg.Inch, "plot.png"); err != nil {
|
if err = plotutil.AddLines(p, "fft", pts); err != nil {
|
||||||
panic(err)
|
panic(err)
|
||||||
|
}
|
||||||
|
|
||||||
|
//p.Y.Scale = plot.LogScale{}
|
||||||
|
//p.X.Scale = plot.LogScale{}
|
||||||
|
|
||||||
|
log.Printf("Saving...")
|
||||||
|
filename := fmt.Sprintf("plot-%03.2f.png", harmonic)
|
||||||
|
if err := p.Save(8*vg.Inch, 8*vg.Inch, filename); err != nil {
|
||||||
|
panic(err)
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
Reference in New Issue