From 80ab5e476877e5f673e25906c97727913cfeb2e9 Mon Sep 17 00:00:00 2001 From: ausocean-david Date: Wed, 14 Dec 2022 01:53:10 +1030 Subject: [PATCH 1/9] Audiofiltering: Create Lowpass filter with frequency control, with efficient convolution algorithm. --- cmd/audiofiltering/go.mod | 5 + cmd/audiofiltering/go.sum | 2 + cmd/audiofiltering/main.go | 225 +++++++++++++++++++++++++++++++++++++ 3 files changed, 232 insertions(+) create mode 100644 cmd/audiofiltering/go.mod create mode 100644 cmd/audiofiltering/go.sum create mode 100644 cmd/audiofiltering/main.go diff --git a/cmd/audiofiltering/go.mod b/cmd/audiofiltering/go.mod new file mode 100644 index 00000000..65df9827 --- /dev/null +++ b/cmd/audiofiltering/go.mod @@ -0,0 +1,5 @@ +module bitbucket.org/ausocean/av/cmd/audiofiltering + +go 1.19 + +require github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 diff --git a/cmd/audiofiltering/go.sum b/cmd/audiofiltering/go.sum new file mode 100644 index 00000000..c029ccc0 --- /dev/null +++ b/cmd/audiofiltering/go.sum @@ -0,0 +1,2 @@ +github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 h1:dd7vnTDfjtwCETZDrRe+GPYNLA1jBtbZeyfyE8eZCyk= +github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12/go.mod h1:i/KKcxEWEO8Yyl11DYafRPKOPVYTrhxiTRigjtEEXZU= diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go new file mode 100644 index 00000000..5fec888e --- /dev/null +++ b/cmd/audiofiltering/main.go @@ -0,0 +1,225 @@ +package main + +import( + "math" + "fmt" + "os" + "encoding/binary" + "github.com/mjibson/go-dsp/fft" + // "github.com/mjibson/go-dsp/window" + "math/cmplx" + "time" +) + +// define constants used in the generation of the sound waves +const( + SampleRate float64 = 44100 + Duration = 2 + tau = math.Pi * 2 + length int = 88200 +) + +func main() { + + start := time.Now() + + // generate two sine waves with different frequencies to test frequency response + n := 2 + audio := make([][]float64, n) + // for i:=0; i runMax[0]: + for i:=0; i<4; i++ { + runMax[4-i] = runMax[4-(1+i)] + indices[4-i] = indices[4-(1+i)] + } + runMax[0] = a[i] + indices[0] = i + case a[i] > runMax[1]: + for i:=0; i<3; i++ { + runMax[4-i] = runMax[4-(1+i)] + indices[4-i] = indices[4-(1+i)] + } + runMax[1] = a[i] + indices[1] = i + case a[i] > runMax[2]: + for i:=0; i<2; i++ { + runMax[4-i] = runMax[4-(1+i)] + indices[4-i] = indices[4-(1+i)] + } + runMax[2] = a[i] + indices[2] = i + case a[i] > runMax[3]: + for i:=0; i<1; i++ { + runMax[4-i] = runMax[4-(1+i)] + indices[4-i] = indices[4-(1+i)] + } + runMax[3] = a[i] + indices[3] = i + } + } + + return runMax, indices + +} + +func Max (a []float64) float64 { + + var runMax float64 = -1 + for i:= range a { + if math.Abs(a[i]) > runMax { + runMax = math.Abs(a[i]) + } + } + + return runMax + +} + +func LowPass (n int, fc float64) (filter []float64) { + + // n is number of points on either side of 0 + // determine digital frequency equivalent for fc + fd := fc/(2*SampleRate) + // create sinc function + return Sinc(n, fd) + + +} + +func Convolve (x, h []float64) []float64 { + + convLen := len(x)+len(h) + y := make([]float64, convLen) + for n:=0; n= 0 && n-k < len(h) { + sum += x[k]*h[n-k] + } + } + y[n] = sum +} + + +func SaveAudioData (signal []float64, fileName string) { + + // compute fft of signal + FFTaudio := fft.FFTReal(signal) + + // normalise and save signal + spectrum := make([]float64, len(signal)) + for i := range FFTaudio { + spectrum[i] = cmplx.Abs(FFTaudio[i])/float64(length) + } + maximum := Max(spectrum) + for i := range spectrum { + spectrum[i] = spectrum[i]/maximum + } + // spectrumAlt := spectrum[0:length/2 + 1] + + // SAVE + file := fileName + ".txt" + f, _ := os.Create(file) + for i:=0; i<20000; i++ { + fmt.Fprintf(f, "%v\n", /*10*math.Log10*/(spectrum[i])) + } + fmt.Printf("Saved spectrum values to: %s\n", fileName) + + file = fileName + ".bin" + f, _ = os.Create(file) + var buf [8]byte + for i:=0; i Date: Thu, 15 Dec 2022 00:42:44 +1030 Subject: [PATCH 2/9] Audiofiltering: Increase the efficieny and stability of algorithms by making use of waitgroups and goroutines. --- cmd/audiofiltering/main.go | 185 ++++++++++++++----------------------- 1 file changed, 70 insertions(+), 115 deletions(-) diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index 5fec888e..8a66ed22 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -6,9 +6,9 @@ import( "os" "encoding/binary" "github.com/mjibson/go-dsp/fft" - // "github.com/mjibson/go-dsp/window" "math/cmplx" "time" + "sync" ) // define constants used in the generation of the sound waves @@ -21,18 +21,27 @@ const( func main() { + var wg1, wg2 sync.WaitGroup start := time.Now() - // generate two sine waves with different frequencies to test frequency response - n := 2 + // generate sine waves with different frequencies to test frequency response + n := 100 + wg1.Add(n) audio := make([][]float64, n) - // for i:=0; i runMax[0]: - for i:=0; i<4; i++ { - runMax[4-i] = runMax[4-(1+i)] - indices[4-i] = indices[4-(1+i)] - } - runMax[0] = a[i] - indices[0] = i - case a[i] > runMax[1]: - for i:=0; i<3; i++ { - runMax[4-i] = runMax[4-(1+i)] - indices[4-i] = indices[4-(1+i)] - } - runMax[1] = a[i] - indices[1] = i - case a[i] > runMax[2]: - for i:=0; i<2; i++ { - runMax[4-i] = runMax[4-(1+i)] - indices[4-i] = indices[4-(1+i)] - } - runMax[2] = a[i] - indices[2] = i - case a[i] > runMax[3]: - for i:=0; i<1; i++ { - runMax[4-i] = runMax[4-(1+i)] - indices[4-i] = indices[4-(1+i)] - } - runMax[3] = a[i] - indices[3] = i - } - } - - return runMax, indices - -} - -func Max (a []float64) float64 { - - var runMax float64 = -1 - for i:= range a { - if math.Abs(a[i]) > runMax { - runMax = math.Abs(a[i]) - } - } - - return runMax - -} - -func LowPass (n int, fc float64) (filter []float64) { +func LowPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { // n is number of points on either side of 0 // determine digital frequency equivalent for fc - fd := fc/(2*SampleRate) + fd := (fc*4)/(2*SampleRate) // create sinc function - return Sinc(n, fd) - + ch <- Sinc(n, fd) + wg.Done() } -func Convolve (x, h []float64) []float64 { +func Convolve (x, h []float64, wg *sync.WaitGroup) []float64 { convLen := len(x)+len(h) y := make([]float64, convLen) for n:=0; n= 0 && n-k < len(h) { + sum += x[k]*h[n-k] + } + } + y[n] = sum + wg.Done() + }(n, x, h, y, wg) } - + wg.Done() return y } -func SubConvolve (n int, x, h, y []float64) { - var sum float64 = 0 - for k:=0; k= 0 && n-k < len(h) { - sum += x[k]*h[n-k] - } - } - y[n] = sum -} +func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { - -func SaveAudioData (signal []float64, fileName string) { - // compute fft of signal FFTaudio := fft.FFTReal(signal) @@ -188,23 +134,32 @@ func SaveAudioData (signal []float64, fileName string) { for i := range spectrum { spectrum[i] = spectrum[i]/maximum } - // spectrumAlt := spectrum[0:length/2 + 1] // SAVE - file := fileName + ".txt" - f, _ := os.Create(file) - for i:=0; i<20000; i++ { - fmt.Fprintf(f, "%v\n", /*10*math.Log10*/(spectrum[i])) - } - fmt.Printf("Saved spectrum values to: %s\n", fileName) + wg1.Add(2) + go func (fileName string, length int, spectrum []float64, wg1 *sync.WaitGroup) { + file := fileName + ".txt" + f, _ := os.Create(file) + for i:=0; i<20000; i++ { + fmt.Fprintf(f, "%v\n", /*10*math.Log10*/(spectrum[i])) + } + fmt.Printf("Saved spectrum values to: %s\n", fileName) + wg1.Done() + }(fileName, 44100, spectrum, wg1) + + go func (fileName string, signal []float64, wg1 *sync.WaitGroup) { + file := fileName + ".bin" + f, _ := os.Create(file) + var buf [8]byte + for i:=0; i Date: Thu, 15 Dec 2022 03:04:26 +1030 Subject: [PATCH 3/9] Audiofiltering: Add highpass filter functionality. --- cmd/audiofiltering/main.go | 38 ++++++++++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index 8a66ed22..c4d99939 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -36,10 +36,10 @@ func main() { } // create filter - filterLen := 500 + filterLen := 25 wg2.Add(1) ch := make(chan []float64, 1) - go LowPass(filterLen, 5000, &wg2, ch) + go HighPass(filterLen, 2000, &wg2, ch) wg1.Wait() // combine audio @@ -60,7 +60,7 @@ func main() { wg2.Add(1) filteredAudio := Convolve(combinedAudio, filter, &wg2) wg2.Wait() - go SaveAudioData(filteredAudio, "lowpass", &wg1) + go SaveAudioData(filteredAudio, "highpass", &wg1) wg1.Wait() fmt.Println(time.Since(start)) @@ -87,6 +87,19 @@ func generate(Frequency float64) []float64 { } +func Max (a []float64) float64 { + + var runMax float64 = -1 + for i:= range a { + if math.Abs(a[i]) > runMax { + runMax = math.Abs(a[i]) + } + } + + return runMax + +} + func LowPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { // n is number of points on either side of 0 @@ -98,6 +111,23 @@ func LowPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { } +func HighPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { + + // n is number of points on either side of 0 + // determine digital frequency equivalent for fc + fd := 0.5-((2*fc)/(SampleRate)) + fmt.Println(fd) + // create sinc function + filter := Sinc(n, fd) + // convert lowpass to highpass filter + for i:=-n; i<=n; i++ { + filter[i+n] = filter[i+n] * math.Pow(-1, float64(i)) + } + ch <- filter + wg.Done() + +} + func Convolve (x, h []float64, wg *sync.WaitGroup) []float64 { convLen := len(x)+len(h) @@ -141,7 +171,7 @@ func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { file := fileName + ".txt" f, _ := os.Create(file) for i:=0; i<20000; i++ { - fmt.Fprintf(f, "%v\n", /*10*math.Log10*/(spectrum[i])) + fmt.Fprintf(f, "%v\n", /*20*math.Log10*/(spectrum[i])) } fmt.Printf("Saved spectrum values to: %s\n", fileName) wg1.Done() From 49db1041d3d23b5bac74a735f6a7e225d72cbbd1 Mon Sep 17 00:00:00 2001 From: ausocean-david Date: Sat, 17 Dec 2022 02:33:57 +1030 Subject: [PATCH 4/9] Audiofiltering: Improve filter performance by reducing frequency leakage. --- cmd/audiofiltering/main.go | 114 +++++++++++++++++++++++-------------- 1 file changed, 70 insertions(+), 44 deletions(-) diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index c4d99939..c8a17d28 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -6,6 +6,7 @@ import( "os" "encoding/binary" "github.com/mjibson/go-dsp/fft" + "github.com/mjibson/go-dsp/window" "math/cmplx" "time" "sync" @@ -36,7 +37,7 @@ func main() { } // create filter - filterLen := 25 + filterLen := 500 wg2.Add(1) ch := make(chan []float64, 1) go HighPass(filterLen, 2000, &wg2, ch) @@ -55,12 +56,12 @@ func main() { filter := <- ch // convolve sinc with audio - wg1.Add(2) - go SaveAudioData(combinedAudio, "unfiltered", &wg1) + wg1.Add(1) + // go SaveAudioData(combinedAudio, "unfiltered", &wg1) wg2.Add(1) filteredAudio := Convolve(combinedAudio, filter, &wg2) wg2.Wait() - go SaveAudioData(filteredAudio, "highpass", &wg1) + go SaveAudioData(filteredAudio, "highpass-2KHz", &wg1) wg1.Wait() fmt.Println(time.Since(start)) @@ -100,29 +101,41 @@ func Max (a []float64) float64 { } -func LowPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { +func LowPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { - // n is number of points on either side of 0 - // determine digital frequency equivalent for fc - fd := (fc*4)/(2*SampleRate) // create sinc function - ch <- Sinc(n, fd) + size := taps + 1 + fd := fc/SampleRate + b := (2*math.Pi) * fd + filter := make([]float64, size) + winData := window.FlatTop(size) + for n:=0; n<(taps/2); n++ { + c := float64(n) - float64(taps)/2 + y := math.Sin(c*b) / (math.Pi * c) + filter[n] = y * winData[n] + filter[size-1-n] = filter[n] + } + filter[taps/2] = 2*fd*winData[taps/2] + ch <- filter wg.Done() } -func HighPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { +func HighPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { - // n is number of points on either side of 0 - // determine digital frequency equivalent for fc - fd := 0.5-((2*fc)/(SampleRate)) - fmt.Println(fd) // create sinc function - filter := Sinc(n, fd) - // convert lowpass to highpass filter - for i:=-n; i<=n; i++ { - filter[i+n] = filter[i+n] * math.Pow(-1, float64(i)) + size := taps + 1 + fd := fc/SampleRate + b := (2*math.Pi) * fd + filter := make([]float64, size) + winData := window.FlatTop(size) + for n:=0; n<(taps/2); n++ { + c := float64(n) - float64(taps)/2 + y := math.Sin(c*b) / (math.Pi * c) + filter[n] = -y * winData[n] + filter[size-1-n] = filter[n] } + filter[taps/2] = (1 - 2*fd) * winData[taps/2] ch <- filter wg.Done() @@ -130,20 +143,21 @@ func HighPass (n int, fc float64, wg *sync.WaitGroup, ch chan []float64) { func Convolve (x, h []float64, wg *sync.WaitGroup) []float64 { - convLen := len(x)+len(h) + // conv waitgroup + // var convwg sync.WaitGroup + convLen := len(x)+len(h)-1 y := make([]float64, convLen) for n:=0; n= 0 && n-k < len(h) { - sum += x[k]*h[n-k] + go func(n int, y []float64) { + var sum float64 = 0 + for k:=0; k= 0 && n-k < len(h) { + sum += x[k]*h[n-k] + } } - } y[n] = sum - wg.Done() - }(n, x, h, y, wg) + } (n, y) + } wg.Done() return y @@ -160,6 +174,10 @@ func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { for i := range FFTaudio { spectrum[i] = cmplx.Abs(FFTaudio[i])/float64(length) } + // maximum := Max(signal) + // for i := range signal { + // signal[i] = signal[i]/maximum + // } maximum := Max(spectrum) for i := range spectrum { spectrum[i] = spectrum[i]/maximum @@ -170,8 +188,8 @@ func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { go func (fileName string, length int, spectrum []float64, wg1 *sync.WaitGroup) { file := fileName + ".txt" f, _ := os.Create(file) - for i:=0; i<20000; i++ { - fmt.Fprintf(f, "%v\n", /*20*math.Log10*/(spectrum[i])) + for i:=0; i Date: Sat, 17 Dec 2022 03:27:57 +1030 Subject: [PATCH 5/9] Audiofiltering: Implement bandpass filter by combining lowpass and highpass filters. --- cmd/audiofiltering/main.go | 133 ++++++++++++++++++++++--------------- 1 file changed, 78 insertions(+), 55 deletions(-) diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index c8a17d28..fb56e774 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -12,7 +12,7 @@ import( "sync" ) -// define constants used in the generation of the sound waves +// Define the constants to be used in the generation of the sound files. const( SampleRate float64 = 44100 Duration = 2 @@ -20,12 +20,17 @@ const( length int = 88200 ) +// main is a driver function which generates a sound file which contains sine waves from across the spectrum. +// Filters can then be generated and applied to the audio, in order to test their frequency response. func main() { + // Create waitgroups to ensure program runs correctly. var wg1, wg2 sync.WaitGroup + + // Take a measurement of the starting point of the program execution to use for execution timing. start := time.Now() - // generate sine waves with different frequencies to test frequency response + // Generate sine waves with different frequencies to test frequency response. n := 100 wg1.Add(n) audio := make([][]float64, n) @@ -36,14 +41,15 @@ func main() { } (i, audio, &wg1) } - // create filter + // Create the filter. filterLen := 500 wg2.Add(1) ch := make(chan []float64, 1) - go HighPass(filterLen, 2000, &wg2, ch) + go BandPass(filterLen, 2000, 5000, &wg2, ch) + + // Combine all the generated sine waves into a single audio slice. wg1.Wait() - // combine audio combinedAudio := make([]float64, length) for i := range audio[0] { combinedAudio[i] = 0 @@ -52,34 +58,41 @@ func main() { } } + // Get the filter from the filter goroutine called earlier. wg2.Wait() filter := <- ch - // convolve sinc with audio - wg1.Add(1) - // go SaveAudioData(combinedAudio, "unfiltered", &wg1) + wg1.Add(2) + // Save the unfiltered audio to compare the filtered audio against. + go SaveAudioData(combinedAudio, "unfiltered", &wg1) + + // Apply the filter to the combined audio file by computing the convolution of the 2 slices. wg2.Add(1) filteredAudio := Convolve(combinedAudio, filter, &wg2) - wg2.Wait() - go SaveAudioData(filteredAudio, "highpass-2KHz", &wg1) + // Save the filtered audio. + wg2.Wait() + go SaveAudioData(filteredAudio, "bandpass-2-5KHz", &wg1) + + // Print the total time of execution, to compare different algorithms. wg1.Wait() fmt.Println(time.Since(start)) } +// generate is used to generate a sine wave of the given frequency, calls on the global variables of Duration and SampleRate func generate(Frequency float64) []float64 { - // deteremine number of samples based off duration and sample rate + // Deteremine the number of samples required based off Duration and SampleRate. nsamps := Duration * SampleRate - // generate x-values + // Generate the x-values to plot against. var angle float64 = tau / float64(nsamps) - // create sample array + // Create an array for the generated samples. samp := make([]float64, int(nsamps)) - // generate samples and write to file + // Generate the samples. for i := 0; i < int(nsamps); i++ { samp[i] = math.Sin(angle * Frequency * float64(4*i)) } @@ -88,7 +101,8 @@ func generate(Frequency float64) []float64 { } -func Max (a []float64) float64 { +// max returns the absolute highest value in a given array. +func max(a []float64) float64 { var runMax float64 = -1 for i:= range a { @@ -101,9 +115,10 @@ func Max (a []float64) float64 { } +// LowPass returns an array of coeffiecients over the given channel that can be used in a convolution with a signal to perform +// a lowpass filtering with the given cut-off frequency fc. func LowPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { - // create sinc function size := taps + 1 fd := fc/SampleRate b := (2*math.Pi) * fd @@ -121,9 +136,10 @@ func LowPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { } +// HighPass returns an array of coeffiecients over the given channel that can be used in a convolution with a signal to perform +// a highpass filtering with the given cut-off frequency fc. func HighPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { - // create sinc function size := taps + 1 fd := fc/SampleRate b := (2*math.Pi) * fd @@ -141,14 +157,47 @@ func HighPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { } +// BandPass returns an array of coeffiecients over the given channel that can be used in a convolution with a signal to perform +// a bandpass filtering with the given lower and upper cutoff frequencies (using the low and highpass filter functions). +func BandPass (taps int, lower, upper float64, wg *sync.WaitGroup, ch chan []float64) { + + // Create a waitgroup and channel to be used internally within the function. + var bpwg sync.WaitGroup + bpch := make(chan []float64, 2) + + // Call low and highpass filter functions. + filters := make([][]float64, 2) + bpwg.Add(2) + go LowPass(taps, upper, &bpwg, bpch) + go HighPass(taps, lower, &bpwg, bpch) + bpwg.Wait() + for i:=0; i<2; i++ { + filters[i] = <- bpch + } + + // Convolve lowpass filter with highpass filter to get bandpass filter. + bpwg.Add(1) + bpFilter := Convolve(filters[0], filters[1], &bpwg) + ch <- bpFilter + + wg.Done() + +} + + + +// Convolve returns the real convolution of the 2 given arrays. func Convolve (x, h []float64, wg *sync.WaitGroup) []float64 { - // conv waitgroup - // var convwg sync.WaitGroup + // Create a waitgroup to be used in goroutines called by Convolution + var convwg sync.WaitGroup + + // Compute the convolution convLen := len(x)+len(h)-1 y := make([]float64, convLen) for n:=0; n= 0 && n-k < len(h) { @@ -156,34 +205,31 @@ func Convolve (x, h []float64, wg *sync.WaitGroup) []float64 { } } y[n] = sum - } (n, y) - + wg.Done() + } (n, y, &convwg) } wg.Done() return y } +// SaveAudioData takes an input signal and filename and saves the frequency spectrum (to .txt) and the audio signal (to .bin). func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { - // compute fft of signal + // Compute the FFT of the signal FFTaudio := fft.FFTReal(signal) - // normalise and save signal + // Normalise the frequency response (-1:1). spectrum := make([]float64, len(signal)) for i := range FFTaudio { spectrum[i] = cmplx.Abs(FFTaudio[i])/float64(length) } - // maximum := Max(signal) - // for i := range signal { - // signal[i] = signal[i]/maximum - // } - maximum := Max(spectrum) + maximum := max(spectrum) for i := range spectrum { spectrum[i] = spectrum[i]/maximum } - // SAVE + // Save the frequency response/spectrum. wg1.Add(2) go func (fileName string, length int, spectrum []float64, wg1 *sync.WaitGroup) { file := fileName + ".txt" @@ -195,6 +241,7 @@ func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { wg1.Done() }(fileName, 44100, spectrum, wg1) + // Save the audio file. go func (fileName string, signal []float64, wg1 *sync.WaitGroup) { file := fileName + ".bin" f, _ := os.Create(file) @@ -209,28 +256,4 @@ func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { wg1.Done() -} - -// func SincRedundant (nsamps int, fc float64) []float64 { - -// // n is number of samples either side of n = 0 -// h := make([]float64, nsamps*2 + 1) -// for n:=-nsamps; n<=nsamps; n++ { -// if n == 0 { -// h[n+nsamps] = 2*fc -// } else { -// h[n+nsamps] = math.Sin(2*math.Pi*fc*float64(n))/(math.Pi*float64(n)) -// } -// } - -// return h - -// } - -// func Sinc (x float64) float64 { -// if x == 0 { -// return -// } else { -// h[n+nsamps] = math.Sin(2*math.Pi*fc*float64(n))/(math.Pi*float64(n)) -// } -// } \ No newline at end of file +} \ No newline at end of file From d029038db9db3ec599b210bb78671453c005f970 Mon Sep 17 00:00:00 2001 From: ausocean-david Date: Mon, 19 Dec 2022 20:46:20 +1030 Subject: [PATCH 6/9] Audiofiltering: Interface the function with pre-existing data structures. The filters can now be generated into a filter type using a generic Generate function. This generates a filter off of the specifications within the filter struct. There is a generic Apply function which takes in a buffer of PCM data (defined in pcm.go), and outputs to a []byte. --- cmd/audiofiltering/go.mod | 5 - cmd/audiofiltering/go.sum | 2 - cmd/audiofiltering/main.go | 273 +++++-------------------------------- codec/pcm/filters.go | 198 +++++++++++++++++++++++++++ codec/pcm/pcm.go | 6 +- go.mod | 5 +- go.sum | 32 ++--- 7 files changed, 249 insertions(+), 272 deletions(-) delete mode 100644 cmd/audiofiltering/go.mod delete mode 100644 cmd/audiofiltering/go.sum create mode 100644 codec/pcm/filters.go diff --git a/cmd/audiofiltering/go.mod b/cmd/audiofiltering/go.mod deleted file mode 100644 index 65df9827..00000000 --- a/cmd/audiofiltering/go.mod +++ /dev/null @@ -1,5 +0,0 @@ -module bitbucket.org/ausocean/av/cmd/audiofiltering - -go 1.19 - -require github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 diff --git a/cmd/audiofiltering/go.sum b/cmd/audiofiltering/go.sum deleted file mode 100644 index c029ccc0..00000000 --- a/cmd/audiofiltering/go.sum +++ /dev/null @@ -1,2 +0,0 @@ -github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 h1:dd7vnTDfjtwCETZDrRe+GPYNLA1jBtbZeyfyE8eZCyk= -github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12/go.mod h1:i/KKcxEWEO8Yyl11DYafRPKOPVYTrhxiTRigjtEEXZU= diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index fb56e774..5ff73afe 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -1,259 +1,56 @@ package main -import( - "math" +import ( "fmt" + "math" "os" - "encoding/binary" - "github.com/mjibson/go-dsp/fft" - "github.com/mjibson/go-dsp/window" - "math/cmplx" "time" - "sync" + + "bitbucket.org/ausocean/av/codec/pcm" ) -// Define the constants to be used in the generation of the sound files. -const( - SampleRate float64 = 44100 - Duration = 2 - tau = math.Pi * 2 - length int = 88200 -) - -// main is a driver function which generates a sound file which contains sine waves from across the spectrum. -// Filters can then be generated and applied to the audio, in order to test their frequency response. +// main is a driver function for testing the filters defined in codec/pcm/filters.go func main() { - // Create waitgroups to ensure program runs correctly. - var wg1, wg2 sync.WaitGroup - - // Take a measurement of the starting point of the program execution to use for execution timing. + // Define start time for execution timing. start := time.Now() - - // Generate sine waves with different frequencies to test frequency response. - n := 100 - wg1.Add(n) - audio := make([][]float64, n) - for i:=0; i runMax { - runMax = math.Abs(a[i]) - } - } - - return runMax - -} - -// LowPass returns an array of coeffiecients over the given channel that can be used in a convolution with a signal to perform -// a lowpass filtering with the given cut-off frequency fc. -func LowPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { - - size := taps + 1 - fd := fc/SampleRate - b := (2*math.Pi) * fd - filter := make([]float64, size) - winData := window.FlatTop(size) - for n:=0; n<(taps/2); n++ { - c := float64(n) - float64(taps)/2 - y := math.Sin(c*b) / (math.Pi * c) - filter[n] = y * winData[n] - filter[size-1-n] = filter[n] - } - filter[taps/2] = 2*fd*winData[taps/2] - ch <- filter - wg.Done() - -} - -// HighPass returns an array of coeffiecients over the given channel that can be used in a convolution with a signal to perform -// a highpass filtering with the given cut-off frequency fc. -func HighPass (taps int, fc float64, wg *sync.WaitGroup, ch chan []float64) { - - size := taps + 1 - fd := fc/SampleRate - b := (2*math.Pi) * fd - filter := make([]float64, size) - winData := window.FlatTop(size) - for n:=0; n<(taps/2); n++ { - c := float64(n) - float64(taps)/2 - y := math.Sin(c*b) / (math.Pi * c) - filter[n] = -y * winData[n] - filter[size-1-n] = filter[n] - } - filter[taps/2] = (1 - 2*fd) * winData[taps/2] - ch <- filter - wg.Done() - -} - -// BandPass returns an array of coeffiecients over the given channel that can be used in a convolution with a signal to perform -// a bandpass filtering with the given lower and upper cutoff frequencies (using the low and highpass filter functions). -func BandPass (taps int, lower, upper float64, wg *sync.WaitGroup, ch chan []float64) { - - // Create a waitgroup and channel to be used internally within the function. - var bpwg sync.WaitGroup - bpch := make(chan []float64, 2) - - // Call low and highpass filter functions. - filters := make([][]float64, 2) - bpwg.Add(2) - go LowPass(taps, upper, &bpwg, bpch) - go HighPass(taps, lower, &bpwg, bpch) - bpwg.Wait() - for i:=0; i<2; i++ { - filters[i] = <- bpch - } - - // Convolve lowpass filter with highpass filter to get bandpass filter. - bpwg.Add(1) - bpFilter := Convolve(filters[0], filters[1], &bpwg) - ch <- bpFilter - - wg.Done() - -} - - - -// Convolve returns the real convolution of the 2 given arrays. -func Convolve (x, h []float64, wg *sync.WaitGroup) []float64 { - - // Create a waitgroup to be used in goroutines called by Convolution - var convwg sync.WaitGroup - - // Compute the convolution - convLen := len(x)+len(h)-1 - y := make([]float64, convLen) - for n:=0; n= 0 && n-k < len(h) { - sum += x[k]*h[n-k] - } - } - y[n] = sum - wg.Done() - } (n, y, &convwg) - } - wg.Done() - return y - -} - -// SaveAudioData takes an input signal and filename and saves the frequency spectrum (to .txt) and the audio signal (to .bin). -func SaveAudioData (signal []float64, fileName string, wg1 *sync.WaitGroup) { - - // Compute the FFT of the signal - FFTaudio := fft.FFTReal(signal) - - // Normalise the frequency response (-1:1). - spectrum := make([]float64, len(signal)) - for i := range FFTaudio { - spectrum[i] = cmplx.Abs(FFTaudio[i])/float64(length) - } - maximum := max(spectrum) - for i := range spectrum { - spectrum[i] = spectrum[i]/maximum - } - - // Save the frequency response/spectrum. - wg1.Add(2) - go func (fileName string, length int, spectrum []float64, wg1 *sync.WaitGroup) { - file := fileName + ".txt" - f, _ := os.Create(file) - for i:=0; i 32767 { + inputAsFloat[i] -= 32768 * 2 + } + inputAsFloat[i] /= (32767) + } + + // Convolve input with filter. + var wg sync.WaitGroup + ch := make(chan []float64, 1) + wg.Add(1) + go Convolve(inputAsFloat, filter.Coeffs, &wg, ch) + wg.Wait() + convolution := <-ch + + // Convert convolution output back to bytes. + var output []byte + buf := make([]byte, 2) + for i := range convolution { + convolution[i] = convolution[i] * 32767 + if convolution[i] < 0 { + convolution[i] = convolution[i] + 32768*2 + } + binary.LittleEndian.PutUint16(buf[:], uint16(convolution[i])) + output = append(output, buf[0], buf[1]) + } + + return output + +} + +func Convolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { + + // Create a waitgroup to be used in goroutines called by Convolution + var convwg sync.WaitGroup + + // Compute the convolution + convLen := len(x) + len(h) - 1 + y := make([]float64, convLen) + for n := 0; n < convLen; n++ { + convwg.Add(1) + go func(n int, y []float64, convwg *sync.WaitGroup) { + var sum float64 = 0 + for k := 0; k < len(x); k++ { + if n-k >= 0 && n-k < len(h) { + sum += x[k] * h[n-k] + } else if n-k < 0 { + break + } + } + y[n] = sum + convwg.Done() + }(n, y, &convwg) + } + ch <- y + close(ch) + wg.Done() +} + +// LEGACY FUNCTION FOR TESTING +// Max returns the absolute highest value in a given array. +func Max(a []float64) float64 { + + var runMax float64 = -1 + for i := range a { + if math.Abs(a[i]) > runMax { + runMax = math.Abs(a[i]) + } + } + + return runMax + +} diff --git a/codec/pcm/pcm.go b/codec/pcm/pcm.go index eb93bce7..dda88e93 100644 --- a/codec/pcm/pcm.go +++ b/codec/pcm/pcm.go @@ -73,9 +73,9 @@ func DataSize(rate, channels, bitDepth uint, period float64) int { // Resample takes Buffer c and resamples the pcm audio data to 'rate' Hz and returns a Buffer with the resampled data. // Notes: -// - Currently only downsampling is implemented and c's rate must be divisible by 'rate' or an error will occur. -// - If the number of bytes in c.Data is not divisible by the decimation factor (ratioFrom), the remaining bytes will -// not be included in the result. Eg. input of length 480002 downsampling 6:1 will result in output length 80000. +// - Currently only downsampling is implemented and c's rate must be divisible by 'rate' or an error will occur. +// - If the number of bytes in c.Data is not divisible by the decimation factor (ratioFrom), the remaining bytes will +// not be included in the result. Eg. input of length 480002 downsampling 6:1 will result in output length 80000. func Resample(c Buffer, rate uint) (Buffer, error) { if c.Format.Rate == rate { return c, nil diff --git a/go.mod b/go.mod index 8ba7ec85..7a7b1257 100644 --- a/go.mod +++ b/go.mod @@ -12,10 +12,11 @@ require ( github.com/google/go-cmp v0.4.1 github.com/kidoman/embd v0.0.0-20170508013040-d3d8c0c5c68d github.com/mewkiz/flac v1.0.5 + github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 github.com/pkg/errors v0.9.1 github.com/yobert/alsa v0.0.0-20180630182551-d38d89fa843e gocv.io/x/gocv v0.29.0 - gonum.org/v1/gonum v0.9.3 - gonum.org/v1/plot v0.10.0 + gonum.org/v1/gonum v0.8.2 + gonum.org/v1/plot v0.9.0 gopkg.in/natefinch/lumberjack.v2 v2.0.0 ) diff --git a/go.sum b/go.sum index d340b280..cc2a47a7 100644 --- a/go.sum +++ b/go.sum @@ -16,13 +16,11 @@ github.com/Shopify/toxiproxy v2.1.4+incompatible h1:TKdv8HiTLgE5wdJuEML90aBgNWso github.com/Shopify/toxiproxy v2.1.4+incompatible/go.mod h1:OXgGpZ6Cli1/URJOF1DMxUHB2q5Ap20/P/eIdh4G0pI= github.com/adrianmo/go-nmea v1.1.1-0.20190109062325-c448653979f7/go.mod h1:HHPxPAm2kmev+61qmkZh7xgZF/7qHtSpsWppip2Ipv8= github.com/aerth/playwav v0.0.0-20170324024803-17dfe21a406f/go.mod h1:aTANbm/GXj1ilCRMsIiSpX0i7LUkqjAPj6R0fpWbLNA= +github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af h1:wVe6/Ea46ZMeNkQjjBW6xcqyQA/j5e0D6GytH95g0gQ= github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af/go.mod h1:K08gAheRH3/J6wwsYMMT4xOr94bZjxIelGM0+d/wbFw= -github.com/ajstarks/svgo v0.0.0-20210923152817-c3b6e2f0c527 h1:NImof/JkF93OVWZY+PINgl6fPtQyF6f+hNUtZ0QZA1c= -github.com/ajstarks/svgo v0.0.0-20210923152817-c3b6e2f0c527/go.mod h1:K08gAheRH3/J6wwsYMMT4xOr94bZjxIelGM0+d/wbFw= github.com/andreyvit/diff v0.0.0-20170406064948-c7f18ee00883 h1:bvNMNQO63//z+xNgfBlViaCIJKLlCJ6/fmUseuG0wVQ= github.com/andreyvit/diff v0.0.0-20170406064948-c7f18ee00883/go.mod h1:rCTlJbsFo29Kk6CurOXKm700vrz8f0KW0JNfpkRJY/8= github.com/boombuler/barcode v1.0.0/go.mod h1:paBWMcWSl3LHKBqUq+rly7CNSldXjb2rDl3JlRe0mD8= -github.com/boombuler/barcode v1.0.1/go.mod h1:paBWMcWSl3LHKBqUq+rly7CNSldXjb2rDl3JlRe0mD8= github.com/cheekybits/is v0.0.0-20150225183255-68e9c0620927/go.mod h1:h/aW8ynjgkuj+NQRlZcDbAbM1ORAbXjXX77sX7T289U= github.com/cocoonlife/goalsa v0.0.0-20160812085113-b711ae6f3eff/go.mod h1:5aLO409bJnd+jCw0t/SB/DhHkVBhPAE31lnHJnYQxy0= github.com/cocoonlife/testify v0.0.0-20160218172820-792cc1faeb64/go.mod h1:LoCAz53rbPcqs8Da2BjB/yDy4gxMtiSQmqnYI/DGH+U= @@ -46,18 +44,13 @@ github.com/go-audio/wav v0.0.0-20181013172942-de841e69b884 h1:2TaXIaVA4ff/MHHezO github.com/go-audio/wav v0.0.0-20181013172942-de841e69b884/go.mod h1:UiqzUyfX0zs3pJ/DPyvS5v8sN6s5bXPUDDIVA5v8dks= github.com/go-fonts/dejavu v0.1.0 h1:JSajPXURYqpr+Cu8U9bt8K+XcACIHWqWrvWCKyeFmVQ= github.com/go-fonts/dejavu v0.1.0/go.mod h1:4Wt4I4OU2Nq9asgDCteaAaWZOV24E+0/Pwo0gppep4g= -github.com/go-fonts/latin-modern v0.2.0 h1:5/Tv1Ek/QCr20C6ZOz15vw3g7GELYL98KWr8Hgo+3vk= github.com/go-fonts/latin-modern v0.2.0/go.mod h1:rQVLdDMK+mK1xscDwsqM5J8U2jrRa3T0ecnM9pNujks= +github.com/go-fonts/liberation v0.1.1 h1:wBrPaMkrXFBW3qXpXAjiKljdVUMxn9bX2ia3XjPHoik= github.com/go-fonts/liberation v0.1.1/go.mod h1:K6qoJYypsmfVjWg8KOVDQhLc8UDgIK2HYqyqAO9z7GY= -github.com/go-fonts/liberation v0.2.0 h1:jAkAWJP4S+OsrPLZM4/eC9iW7CtHy+HBXrEwZXWo5VM= -github.com/go-fonts/liberation v0.2.0/go.mod h1:K6qoJYypsmfVjWg8KOVDQhLc8UDgIK2HYqyqAO9z7GY= github.com/go-fonts/stix v0.1.0/go.mod h1:w/c1f0ldAUlJmLBvlbkvVXLAD+tAMqobIIQpmnUIzUY= github.com/go-gl/glfw v0.0.0-20190409004039-e6da0acd62b1/go.mod h1:vR7hzQXu2zJy9AVAgeJqvqgH9Q5CA+iKCZ2gyEVpxRU= +github.com/go-latex/latex v0.0.0-20210118124228-b3d85cf34e07 h1:OTlfMvwR1rLyf9goVmXfuS5AJn80+Vmj4rTf4n46SOs= github.com/go-latex/latex v0.0.0-20210118124228-b3d85cf34e07/go.mod h1:CO1AlKB2CSIqUrmQPqA0gdRIlnLEY0gK5JGjh37zN5U= -github.com/go-latex/latex v0.0.0-20210823091927-c0d11ff05a81 h1:6zl3BbBhdnMkpSj2YY30qV3gDcVBGtFgVsV3+/i+mKQ= -github.com/go-latex/latex v0.0.0-20210823091927-c0d11ff05a81/go.mod h1:SX0U8uGpxhq9o2S/CELCSUxEWWAuoCUcVCQWv7G2OCk= -github.com/go-pdf/fpdf v0.5.0 h1:GHpcYsiDV2hdo77VTOuTF9k1sN8F8IY7NjnCo9x+NPY= -github.com/go-pdf/fpdf v0.5.0/go.mod h1:HzcnA+A23uwogo0tp9yU+l3V+KXhiESpt1PMayhOh5M= github.com/go-ping/ping v0.0.0-20201115131931-3300c582a663/go.mod h1:35JbSyV/BYqHwwRA6Zr1uVDm1637YlNOU61wI797NPI= github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 h1:DACJavvAHhabrF08vX0COfcOBJRhZ8lUbR+ZWIs0Y5g= github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0/go.mod h1:E/TSTwGwJL78qG/PmXZO1EjYhfJinVAhrmmHX6Z8B9k= @@ -79,17 +72,18 @@ github.com/mattetti/audio v0.0.0-20180912171649-01576cde1f21 h1:Hc1iKlyxNHp3CV59 github.com/mattetti/audio v0.0.0-20180912171649-01576cde1f21/go.mod h1:LlQmBGkOuV/SKzEDXBPKauvN2UqCgzXO2XjecTGj40s= github.com/mewkiz/flac v1.0.5 h1:dHGW/2kf+/KZ2GGqSVayNEhL9pluKn/rr/h/QqD9Ogc= github.com/mewkiz/flac v1.0.5/go.mod h1:EHZNU32dMF6alpurYyKHDLYpW1lYpBZ5WrXi/VuNIGs= +github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 h1:dd7vnTDfjtwCETZDrRe+GPYNLA1jBtbZeyfyE8eZCyk= +github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12/go.mod h1:i/KKcxEWEO8Yyl11DYafRPKOPVYTrhxiTRigjtEEXZU= github.com/pascaldekloe/goe v0.1.0/go.mod h1:lzWF7FIEvWOWxwDKqyGYQf6ZUaNfKdP144TG7ZOy1lc= +github.com/phpdave11/gofpdf v1.4.2 h1:KPKiIbfwbvC/wOncwhrpRdXVj2CZTCFlw4wnoyjtHfQ= github.com/phpdave11/gofpdf v1.4.2/go.mod h1:zpO6xFn9yxo3YLyMvW8HcKWVdbNqgIfOOp2dXMnm1mY= github.com/phpdave11/gofpdi v1.0.12/go.mod h1:vBmVV0Do6hSBHC8uKUQ71JGW+ZGQq74llk/7bXwjDoI= -github.com/phpdave11/gofpdi v1.0.13/go.mod h1:vBmVV0Do6hSBHC8uKUQ71JGW+ZGQq74llk/7bXwjDoI= github.com/pkg/errors v0.8.1/go.mod h1:bwawxfHBFNV+L2hUp1rHADufV3IMtnDRdf1r5NINEl0= github.com/pkg/errors v0.9.1 h1:FEBLx1zS214owpjy7qsBeixbURkuhQAwrK5UwLGTwt4= github.com/pkg/errors v0.9.1/go.mod h1:bwawxfHBFNV+L2hUp1rHADufV3IMtnDRdf1r5NINEl0= github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM= github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4= github.com/ruudk/golang-pdf417 v0.0.0-20181029194003-1af4ab5afa58/go.mod h1:6lfFZQK844Gfx8o5WFuvpxWRwnSoipWe/p622j1v06w= -github.com/ruudk/golang-pdf417 v0.0.0-20201230142125-a7e3863a1245/go.mod h1:pQAZKsJ8yyVxGRWYNEm9oFB8ieLgKFnamEyDmSA0BRk= github.com/sergi/go-diff v1.0.0 h1:Kpca3qRNrduNnOQeazBd0ysaKrUJiIuISHxogkT9RPQ= github.com/sergi/go-diff v1.0.0/go.mod h1:0CfEIISq7TuYL3j771MWULgwwjU+GofnZX9QAmXWZgo= github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= @@ -129,10 +123,8 @@ golang.org/x/image v0.0.0-20200119044424-58c23975cae1/go.mod h1:FeLwcggjj3mMvU+o golang.org/x/image v0.0.0-20200430140353-33d19683fad8/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0= golang.org/x/image v0.0.0-20200618115811-c13761719519/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0= golang.org/x/image v0.0.0-20201208152932-35266b937fa6/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0= +golang.org/x/image v0.0.0-20210216034530-4410531fe030 h1:lP9pYkih3DUSC641giIXa2XqfTIbbbRr0w2EOTA7wHA= golang.org/x/image v0.0.0-20210216034530-4410531fe030/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0= -golang.org/x/image v0.0.0-20210607152325-775e3b0c77b9/go.mod h1:023OzeP/+EPmXeapQh35lcL3II3LrY8Ic+EFFKVhULM= -golang.org/x/image v0.0.0-20210628002857-a66eb6448b8d h1:RNPAfi2nHY7C2srAV8A49jpsYr0ADedCk1wq6fTMTvs= -golang.org/x/image v0.0.0-20210628002857-a66eb6448b8d/go.mod h1:023OzeP/+EPmXeapQh35lcL3II3LrY8Ic+EFFKVhULM= golang.org/x/mobile v0.0.0-20190719004257-d2bd2a29d028/go.mod h1:E/iHnbuqvinMTCcRqshq8CkpyQDoeVncDDYHnLhea+o= golang.org/x/mod v0.1.0/go.mod h1:0QHyrYULN0/3qlju5TqG8bIK38QM8yzMo5ekMj3DlcY= golang.org/x/net v0.0.0-20190404232315-eb5bcb51f2a3/go.mod h1:t9HGtf8HONx5eT2rtn7q6eTqICYqUVnKs3thJo3Qplg= @@ -152,9 +144,8 @@ golang.org/x/sys v0.0.0-20210909193231-528a39cd75f3/go.mod h1:oPkhp1MJrh7nUepCBc golang.org/x/term v0.0.0-20201126162022-7de9c90e9dd1/go.mod h1:bj7SfCRtBDWHUb9snDiAeCFNEtKQo2Wmx5Cou7ajbmo= golang.org/x/text v0.3.0/go.mod h1:NqM8EUOU14njkJ3fqMW+pc6Ldnwhi/IjpwHt7yyuwOQ= golang.org/x/text v0.3.3/go.mod h1:5Zoc/QRtKVWzQhOtBMvqHzDpF6irO9z98xDceosuGiQ= +golang.org/x/text v0.3.5 h1:i6eZZ+zk0SOf0xgBpEpPD18qWcJda6q1sxt3S0kzyUQ= golang.org/x/text v0.3.5/go.mod h1:5Zoc/QRtKVWzQhOtBMvqHzDpF6irO9z98xDceosuGiQ= -golang.org/x/text v0.3.6 h1:aRYxNxv6iGQlyVaZmk6ZgYEDa+Jg18DxebPSrd6bg1M= -golang.org/x/text v0.3.6/go.mod h1:5Zoc/QRtKVWzQhOtBMvqHzDpF6irO9z98xDceosuGiQ= golang.org/x/tools v0.0.0-20180525024113-a5b4c53f6e8b/go.mod h1:n7NCudcB/nEzxVGmLbDWY5pfWTLqBcC2KZ6jyYvM4mQ= golang.org/x/tools v0.0.0-20180917221912-90fa682c2a6e/go.mod h1:n7NCudcB/nEzxVGmLbDWY5pfWTLqBcC2KZ6jyYvM4mQ= golang.org/x/tools v0.0.0-20190206041539-40960b6deb8e/go.mod h1:n7NCudcB/nEzxVGmLbDWY5pfWTLqBcC2KZ6jyYvM4mQ= @@ -163,20 +154,17 @@ golang.org/x/xerrors v0.0.0-20190717185122-a985d3407aa7/go.mod h1:I/5z698sn9Ka8T golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543 h1:E7g+9GITq07hpfrRu66IVDexMakfv52eLZ2CXBWiKr4= golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= gonum.org/v1/gonum v0.0.0-20180816165407-929014505bf4/go.mod h1:Y+Yx5eoAFn32cQvJDxZx5Dpnq+c3wtXuadVZAcxbbBo= +gonum.org/v1/gonum v0.8.2 h1:CCXrcPKiGGotvnN6jfUsKk4rRqm7q09/YbKb5xCEvtM= gonum.org/v1/gonum v0.8.2/go.mod h1:oe/vMfY3deqTw+1EZJhuvEW2iwGF1bW9wwu7XCu0+v0= -gonum.org/v1/gonum v0.9.3 h1:DnoIG+QAMaF5NvxnGe/oKsgKcAc6PcUyl8q0VetfQ8s= -gonum.org/v1/gonum v0.9.3/go.mod h1:TZumC3NeyVQskjXqmyWt4S3bINhy7B4eYwW69EbyX+0= gonum.org/v1/netlib v0.0.0-20190313105609-8cb42192e0e0 h1:OE9mWmgKkjJyEmDAAtGMPjXu+YNeGvK9VTSHY6+Qihc= gonum.org/v1/netlib v0.0.0-20190313105609-8cb42192e0e0/go.mod h1:wa6Ws7BG/ESfp6dHfk7C6KdzKA7wR7u/rKwOGE66zvw= gonum.org/v1/plot v0.0.0-20190515093506-e2840ee46a6b/go.mod h1:Wt8AAjI+ypCyYX3nZBvf6cAIx93T+c/OS2HFAYskSZc= +gonum.org/v1/plot v0.9.0 h1:3sEo36Uopv1/SA/dMFFaxXoL5XyikJ9Sf2Vll/k6+2E= gonum.org/v1/plot v0.9.0/go.mod h1:3Pcqqmp6RHvJI72kgb8fThyUnav364FOsdDo2aGW5lY= -gonum.org/v1/plot v0.10.0 h1:ymLukg4XJlQnYUJCp+coQq5M7BsUJFk6XQE4HPflwdw= -gonum.org/v1/plot v0.10.0/go.mod h1:JWIHJ7U20drSQb/aDpTetJzfC1KlAPldJLpkSy88dvQ= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= gopkg.in/check.v1 v1.0.0-20180628173108-788fd7840127/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= gopkg.in/natefinch/lumberjack.v2 v2.0.0 h1:1Lc07Kr7qY4U2YPouBjpCLxpiyxIVoxqXgkXLknAOE8= gopkg.in/natefinch/lumberjack.v2 v2.0.0/go.mod h1:l0ndWWf7gzL7RNwBG7wST/UCcT4T24xpD6X8LsfU/+k= gopkg.in/yaml.v2 v2.2.2 h1:ZCJp+EgiOT7lHqUV2J862kp8Qj64Jo6az82+3Td9dZw= gopkg.in/yaml.v2 v2.2.2/go.mod h1:hI93XBmqTisBFMUTm0b8Fm+jr3Dg1NNxqwp+5A1VGuI= -rsc.io/pdf v0.1.1 h1:k1MczvYDUvJBe93bYd7wrZLLUEcLZAuF824/I4e5Xr4= rsc.io/pdf v0.1.1/go.mod h1:n8OzWcQ6Sp37PL01nO98y4iUCRdTGarVfzxY20ICaU4= From b2d2a41fdcf6d3ceb88ba75d61cb989ce6cb4e9e Mon Sep 17 00:00:00 2001 From: David Sutton Date: Wed, 28 Dec 2022 00:30:30 +1030 Subject: [PATCH 7/9] Audiofiltering: Add amplification filter which uses filter.Upper as factor for amplification --- cmd/audiofiltering/main.go | 9 ++--- codec/pcm/filters.go | 78 +++++++++++++++++++++++--------------- 2 files changed, 51 insertions(+), 36 deletions(-) diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index 5ff73afe..1d1c129c 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -16,22 +16,21 @@ func main() { start := time.Now() // Read the audio data from the file. - input, _ := os.ReadFile("spike.pcm") + input, _ := os.ReadFile("sine.pcm") // Create Buffer in struct format defined in pcm.go format := pcm.BufferFormat{Rate: 44100, Channels: 1, SFormat: pcm.S16_LE} buf := pcm.Buffer{Format: format, Data: input} // Create a filter. - lp := pcm.Filter{BuffInfo: buf.Format, Type: pcm.LP, Upper: 10000, Taps: 50} - lp.Generate() + amp := pcm.Filter{BuffInfo: buf.Format, Type: pcm.AMPLIFIER, Upper: 3, Taps: 50} // Apply the lowpass filter to the buffer. - output := lp.Apply(buf) + output := amp.Apply(buf) // Save the transformed audio. f, _ := os.Create("output.pcm") - f.Write(output[50 : len(output)-50]) + f.Write(output) // Display execution time. fmt.Println("Finished execution. Total time:", time.Since(start)) diff --git a/codec/pcm/filters.go b/codec/pcm/filters.go index 8e61e53c..4e57d7c2 100644 --- a/codec/pcm/filters.go +++ b/codec/pcm/filters.go @@ -13,10 +13,11 @@ type FilterType int // Currently implemented filter types. const ( - LP FilterType = iota - HP - BP - BS + LOWPASS FilterType = iota + HIGHPASS + BANDPASS + BANDSTOP + AMPLIFIER ) // Filter contains the specifications of the filter, as well as the coefficients to the filter function itself. @@ -48,7 +49,7 @@ func (filter *Filter) Generate() { // Determine the type of filter to generate (based off filter.Type). switch filter.Type { - case LP: + case LOWPASS: // Create a lowpass filter with characteristics from struct. size := filter.Taps + 1 @@ -64,9 +65,9 @@ func (filter *Filter) Generate() { } filter.Coeffs[filter.Taps/2] = 2 * fd * winData[filter.Taps/2] - case HP: + case HIGHPASS: - // Create a Highpass filter with characteristics from struct. + // Create a HighIGHPASSass filter with characteristics from struct. size := filter.Taps + 1 filter.Coeffs = make([]float64, size, size) fd := (filter.Lower + 3000) / float64(filter.SampleRate) @@ -80,35 +81,35 @@ func (filter *Filter) Generate() { } filter.Coeffs[filter.Taps/2] = (1 - 2*fd) * winData[filter.Taps/2] - case BP: + case BANDPASS: - // Make Low and Highpass filters. - lp := Filter{Type: LP, SampleRate: filter.SampleRate, Upper: filter.Upper, Taps: filter.Taps} - hp := Filter{Type: HP, SampleRate: filter.SampleRate, Lower: filter.Lower, Taps: filter.Taps} - lp.Generate() - hp.Generate() + // Make Low and HighIGHPASSass filters. + lowpass := Filter{Type: LOWPASS, SampleRate: filter.SampleRate, Upper: filter.Upper, Taps: filter.Taps} + hIGHPASS := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Lower, Taps: filter.Taps} + lowpass.Generate() + hIGHPASS.Generate() - // Convolve lowpass filter with highpass filter to get bandpass filter. + // Convolve lowpass filter with highIGHPASSass filter to get bandpass filter. var wg sync.WaitGroup ch := make(chan []float64) wg.Add(1) - go Convolve(lp.Coeffs, hp.Coeffs, &wg, ch) + go Convolve(lowpass.Coeffs, hIGHPASS.Coeffs, &wg, ch) wg.Wait() filter.Coeffs = <-ch - case BS: + case BANDSTOP: - // Make Low and Highpass filters. - lp := Filter{Type: LP, SampleRate: filter.SampleRate, Upper: filter.Lower, Taps: filter.Taps} - hp := Filter{Type: HP, SampleRate: filter.SampleRate, Lower: filter.Upper, Taps: filter.Taps} - lp.Generate() - hp.Generate() + // Make Low and HighIGHPASSass filters. + lowpass := Filter{Type: LOWPASS, SampleRate: filter.SampleRate, Upper: filter.Lower, Taps: filter.Taps} + hIGHPASS := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Upper, Taps: filter.Taps} + lowpass.Generate() + hIGHPASS.Generate() - // Add lowpass filter to highpass filter to get bandstop filter. + // Add lowpass filter to highIGHPASSass filter to get bandstop filter. size := filter.Taps + 1 filter.Coeffs = make([]float64, size, size) - for i := range lp.Coeffs { - filter.Coeffs[i] = lp.Coeffs[i] + hp.Coeffs[i] + for i := range lowpass.Coeffs { + filter.Coeffs[i] = lowpass.Coeffs[i] + hIGHPASS.Coeffs[i] } } @@ -130,13 +131,28 @@ func (filter *Filter) Apply(b Buffer) []byte { inputAsFloat[i] /= (32767) } - // Convolve input with filter. - var wg sync.WaitGroup - ch := make(chan []float64, 1) - wg.Add(1) - go Convolve(inputAsFloat, filter.Coeffs, &wg, ch) - wg.Wait() - convolution := <-ch + var convolution []float64 + // Check if filter type is frequency filter or amplifier. + if filter.Type == AMPLIFIER { + convolution = make([]float64, len(inputAsFloat)) + for i := range inputAsFloat { + convolution[i] = inputAsFloat[i] * filter.Upper + // Stop audio artifacting by clipping outputs + if convolution[i] > 1 { + convolution[i] = 1 + } else if convolution[i] < -1 { + convolution[i] = -1 + } + } + } else { + // Convolve input with filter. + var wg sync.WaitGroup + ch := make(chan []float64, 1) + wg.Add(1) + go Convolve(inputAsFloat, filter.Coeffs, &wg, ch) + wg.Wait() + convolution = <-ch + } // Convert convolution output back to bytes. var output []byte From 52a56f3a52b0460deabf3ba77c179faa6f648a9f Mon Sep 17 00:00:00 2001 From: David Sutton Date: Wed, 28 Dec 2022 14:32:24 +1030 Subject: [PATCH 8/9] Audiofiltering: Increase efficiency of convolution algorithm. Approximately 150x faster. (takes ~4.5s to lowpass filter ~1min of audio) --- cmd/audiofiltering/main.go | 51 ++++++++++++++++++--- codec/pcm/filters.go | 90 +++++++++++++++++++++++++++++++------- 2 files changed, 120 insertions(+), 21 deletions(-) diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go index 1d1c129c..ff9a2d23 100644 --- a/cmd/audiofiltering/main.go +++ b/cmd/audiofiltering/main.go @@ -16,21 +16,60 @@ func main() { start := time.Now() // Read the audio data from the file. - input, _ := os.ReadFile("sine.pcm") + input, _ := os.ReadFile("whitenoise.pcm") // Create Buffer in struct format defined in pcm.go format := pcm.BufferFormat{Rate: 44100, Channels: 1, SFormat: pcm.S16_LE} buf := pcm.Buffer{Format: format, Data: input} // Create a filter. - amp := pcm.Filter{BuffInfo: buf.Format, Type: pcm.AMPLIFIER, Upper: 3, Taps: 50} + bs := pcm.Filter{BuffInfo: buf.Format, Type: pcm.BANDSTOP, Lower: 1000, Upper: 2000, Taps: 500} - // Apply the lowpass filter to the buffer. - output := amp.Apply(buf) + // Apply different filters to save and compare. + bs.Generate() + bs_1khz := bs.Apply(buf) + fmt.Println("Applied 1Khz Bandstop filter") + + bs.Lower, bs.Upper = 2000, 5000 + bs.Generate() + bs_2khz := bs.Apply(buf) + fmt.Println("Applied 2Khz Bandstop filter") + + bs.Lower, bs.Upper = 5000, 10000 + bs.Generate() + bs_5khz := bs.Apply(buf) + fmt.Println("Applied 5Khz Bandstop filter") + + bs.Lower, bs.Upper = 10000, 15000 + bs.Generate() + bs_10khz := bs.Apply(buf) + fmt.Println("Applied 10Khz Bandstop filter") + + bs.Lower, bs.Upper = 15000, 18000 + bs.Generate() + bs_15khz := bs.Apply(buf) + fmt.Println("Applied 15Khz Bandstop filter") // Save the transformed audio. - f, _ := os.Create("output.pcm") - f.Write(output) + f, _ := os.Create("bs_1khz.pcm") + f.Write(bs_1khz) + fmt.Println("Wrote audio to bs_1khz.pcm") + + f, _ = os.Create("bs_2khz.pcm") + f.Write(bs_2khz) + fmt.Println("Wrote audio to bs_2khz.pcm") + + f, _ = os.Create("bs_5khz.pcm") + f.Write(bs_5khz) + fmt.Println("Wrote audio to bs_5khz.pcm") + + f, _ = os.Create("bs_10khz.pcm") + f.Write(bs_10khz) + fmt.Println("Wrote audio to bs_10khz.pcm") + + f, _ = os.Create("bs_15khz.pcm") + f.Write(bs_15khz) + fmt.Println("Wrote audio to bs_15khz.pcm") // Display execution time. fmt.Println("Finished execution. Total time:", time.Since(start)) diff --git a/codec/pcm/filters.go b/codec/pcm/filters.go index 4e57d7c2..439fbe9c 100644 --- a/codec/pcm/filters.go +++ b/codec/pcm/filters.go @@ -2,9 +2,11 @@ package pcm import ( "encoding/binary" + "fmt" "math" "sync" + "github.com/mjibson/go-dsp/fft" "github.com/mjibson/go-dsp/window" ) @@ -54,7 +56,7 @@ func (filter *Filter) Generate() { // Create a lowpass filter with characteristics from struct. size := filter.Taps + 1 filter.Coeffs = make([]float64, size, size) - fd := (filter.Upper * 0.7) / float64(filter.SampleRate) + fd := (filter.Upper) / float64(filter.SampleRate) b := (2 * math.Pi) * fd winData := window.FlatTop(size) for n := 0; n < (filter.Taps / 2); n++ { @@ -70,7 +72,7 @@ func (filter *Filter) Generate() { // Create a HighIGHPASSass filter with characteristics from struct. size := filter.Taps + 1 filter.Coeffs = make([]float64, size, size) - fd := (filter.Lower + 3000) / float64(filter.SampleRate) + fd := (filter.Lower) / float64(filter.SampleRate) b := (2 * math.Pi) * fd winData := window.FlatTop(size) for n := 0; n < (filter.Taps / 2); n++ { @@ -85,31 +87,31 @@ func (filter *Filter) Generate() { // Make Low and HighIGHPASSass filters. lowpass := Filter{Type: LOWPASS, SampleRate: filter.SampleRate, Upper: filter.Upper, Taps: filter.Taps} - hIGHPASS := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Lower, Taps: filter.Taps} + highpass := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Lower, Taps: filter.Taps} lowpass.Generate() - hIGHPASS.Generate() + highpass.Generate() // Convolve lowpass filter with highIGHPASSass filter to get bandpass filter. var wg sync.WaitGroup - ch := make(chan []float64) + ch := make(chan []float64, 1) wg.Add(1) - go Convolve(lowpass.Coeffs, hIGHPASS.Coeffs, &wg, ch) + go FastConvolve(lowpass.Coeffs, highpass.Coeffs, &wg, ch) wg.Wait() filter.Coeffs = <-ch case BANDSTOP: - // Make Low and HighIGHPASSass filters. + // Make Low and Highpass filters. lowpass := Filter{Type: LOWPASS, SampleRate: filter.SampleRate, Upper: filter.Lower, Taps: filter.Taps} - hIGHPASS := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Upper, Taps: filter.Taps} + highpass := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Upper, Taps: filter.Taps} lowpass.Generate() - hIGHPASS.Generate() + highpass.Generate() - // Add lowpass filter to highIGHPASSass filter to get bandstop filter. + // Add lowpass filter to highpass filter to get bandstop filter. size := filter.Taps + 1 filter.Coeffs = make([]float64, size, size) for i := range lowpass.Coeffs { - filter.Coeffs[i] = lowpass.Coeffs[i] + hIGHPASS.Coeffs[i] + filter.Coeffs[i] = lowpass.Coeffs[i] + highpass.Coeffs[i] } } @@ -149,7 +151,7 @@ func (filter *Filter) Apply(b Buffer) []byte { var wg sync.WaitGroup ch := make(chan []float64, 1) wg.Add(1) - go Convolve(inputAsFloat, filter.Coeffs, &wg, ch) + go FastConvolve(inputAsFloat, filter.Coeffs, &wg, ch) wg.Wait() convolution = <-ch } @@ -158,9 +160,14 @@ func (filter *Filter) Apply(b Buffer) []byte { var output []byte buf := make([]byte, 2) for i := range convolution { + if convolution[i] >= 1 { + convolution[i] = 0.9999 + } else if convolution[i] <= -1 { + convolution[i] = -0.9999 + } convolution[i] = convolution[i] * 32767 if convolution[i] < 0 { - convolution[i] = convolution[i] + 32768*2 + convolution[i] = convolution[i] + 32767*2 } binary.LittleEndian.PutUint16(buf[:], uint16(convolution[i])) output = append(output, buf[0], buf[1]) @@ -170,6 +177,7 @@ func (filter *Filter) Apply(b Buffer) []byte { } +// Convolve takes in a signal and an FIR filter and computes the convolution. (runs in O(n^2) time) func Convolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { // Create a waitgroup to be used in goroutines called by Convolution @@ -178,9 +186,10 @@ func Convolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { // Compute the convolution convLen := len(x) + len(h) - 1 y := make([]float64, convLen) + var progress int for n := 0; n < convLen; n++ { convwg.Add(1) - go func(n int, y []float64, convwg *sync.WaitGroup) { + go func(n int, y []float64, convwg *sync.WaitGroup, progress *int) { var sum float64 = 0 for k := 0; k < len(x); k++ { if n-k >= 0 && n-k < len(h) { @@ -190,14 +199,65 @@ func Convolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { } } y[n] = sum + *progress++ convwg.Done() - }(n, y, &convwg) + }(n, y, &convwg, &progress) + fmt.Println(float64(progress) * 100 / float64(convLen)) } + convwg.Wait() ch <- y close(ch) wg.Done() } +// FastConvolve takes in a signal and an FIR filter and computes the convolution. (runs in O(nlog(n)) time) +func FastConvolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { + + // Calculate the length of the linear convolution + convLen := len(x) + len(h) - 1 + + // Pad signals to the next largest power of 2 larger than convLen + padLen := int(math.Pow(2, math.Ceil(math.Log2(float64(convLen))))) + zeros := make([]float64, padLen-len(x), padLen-len(h)) + x = append(x, zeros...) + zeros = make([]float64, padLen-len(h)) + h = append(h, zeros...) + + // Compute DFFTs + X, H := fft.FFTReal(x), fft.FFTReal(h) + + // Compute the multiplication of the two signals in the freq domain + var convWG sync.WaitGroup + Y := make([]complex128, padLen) + for i := range X { + convWG.Add(1) + go func(a, b complex128, y *complex128, convWG *sync.WaitGroup, i int) { + *y = a * b + convWG.Done() + }(X[i], H[i], &Y[i], &convWG, i) + } + + convWG.Wait() + + // Compute the IDFFT + iy := fft.IFFT(Y) + + // Convert to []float64 + y := make([]float64, padLen) + for i := range iy { + convWG.Add(1) + go func(a complex128, y *float64, convWG *sync.WaitGroup) { + *y = real(a) + convWG.Done() + }(iy[i], &y[i], &convWG) + } + convWG.Wait() + + // Trim to length of linear convolution + ch <- y[0:convLen] + wg.Done() +} + // LEGACY FUNCTION FOR TESTING // Max returns the absolute highest value in a given array. func Max(a []float64) float64 { From 70afcdb8164e29f9882616466cf43cdeb6178ed1 Mon Sep 17 00:00:00 2001 From: ausocean-david Date: Wed, 28 Dec 2022 20:37:39 +1030 Subject: [PATCH 9/9] Audiofiltering: Add amplifying capabilities, using Filter.Upper as the factor for amplification --- cmd/audiofiltering/main.go | 94 ------- codec/pcm/filters.go | 500 ++++++++++++++++++++----------------- codec/pcm/filters_test.go | 368 +++++++++++++++++++++++++++ codec/pcm/pcm.go | 6 +- go.mod | 1 + go.sum | 22 ++ 6 files changed, 671 insertions(+), 320 deletions(-) delete mode 100644 cmd/audiofiltering/main.go create mode 100644 codec/pcm/filters_test.go diff --git a/cmd/audiofiltering/main.go b/cmd/audiofiltering/main.go deleted file mode 100644 index ff9a2d23..00000000 --- a/cmd/audiofiltering/main.go +++ /dev/null @@ -1,94 +0,0 @@ -package main - -import ( - "fmt" - "math" - "os" - "time" - - "bitbucket.org/ausocean/av/codec/pcm" -) - -// main is a driver function for testing the filters defined in codec/pcm/filters.go -func main() { - - // Define start time for execution timing. - start := time.Now() - - // Read the audio data from the file. - input, _ := os.ReadFile("whitenoise.pcm") - - // Create Buffer in struct format defined in pcm.go - format := pcm.BufferFormat{Rate: 44100, Channels: 1, SFormat: pcm.S16_LE} - buf := pcm.Buffer{Format: format, Data: input} - - // Create a filter. - bs := pcm.Filter{BuffInfo: buf.Format, Type: pcm.BANDSTOP, Lower: 1000, Upper: 2000, Taps: 500} - - // Apply different filters to save and compare. - bs.Generate() - bs_1khz := bs.Apply(buf) - fmt.Println("Applied 1Khz Bandstop filter") - - bs.Lower, bs.Upper = 2000, 5000 - bs.Generate() - bs_2khz := bs.Apply(buf) - fmt.Println("Applied 2Khz Bandstop filter") - - bs.Lower, bs.Upper = 5000, 10000 - bs.Generate() - bs_5khz := bs.Apply(buf) - fmt.Println("Applied 5Khz Bandstop filter") - - bs.Lower, bs.Upper = 10000, 15000 - bs.Generate() - bs_10khz := bs.Apply(buf) - fmt.Println("Applied 10Khz Bandstop filter") - - bs.Lower, bs.Upper = 15000, 18000 - bs.Generate() - bs_15khz := bs.Apply(buf) - fmt.Println("Applied 15Khz Bandstop filter") - - // Save the transformed audio. - f, _ := os.Create("bs_1khz.pcm") - f.Write(bs_1khz) - fmt.Println("Wrote audio to bs_1khz.pcm") - - f, _ = os.Create("bs_2khz.pcm") - f.Write(bs_2khz) - fmt.Println("Wrote audio to bs_2khz.pcm") - - f, _ = os.Create("bs_5khz.pcm") - f.Write(bs_5khz) - fmt.Println("Wrote audio to bs_5khz.pcm") - - f, _ = os.Create("bs_10khz.pcm") - f.Write(bs_10khz) - fmt.Println("Wrote audio to bs_10khz.pcm") - - f, _ = os.Create("bs_15khz.pcm") - f.Write(bs_15khz) - fmt.Println("Wrote audio to bs_15khz.pcm") - - // Display execution time. - fmt.Println("Finished execution. Total time:", time.Since(start)) - -} - -// LEGACY FUNCTION FOR TESTING -// generate is used to generate a sine wave of the given frequency, calls on the global variables of Duration and SampleRate (bytes) -func generate(freq float64, SampleRate, duration int) []byte { - - // Create slices to store values. - signal := make([]byte, 2*duration*SampleRate) - t := make([]float64, 2*duration*SampleRate) - - // Generate the values to plot. - for n := range t { - t[n] = math.Pi * float64(n) / float64(SampleRate) - signal[n] = byte(math.Round(127 * math.Sin(freq*t[n]))) - } - - return signal -} diff --git a/codec/pcm/filters.go b/codec/pcm/filters.go index 439fbe9c..7c23a5aa 100644 --- a/codec/pcm/filters.go +++ b/codec/pcm/filters.go @@ -1,274 +1,328 @@ +/* +NAME + filters.go + +DESCRIPTION + filter.go contains functions for filtering PCM audio. + +AUTHOR + David Sutton + +LICENSE + filters_test.go is Copyright (C) 2023 the Australian Ocean Lab (AusOcean) + + It is free software: you can redistribute it and/or modify them + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + It is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License in gpl.txt. + If not, see [GNU licenses](http://www.gnu.org/licenses). +*/ + package pcm import ( + "bytes" "encoding/binary" + "errors" "fmt" "math" - "sync" "github.com/mjibson/go-dsp/fft" "github.com/mjibson/go-dsp/window" ) -// FilterType is the type of filter which can be generated. -type FilterType int - -// Currently implemented filter types. -const ( - LOWPASS FilterType = iota - HIGHPASS - BANDPASS - BANDSTOP - AMPLIFIER -) - -// Filter contains the specifications of the filter, as well as the coefficients to the filter function itself. -type Filter struct { - Coeffs []float64 - Type FilterType - SampleRate uint - Lower, Upper float64 - Taps int - BuffInfo BufferFormat -} - -// AudioFilter interface contains Generate and Apply. Generate is used to generate the coefficients of the filter based off -// the specifications within the Filter struct. +// AudioFilter is an interface which contains an Apply function. // Apply is used to apply the filter to the given buffer of PCM data (b.Data). type AudioFilter interface { - Generate() Apply(b Buffer) } -// Generate is used to generate the coefficients of the filter function used for convolution with a signal to -// perform specified filter. -func (filter *Filter) Generate() { - - // Update the sample rate from the buffer info (if supplied). - if filter.BuffInfo.Rate != 0 { - filter.SampleRate = filter.BuffInfo.Rate - } - - // Determine the type of filter to generate (based off filter.Type). - switch filter.Type { - case LOWPASS: - - // Create a lowpass filter with characteristics from struct. - size := filter.Taps + 1 - filter.Coeffs = make([]float64, size, size) - fd := (filter.Upper) / float64(filter.SampleRate) - b := (2 * math.Pi) * fd - winData := window.FlatTop(size) - for n := 0; n < (filter.Taps / 2); n++ { - c := float64(n) - float64(filter.Taps)/2 - y := math.Sin(c*b) / (math.Pi * c) - filter.Coeffs[n] = (y * winData[n]) - filter.Coeffs[size-1-n] = filter.Coeffs[n] - } - filter.Coeffs[filter.Taps/2] = 2 * fd * winData[filter.Taps/2] - - case HIGHPASS: - - // Create a HighIGHPASSass filter with characteristics from struct. - size := filter.Taps + 1 - filter.Coeffs = make([]float64, size, size) - fd := (filter.Lower) / float64(filter.SampleRate) - b := (2 * math.Pi) * fd - winData := window.FlatTop(size) - for n := 0; n < (filter.Taps / 2); n++ { - c := float64(n) - float64(filter.Taps)/2 - y := math.Sin(c*b) / (math.Pi * c) - filter.Coeffs[n] = -y * winData[n] - filter.Coeffs[size-1-n] = filter.Coeffs[n] - } - filter.Coeffs[filter.Taps/2] = (1 - 2*fd) * winData[filter.Taps/2] - - case BANDPASS: - - // Make Low and HighIGHPASSass filters. - lowpass := Filter{Type: LOWPASS, SampleRate: filter.SampleRate, Upper: filter.Upper, Taps: filter.Taps} - highpass := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Lower, Taps: filter.Taps} - lowpass.Generate() - highpass.Generate() - - // Convolve lowpass filter with highIGHPASSass filter to get bandpass filter. - var wg sync.WaitGroup - ch := make(chan []float64, 1) - wg.Add(1) - go FastConvolve(lowpass.Coeffs, highpass.Coeffs, &wg, ch) - wg.Wait() - filter.Coeffs = <-ch - - case BANDSTOP: - - // Make Low and Highpass filters. - lowpass := Filter{Type: LOWPASS, SampleRate: filter.SampleRate, Upper: filter.Lower, Taps: filter.Taps} - highpass := Filter{Type: HIGHPASS, SampleRate: filter.SampleRate, Lower: filter.Upper, Taps: filter.Taps} - lowpass.Generate() - highpass.Generate() - - // Add lowpass filter to highpass filter to get bandstop filter. - size := filter.Taps + 1 - filter.Coeffs = make([]float64, size, size) - for i := range lowpass.Coeffs { - filter.Coeffs[i] = lowpass.Coeffs[i] + highpass.Coeffs[i] - } - - } +// SelectiveFrequencyFilter is a struct which contains all the filter specifications required for a +// lowpass, highpass, bandpass, or bandstop filter. +type SelectiveFrequencyFilter struct { + coeffs []float64 + cutoff [2]float64 + sampleRate uint + taps int + buffInfo BufferFormat } -// Apply takes in a buffer of PCM audio, applies the filter and returns the filtered audio (in byte slice) -func (filter *Filter) Apply(b Buffer) []byte { +// NewLowPass populates a LowPass struct with the specified data. The function also +// generates a lowpass filter based off the given specifications, and returns a pointer. +func NewLowPass(fc float64, info BufferFormat, length int) (*SelectiveFrequencyFilter, error) { + return newLoHiFilter(fc, info, length, [2]float64{0, fc}) +} - // Convert input to floats. - inputAsFloat := make([]float64, len(b.Data)/2) - temp := make([]byte, 2) +// NewHighPass populates a HighPass struct with the specified data. The function also +// generates a highpass filter based off the given specifications, and returns a pointer. +func NewHighPass(fc float64, info BufferFormat, length int) (*SelectiveFrequencyFilter, error) { + return newLoHiFilter(fc, info, length, [2]float64{fc, 0}) +} + +// NewBandPass populates a BandPass struct with the specified data. The function also +// generates a bandpass filter based off the given specifications, and returns a pointer. +func NewBandPass(fc_lower, fc_upper float64, info BufferFormat, length int) (*SelectiveFrequencyFilter, error) { + newFilter, lp, hp, err := newBandFilter([2]float64{fc_lower, fc_upper}, info, length) + if err != nil { + return nil, fmt.Errorf("could not create new band filter: %w", err) + } + + // Convolve the filters to create a bandpass filter. + newFilter.coeffs, err = fastConvolve(hp.coeffs, lp.coeffs) + if err != nil { + return nil, fmt.Errorf("could not compute fast convolution: %w", err) + } + + // Return a pointer to the filter. + return newFilter, nil +} + +// NewBandStop populates a BandStop struct with the specified data. The function also +// generates a bandstop filter based off the given specifications, and returns a pointer. +func NewBandStop(fc_lower, fc_upper float64, info BufferFormat, length int) (*SelectiveFrequencyFilter, error) { + newFilter, lp, hp, err := newBandFilter([2]float64{fc_upper, fc_lower}, info, length) + if err != nil { + return nil, fmt.Errorf("could not create new band filter: %w", err) + } + size := newFilter.taps + 1 + newFilter.coeffs = make([]float64, size) + for i := range lp.coeffs { + newFilter.coeffs[i] = lp.coeffs[i] + hp.coeffs[i] + } + + // Return a pointer to the filter. + return newFilter, nil +} + +// Apply is the SelectiveFrequencyFilter implementation of the AudioFilter interface. This implementation +// takes the buffer data (b.Data), applies the highpass filter and returns a byte slice of filtered audio. +func (filter *SelectiveFrequencyFilter) Apply(b Buffer) ([]byte, error) { + // Apply the lowpass filter to the given audio buffer. + return convolveFromBytes(b.Data, filter.coeffs) +} + +// Amplifier is a struct which contains the factor of amplification to be used in the application +// of the filter. +type Amplifier struct { + factor float64 +} + +// NewAmplifier defines the factor of amplification for an amplifying filter. +func NewAmplifier(factor float64) Amplifier { + // Return populated Amplifier filter. + // Uses the absolute value of the factor to ensure compatibility. + return Amplifier{factor: math.Abs(factor)} +} + +// Apply implemented for an amplifier takes the buffer data (b.Data), applies +// the amplification and returns a byte slice of filtered audio. +func (amp *Amplifier) Apply(b Buffer) ([]byte, error) { + inputAsFloat, err := bytesToFloats(b.Data) + if err != nil { + return nil, fmt.Errorf("failed to convert to floats: %w", err) + } + + // Multiply every sample by the factor of amplification. + floatOutput := make([]float64, len(inputAsFloat)) for i := range inputAsFloat { - temp[0] = b.Data[2*i] - temp[1] = b.Data[2*i+1] - inputAsFloat[i] = float64(binary.LittleEndian.Uint16(temp)) - if inputAsFloat[i] > 32767 { - inputAsFloat[i] -= 32768 * 2 + floatOutput[i] = inputAsFloat[i] * amp.factor + // Stop audio artifacting by clipping outputs. + if floatOutput[i] > 1 { + floatOutput[i] = 1 + } else if floatOutput[i] < -1 { + floatOutput[i] = -1 } - inputAsFloat[i] /= (32767) + } + outBytes, err := floatsToBytes(floatOutput) + if err != nil { + return nil, fmt.Errorf("failed to convert to bytes: %w", err) + } + return outBytes, nil +} + +// newLoHiFilter is a function which checks for the validity of the input parameters, and calls the newLoHiFilter function +// to return a pointer to either a lowpass or a highpass filter. +func newLoHiFilter(fc float64, info BufferFormat, length int, cutoff [2]float64) (*SelectiveFrequencyFilter, error) { + // Ensure that all input values are valid. + if fc <= 0 || fc >= float64(info.Rate)/2 { + return nil, errors.New("cutoff frequency out of bounds") + } else if length <= 0 { + return nil, errors.New("cannot create filter with length <= 0") } - var convolution []float64 - // Check if filter type is frequency filter or amplifier. - if filter.Type == AMPLIFIER { - convolution = make([]float64, len(inputAsFloat)) - for i := range inputAsFloat { - convolution[i] = inputAsFloat[i] * filter.Upper - // Stop audio artifacting by clipping outputs - if convolution[i] > 1 { - convolution[i] = 1 - } else if convolution[i] < -1 { - convolution[i] = -1 - } - } + // Determine the type of filter to be generated. + var fd float64 + var factor1 float64 + var factor2 float64 + if cutoff[0] == 0 { // For a lowpass filter, cutoff[0] = 0, cutoff[1] = fc. + // The filter must be a lowpass filter. + fd = cutoff[1] / float64(info.Rate) + factor1 = 1 + factor2 = 2 * fd + } else if cutoff[1] == 0 { // For a highpass filter, cutoff[0] = fc, cutoff[1] = 0. + // The filter must be a highpass filter. + fd = cutoff[0] / float64(info.Rate) + factor1 = -1 + factor2 = 1 - 2*fd } else { - // Convolve input with filter. - var wg sync.WaitGroup - ch := make(chan []float64, 1) - wg.Add(1) - go FastConvolve(inputAsFloat, filter.Coeffs, &wg, ch) - wg.Wait() - convolution = <-ch + // Otherwise the filter must be a different type of filter. + return nil, errors.New("tried to use newLoHiFilter to generate bandpass or bandstop filter") } - // Convert convolution output back to bytes. - var output []byte - buf := make([]byte, 2) - for i := range convolution { - if convolution[i] >= 1 { - convolution[i] = 0.9999 - } else if convolution[i] <= -1 { - convolution[i] = -0.9999 - } - convolution[i] = convolution[i] * 32767 - if convolution[i] < 0 { - convolution[i] = convolution[i] + 32767*2 - } - binary.LittleEndian.PutUint16(buf[:], uint16(convolution[i])) - output = append(output, buf[0], buf[1]) + // Create a new filter struct to return, populated with all correct data. + var newFilter = SelectiveFrequencyFilter{cutoff: cutoff, sampleRate: info.Rate, taps: length, buffInfo: info} + + // Create a filter with characteristics from struct. + size := newFilter.taps + 1 + newFilter.coeffs = make([]float64, size) + b := 2 * math.Pi * fd + winData := window.FlatTop(size) + for n := 0; n < (newFilter.taps / 2); n++ { + c := float64(n) - float64(newFilter.taps)/2 + y := math.Sin(c*b) / (math.Pi * c) + newFilter.coeffs[n] = factor1 * y * winData[n] + newFilter.coeffs[size-1-n] = newFilter.coeffs[n] } + newFilter.coeffs[newFilter.taps/2] = factor2 * winData[newFilter.taps/2] - return output - + // Return a pointer to the filter. + return &newFilter, nil } -// Convolve takes in a signal and an FIR filter and computes the convolution. (runs in O(n^2) time) -func Convolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { - - // Create a waitgroup to be used in goroutines called by Convolution - var convwg sync.WaitGroup - - // Compute the convolution - convLen := len(x) + len(h) - 1 - y := make([]float64, convLen) - var progress int - for n := 0; n < convLen; n++ { - convwg.Add(1) - go func(n int, y []float64, convwg *sync.WaitGroup, progress *int) { - var sum float64 = 0 - for k := 0; k < len(x); k++ { - if n-k >= 0 && n-k < len(h) { - sum += x[k] * h[n-k] - } else if n-k < 0 { - break - } - } - y[n] = sum - *progress++ - convwg.Done() - }(n, y, &convwg, &progress) - fmt.Println(float64(progress) * 100 / float64(convLen)) +// newBandFilter creates a ensures the validity of the input parameters, and generates appropriate lowpass and highpass filters +// required for the creation of the specific band filter. +func newBandFilter(cutoff [2]float64, info BufferFormat, length int) (new, lp, hp *SelectiveFrequencyFilter, err error) { + // Ensure that all input values are valid. + if cutoff[0] <= 0 || cutoff[0] >= float64(info.Rate)/2 { + return nil, nil, nil, errors.New("cutoff frequencies out of bounds") + } else if cutoff[1] <= 0 || cutoff[1] >= float64(info.Rate)/2 { + return nil, nil, nil, errors.New("cutoff frequencies out of bounds") + } else if length <= 0 { + return nil, nil, nil, errors.New("cannot create filter with length <= 0") } - convwg.Wait() - ch <- y - close(ch) - wg.Done() + // Create a new filter struct to return, populated with all correct data. + // For a bandpass filter, cutoff[0] = fc_l, cutoff[1] = fc_u. + // For a bandstop filter, cutoff[0] = fc_u, cutoff[1] = fc_l. + var newFilter = SelectiveFrequencyFilter{cutoff: cutoff, sampleRate: info.Rate, taps: length, buffInfo: info} + + // Generate lowpass and highpass filters to create bandpass filter with. + hp, err = NewHighPass(newFilter.cutoff[0], newFilter.buffInfo, newFilter.taps) + if err != nil { + return nil, nil, nil, fmt.Errorf("could not create new highpass filter: %w", err) + } + lp, err = NewLowPass(newFilter.cutoff[1], newFilter.buffInfo, newFilter.taps) + if err != nil { + return nil, nil, nil, fmt.Errorf("could not create new lowpass filter: %w", err) + } + + // Return pointer to new filter. + return &newFilter, hp, lp, nil } -// FastConvolve takes in a signal and an FIR filter and computes the convolution. (runs in O(nlog(n)) time) -func FastConvolve(x, h []float64, wg *sync.WaitGroup, ch chan []float64) { +// convolveFromBytes takes in a byte slice and a float64 slice for a filter, converts to floats, +// convolves the two signals, and converts back to bytes and returns the convolution. +func convolveFromBytes(b []byte, filter []float64) ([]byte, error) { + bufAsFloats, err := bytesToFloats(b) + if err != nil { + return nil, fmt.Errorf("could not convert to floats: %w", err) + } - // Calculate the length of the linear convolution + // Convolve the floats with the filter. + convolution, err := fastConvolve(bufAsFloats, filter) + if err != nil { + return nil, fmt.Errorf("could not compute fast convolution: %w", err) + } + outBytes, err := floatsToBytes(convolution) + if err != nil { + return nil, fmt.Errorf("could not convert convolution to bytes: %w", err) + } + return outBytes, nil +} + +func bytesToFloats(b []byte) ([]float64, error) { + // Ensure the validity of the input. + if len(b) == 0 { + return nil, errors.New("no audio to convert to floats") + } else if len(b)%2 != 0 { + return nil, errors.New("uneven number of bytes (not whole number of samples)") + } + + // Convert bytes to floats. + inputAsFloat := make([]float64, len(b)/2) + inputAsInt := make([]int16, len(b)/2) + bReader := bytes.NewReader(b) + for i := range inputAsFloat { + binary.Read(bReader, binary.LittleEndian, &inputAsInt[i]) + inputAsFloat[i] = float64(inputAsInt[i]) / (math.MaxInt16 + 1) + } + return inputAsFloat, nil +} + +// floatsToBytes converts a slice of float64 PCM data into a slice of signed 16bit PCM data. +// The input float slice should contains values between -1 and 1. The function converts these values +// to a proportionate unsigned value between 0 and 65536. This 16bit integer is split into two bytes, +// then returned in Little Endian notation in a byte slice double the length of the input. +func floatsToBytes(f []float64) ([]byte, error) { + buf := new(bytes.Buffer) + bytes := make([]byte, len(f)*2) + for i := range f { + err := binary.Write(buf, binary.LittleEndian, int16(f[i]*math.MaxInt16)) + if err != nil { + return nil, fmt.Errorf("failed to write ints as bytes: %w", err) + } + } + n, err := buf.Read(bytes) + if err != nil { + return nil, fmt.Errorf("failed to read bytes from buffer: %w", err) + } else if n != len(bytes) { + return nil, fmt.Errorf("buffer and output length mismatch read %d bytes, expected %d: %w", n, len(bytes), err) + } + + return bytes, nil +} + +// fastConvolve takes in a signal and an FIR filter and computes the convolution (runs in O(nlog(n)) time). +func fastConvolve(x, h []float64) ([]float64, error) { + // Ensure valid data to convolve. + if len(x) == 0 || len(h) == 0 { + return nil, errors.New("convolution requires slice of length > 0") + } + + // Calculate the length of the linear convolution. convLen := len(x) + len(h) - 1 - // Pad signals to the next largest power of 2 larger than convLen + // Pad signals to the next largest power of 2 larger than convLen. padLen := int(math.Pow(2, math.Ceil(math.Log2(float64(convLen))))) zeros := make([]float64, padLen-len(x), padLen-len(h)) x = append(x, zeros...) zeros = make([]float64, padLen-len(h)) h = append(h, zeros...) - // Compute DFFTs - X, H := fft.FFTReal(x), fft.FFTReal(h) + // Compute DFFTs. + x_fft, h_fft := fft.FFTReal(x), fft.FFTReal(h) - // Compute the multiplication of the two signals in the freq domain - var convWG sync.WaitGroup - Y := make([]complex128, padLen) - for i := range X { - convWG.Add(1) - go func(a, b complex128, y *complex128, convWG *sync.WaitGroup, i int) { - *y = a * b - convWG.Done() - }(X[i], H[i], &Y[i], &convWG, i) + // Compute the multiplication of the two signals in the freq domain. + y_fft := make([]complex128, padLen) + for i := range x_fft { + y_fft[i] = x_fft[i] * h_fft[i] } - convWG.Wait() + // Compute the IDFFT. + iy := fft.IFFT(y_fft) - // Compute the IDFFT - iy := fft.IFFT(Y) - - // Convert to []float64 + // Convert to []float64. y := make([]float64, padLen) for i := range iy { - convWG.Add(1) - go func(a complex128, y *float64, convWG *sync.WaitGroup) { - *y = real(a) - convWG.Done() - }(iy[i], &y[i], &convWG) - } - convWG.Wait() - - // Trim to length of linear convolution - ch <- y[0:convLen] - wg.Done() -} - -// LEGACY FUNCTION FOR TESTING -// Max returns the absolute highest value in a given array. -func Max(a []float64) float64 { - - var runMax float64 = -1 - for i := range a { - if math.Abs(a[i]) > runMax { - runMax = math.Abs(a[i]) - } + y[i] = real(iy[i]) } - return runMax - + // Trim to length of linear convolution and return. + return y[0:convLen], nil } diff --git a/codec/pcm/filters_test.go b/codec/pcm/filters_test.go new file mode 100644 index 00000000..83cdac29 --- /dev/null +++ b/codec/pcm/filters_test.go @@ -0,0 +1,368 @@ +/* +NAME + filters_test.go + +DESCRIPTION + filter_test.go contains functions for testing functions in filters.go. + +AUTHOR + David Sutton + +LICENSE + filters_test.go is Copyright (C) 2023 the Australian Ocean Lab (AusOcean) + + It is free software: you can redistribute it and/or modify them + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + It is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License in gpl.txt. + If not, see [GNU licenses](http://www.gnu.org/licenses). +*/ + +package pcm + +import ( + "math" + "math/cmplx" + "os" + "testing" + + "github.com/mjibson/go-dsp/fft" +) + +// Set constant values for testing. +const ( + sampleRate = 44100 + filterLength = 500 + freqTest = 1000 +) + +// TestLowPass is used to test the lowpass constructor and application. Testing is done by ensuring frequency response as well as +// comparing against an expected audio file. +func TestLowPass(t *testing.T) { + // Generate an audio buffer to run test on. + genAudio, err := generate() + if err != nil { + t.Fatal(err) + } + var buf = Buffer{Data: genAudio, Format: BufferFormat{SFormat: S16_LE, Rate: sampleRate, Channels: 1}} + + // Create a lowpass filter to test. + const fc = 4500.0 + lp, err := NewLowPass(fc, buf.Format, filterLength) + if err != nil { + t.Fatal(err) + } + + // Filter the audio. + filteredAudio, err := lp.Apply(buf) + if err != nil { + t.Fatal(err) + } + + // Take the FFT of the signal. + filteredFloats, err := bytesToFloats(filteredAudio) + if err != nil { + t.Fatal(err) + } + filteredFFT := fft.FFTReal(filteredFloats) + + // Check if the lowpass filter worked (any high values in filteredFFT above cutoff freq result in fail). + for i := int(fc); i < sampleRate/2; i++ { + mag := math.Pow(cmplx.Abs(filteredFFT[i]), 2) + if mag > freqTest { + t.Error("Lowpass filter failed to meet spec.") + break + } + } + + // Read audio from the test location. + const fileName = "../../../test/test-data/av/input/lp_4500.pcm" + expectedAudio, err := os.ReadFile(fileName) + if err != nil { + t.Fatalf("File for comparison not read.\n\t%s", err) + } + + // Compare the filtered audio against the expected signal. + compare(filteredAudio, expectedAudio, t) +} + +// TestHighPass is used to test the highpass constructor and application. Testing is done by ensuring frequency response as well as +// comparing against an expected audio file. +func TestHighPass(t *testing.T) { + // Generate an audio buffer to run test on. + genAudio, err := generate() + if err != nil { + t.Fatal(err) + } + var buf = Buffer{Data: genAudio, Format: BufferFormat{SFormat: S16_LE, Rate: sampleRate, Channels: 1}} + + // Create a highpass filter to test. + const fc = 4500.0 + hp, err := NewHighPass(fc, buf.Format, filterLength) + if err != nil { + t.Fatal(err) + } + + // Filter the audio. + filteredAudio, err := hp.Apply(buf) + if err != nil { + t.Fatal(err) + } + + // Take the FFT of signal. + filteredFloats, err := bytesToFloats(filteredAudio) + if err != nil { + t.Fatal(err) + } + filteredFFT := fft.FFTReal(filteredFloats) + + // Check if the highpass filter worked (any high values in filteredFFT below cutoff freq result in fail). + for i := 0; i < int(fc); i++ { + mag := math.Pow(cmplx.Abs(filteredFFT[i]), 2) + if mag > freqTest { + t.Error("Highpass Filter doesn't meet Spec", i) + } + } + + // Read audio from the test location. + const fileName = "../../../test/test-data/av/input/hp_4500.pcm" + expectedAudio, err := os.ReadFile(fileName) + if err != nil { + t.Fatalf("File for comparison not read.\n\t%s", err) + } + + // Compare against expected results. + compare(expectedAudio, filteredAudio, t) +} + +// TestBandPass is used to test the bandpass constructor and application. Testing is done by ensuring frequency response as well as +// comparing against an expected audio file. +func TestBandPass(t *testing.T) { + // Generate an audio buffer to run test on. + genAudio, err := generate() + if err != nil { + t.Fatal(err) + } + var buf = Buffer{Data: genAudio, Format: BufferFormat{SFormat: S16_LE, Rate: sampleRate, Channels: 1}} + + // Create a bandpass filter to test. + const ( + fc_l = 4500.0 + fc_u = 9500.0 + ) + hp, err := NewBandPass(fc_l, fc_u, buf.Format, filterLength) + if err != nil { + t.Fatal(err) + } + + // Filter audio with filter. + filteredAudio, err := hp.Apply(buf) + if err != nil { + t.Fatal(err) + } + + // Take FFT of signal. + filteredFloats, err := bytesToFloats(filteredAudio) + if err != nil { + t.Fatal(err) + } + filteredFFT := fft.FFTReal(filteredFloats) + + // Check if the bandpass filter worked (any high values in filteredFFT above cutoff or below cutoff freq result in fail). + for i := 0; i < int(fc_l); i++ { + mag := math.Pow(cmplx.Abs(filteredFFT[i]), 2) + if mag > freqTest { + t.Error("Bandpass Filter doesn't meet Spec", i) + } + } + + for i := int(fc_u); i < sampleRate/2; i++ { + mag := math.Pow(cmplx.Abs(filteredFFT[i]), 2) + if mag > freqTest { + t.Error("Bandpass Filter doesn't meet Spec", i) + } + } + + // Read audio from test location. + const fileName = "../../../test/test-data/av/input/bp_4500-9500.pcm" + expectedAudio, err := os.ReadFile(fileName) + if err != nil { + t.Fatalf("File for comparison not read.\n\t%s", err) + } + + // Compare against the expected audio. + compare(expectedAudio, filteredAudio, t) +} + +// TestBandPass is used to test the bandpass constructor and application. Testing is done by ensuring frequency response as well as +// comparing against an expected audio file. +func TestBandStop(t *testing.T) { + // Generate an audio buffer to run test on. + genAudio, err := generate() + if err != nil { + t.Fatal(err) + } + var buf = Buffer{Data: genAudio, Format: BufferFormat{SFormat: S16_LE, Rate: sampleRate, Channels: 1}} + + // Create a bandpass filter to test. + const ( + fc_l = 4500.0 + fc_u = 9500.0 + ) + bs, err := NewBandStop(fc_l, fc_u, buf.Format, filterLength) + if err != nil { + t.Fatal(err) + } + + // Filter audio with filter. + filteredAudio, err := bs.Apply(buf) + if err != nil { + t.Fatal(err) + } + + // Take FFT of signal. + filteredFloats, err := bytesToFloats(filteredAudio) + if err != nil { + t.Fatal(err) + } + filteredFFT := fft.FFTReal(filteredFloats) + + // Check if the bandpass filter worked (any high values in filteredFFT above cutoff or below cutoff freq result in fail). + for i := int(fc_l); i < int(fc_u); i++ { + mag := math.Pow(cmplx.Abs(filteredFFT[i]), 2) + if mag > freqTest { + t.Error("BandStop Filter doesn't meet Spec", i) + } + } + + // Read audio from test location. + const fileName = "../../../test/test-data/av/input/bs_4500-9500.pcm" + expectedAudio, err := os.ReadFile(fileName) + if err != nil { + t.Fatalf("File for comparison not read.\n\t%s", err) + } + + // Compare against the expected audio. + compare(expectedAudio, filteredAudio, t) +} + +// TestAmplifier is used to test the amplifier constructor and application. Testing is done by checking the maximum value before and +// after application, as well as comparing against an expected audio file. +func TestAmplifier(t *testing.T) { + // Load a simple sine wave with amplitude of 0.1 and load into buffer. + const audioFileName = "../../../test/test-data/av/input/sine.pcm" + lowSine, err := os.ReadFile(audioFileName) + if err != nil { + t.Errorf("File for filtering not read.\n\t%s", err) + t.FailNow() + } + var buf = Buffer{Data: lowSine, Format: BufferFormat{SFormat: S16_LE, Rate: sampleRate, Channels: 1}} + + // Create an amplifier filter. + const factor = 5.0 + amp := NewAmplifier(factor) + + // Apply the amplifier to the audio. + filteredAudio, err := amp.Apply(buf) + if err != nil { + t.Fatal(err) + } + + // Find the maximum sample before and after amplification. + dataFloats, err := bytesToFloats(buf.Data) + if err != nil { + t.Fatal(err) + } + preMax := max(dataFloats) + filteredFloats, err := bytesToFloats(filteredAudio) + if err != nil { + t.Fatal(err) + } + postMax := max(filteredFloats) + + // Compare the values. + if preMax*factor > 1 && postMax > 0.99 { + } else if postMax/preMax > 1.01*factor || postMax/preMax < 0.99*factor { + t.Error("Amplifier failed to meet spec, expected:", factor, " got:", postMax/preMax) + } + + // Load expected audio file. + const compFileName = "../../../test/test-data/av/input/amp_5.pcm" + expectedAudio, err := os.ReadFile(compFileName) + if err != nil { + t.Fatalf("File for comparison not read.\n\t%s", err) + } + + // Compare against the expected audio file. + compare(filteredAudio, expectedAudio, t) +} + +// generate returns a byte slice in the same format that would be read from a PCM file. +// The function generates a sound with a range of frequencies for testing against, +// with a length of 1 second. +func generate() ([]byte, error) { + // Create an slice to generate values across. + t := make([]float64, sampleRate) + s := make([]float64, sampleRate) + // Define spacing of generated frequencies. + const ( + deltaFreq = 1000 + maxFreq = 21000 + amplitude = float64(deltaFreq) / float64((maxFreq - deltaFreq)) + ) + for n := 0; n < sampleRate; n++ { + t[n] = float64(n) / float64(sampleRate) + // Generate sinewaves of different frequencies. + s[n] = 0 + for f := deltaFreq; f < maxFreq; f += deltaFreq { + s[n] += amplitude * math.Sin(float64(f)*2*math.Pi*t[n]) + } + } + // Return the spectrum as bytes (PCM). + bytesOut, err := floatsToBytes(s) + if err != nil { + return nil, err + } + return bytesOut, nil +} + +// compare takes in two audio files (S16_LE), compares them, and returns an error if the +// signals are more than 10% different at any individual sample. +func compare(a, b []byte, t *testing.T) { + // Convert to floats to compare. + aFloats, err := bytesToFloats(a) + if err != nil { + t.Fatal(err) + } + bFloats, err := bytesToFloats(b) + if err != nil { + t.Fatal(err) + } + + // Compare against filtered audio. + for i := range aFloats { + diff := (bFloats[i] - aFloats[i]) + if math.Abs(diff) > 0.1 { + t.Error("Filtered audio is too different to database") + return + } + } +} + +// max takes a float slice and returns the absolute largest value in the slice. +func max(a []float64) float64 { + var runMax float64 = -1 + for i := range a { + if math.Abs(a[i]) > runMax { + runMax = math.Abs(a[i]) + } + } + return runMax +} diff --git a/codec/pcm/pcm.go b/codec/pcm/pcm.go index dda88e93..eb93bce7 100644 --- a/codec/pcm/pcm.go +++ b/codec/pcm/pcm.go @@ -73,9 +73,9 @@ func DataSize(rate, channels, bitDepth uint, period float64) int { // Resample takes Buffer c and resamples the pcm audio data to 'rate' Hz and returns a Buffer with the resampled data. // Notes: -// - Currently only downsampling is implemented and c's rate must be divisible by 'rate' or an error will occur. -// - If the number of bytes in c.Data is not divisible by the decimation factor (ratioFrom), the remaining bytes will -// not be included in the result. Eg. input of length 480002 downsampling 6:1 will result in output length 80000. +// - Currently only downsampling is implemented and c's rate must be divisible by 'rate' or an error will occur. +// - If the number of bytes in c.Data is not divisible by the decimation factor (ratioFrom), the remaining bytes will +// not be included in the result. Eg. input of length 480002 downsampling 6:1 will result in output length 80000. func Resample(c Buffer, rate uint) (Buffer, error) { if c.Format.Rate == rate { return c, nil diff --git a/go.mod b/go.mod index 7a7b1257..e41bb19c 100644 --- a/go.mod +++ b/go.mod @@ -16,6 +16,7 @@ require ( github.com/pkg/errors v0.9.1 github.com/yobert/alsa v0.0.0-20180630182551-d38d89fa843e gocv.io/x/gocv v0.29.0 + golang.org/x/tools v0.5.0 // indirect gonum.org/v1/gonum v0.8.2 gonum.org/v1/plot v0.9.0 gopkg.in/natefinch/lumberjack.v2 v2.0.0 diff --git a/go.sum b/go.sum index cc2a47a7..e1454d15 100644 --- a/go.sum +++ b/go.sum @@ -96,6 +96,7 @@ github.com/tarm/serial v0.0.0-20180830185346-98f6abe2eb07/go.mod h1:kDXzergiv9cb github.com/yobert/alsa v0.0.0-20180630182551-d38d89fa843e h1:3NIzz7weXhh3NToPgbtlQtKiVgerEaG4/nY2skGoGG0= github.com/yobert/alsa v0.0.0-20180630182551-d38d89fa843e/go.mod h1:CaowXBWOiSGWEpBBV8LoVnQTVPV4ycyviC9IBLj8dRw= github.com/yryz/ds18b20 v0.0.0-20180211073435-3cf383a40624/go.mod h1:MqFju5qeLDFh+S9PqxYT7TEla8xeW7bgGr/69q3oki0= +github.com/yuin/goldmark v1.4.13/go.mod h1:6yULJ656Px+3vBD8DxQVa3kxgyrAnzto9xy5taEt/CY= go.uber.org/atomic v1.3.2 h1:2Oa65PReHzfn29GpvgsYwloV9AVFHPDk8tYxt2c2tr4= go.uber.org/atomic v1.3.2/go.mod h1:gD2HeocX3+yG+ygLZcrzQJaqmWj9AIm7n08wl/qW/PE= go.uber.org/multierr v1.1.0 h1:HoEmRHQPVSqub6w2z2d2EOVs2fjyFRGyofhKuyDq0QI= @@ -109,6 +110,7 @@ golang.org/x/crypto v0.0.0-20190308221718-c2843e01d9a2/go.mod h1:djNgcEr1/C05ACk golang.org/x/crypto v0.0.0-20190510104115-cbcb75029529/go.mod h1:yigFU9vqHzYiE8UmvKecakEJjdnWj3jj499lnFckfCI= golang.org/x/crypto v0.0.0-20200622213623-75b288015ac9/go.mod h1:LzIPMQfyMNhhGPhUkYOs5KpL4U8rLKemX1yGLhDgUto= golang.org/x/crypto v0.0.0-20210711020723-a769d52b0f97/go.mod h1:GvvjBRRGRdwPK5ydBHafDWAxML/pGHZbMvKqRZ5+Abc= +golang.org/x/crypto v0.0.0-20210921155107-089bfa567519/go.mod h1:GvvjBRRGRdwPK5ydBHafDWAxML/pGHZbMvKqRZ5+Abc= golang.org/x/exp v0.0.0-20180321215751-8460e604b9de/go.mod h1:CJ0aWSM057203Lf6IL+f9T1iT9GByDxfZKAQTCR3kQA= golang.org/x/exp v0.0.0-20180807140117-3d87b88a115f/go.mod h1:CJ0aWSM057203Lf6IL+f9T1iT9GByDxfZKAQTCR3kQA= golang.org/x/exp v0.0.0-20190125153040-c74c464bbbf2/go.mod h1:CJ0aWSM057203Lf6IL+f9T1iT9GByDxfZKAQTCR3kQA= @@ -127,11 +129,18 @@ golang.org/x/image v0.0.0-20210216034530-4410531fe030 h1:lP9pYkih3DUSC641giIXa2X golang.org/x/image v0.0.0-20210216034530-4410531fe030/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0= golang.org/x/mobile v0.0.0-20190719004257-d2bd2a29d028/go.mod h1:E/iHnbuqvinMTCcRqshq8CkpyQDoeVncDDYHnLhea+o= golang.org/x/mod v0.1.0/go.mod h1:0QHyrYULN0/3qlju5TqG8bIK38QM8yzMo5ekMj3DlcY= +golang.org/x/mod v0.6.0-dev.0.20220419223038-86c51ed26bb4/go.mod h1:jJ57K6gSWd91VN4djpZkiMVwK6gcyfeH4XE8wZrZaV4= +golang.org/x/mod v0.7.0 h1:LapD9S96VoQRhi/GrNTqeBJFrUjs5UHCAtTlgwA5oZA= +golang.org/x/mod v0.7.0/go.mod h1:iBbtSCu2XBx23ZKBPSOrRkjjQPZFPuis4dIYUhu/chs= golang.org/x/net v0.0.0-20190404232315-eb5bcb51f2a3/go.mod h1:t9HGtf8HONx5eT2rtn7q6eTqICYqUVnKs3thJo3Qplg= golang.org/x/net v0.0.0-20190620200207-3b0461eec859/go.mod h1:z5CRVTTTmAJ677TzLLGU+0bjPO0LkuOLi4/5GtJWs/s= golang.org/x/net v0.0.0-20200904194848-62affa334b73/go.mod h1:/O7V0waA8r7cgGh81Ro3o1hOxt32SMVPicZroKQ2sZA= golang.org/x/net v0.0.0-20210226172049-e18ecbb05110/go.mod h1:m0MpNAwzfU5UDzcl9v0D8zg8gWTRqZa9RBIspLL5mdg= +golang.org/x/net v0.0.0-20220722155237-a158d28d115b/go.mod h1:XRhObCWvk6IyKnWLug+ECip1KBveYUHfp+8e9klMJ9c= +golang.org/x/net v0.5.0/go.mod h1:DivGGAXEgPSlEBzxGzZI+ZLohi+xUj054jfeKui00ws= golang.org/x/sync v0.0.0-20190423024810-112230192c58/go.mod h1:RxMgew5VJxzue5/jJTE5uejpjVlOe/izrB70Jof72aM= +golang.org/x/sync v0.0.0-20220722155255-886fb9371eb4/go.mod h1:RxMgew5VJxzue5/jJTE5uejpjVlOe/izrB70Jof72aM= +golang.org/x/sync v0.1.0/go.mod h1:RxMgew5VJxzue5/jJTE5uejpjVlOe/izrB70Jof72aM= golang.org/x/sys v0.0.0-20190215142949-d0b11bdaac8a/go.mod h1:STP8DvDyc/dI5b8T5hshtkjS+E42TnysNCUPdjciGhY= golang.org/x/sys v0.0.0-20190312061237-fead79001313/go.mod h1:h1NjWce9XRLGQEsW7wpKNCjG9DtNlClVuFLEZdDNbEs= golang.org/x/sys v0.0.0-20190412213103-97732733099d/go.mod h1:h1NjWce9XRLGQEsW7wpKNCjG9DtNlClVuFLEZdDNbEs= @@ -141,15 +150,28 @@ golang.org/x/sys v0.0.0-20201119102817-f84b799fce68/go.mod h1:h1NjWce9XRLGQEsW7w golang.org/x/sys v0.0.0-20210304124612-50617c2ba197/go.mod h1:h1NjWce9XRLGQEsW7wpKNCjG9DtNlClVuFLEZdDNbEs= golang.org/x/sys v0.0.0-20210615035016-665e8c7367d1/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg= golang.org/x/sys v0.0.0-20210909193231-528a39cd75f3/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg= +golang.org/x/sys v0.0.0-20220520151302-bc2c85ada10a/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg= +golang.org/x/sys v0.0.0-20220722155257-8c9f86f7a55f/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg= +golang.org/x/sys v0.4.0 h1:Zr2JFtRQNX3BCZ8YtxRE9hNJYC8J6I1MVbMg6owUp18= +golang.org/x/sys v0.4.0/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg= golang.org/x/term v0.0.0-20201126162022-7de9c90e9dd1/go.mod h1:bj7SfCRtBDWHUb9snDiAeCFNEtKQo2Wmx5Cou7ajbmo= +golang.org/x/term v0.0.0-20210927222741-03fcf44c2211/go.mod h1:jbD1KX2456YbFQfuXm/mYQcufACuNUgVhRMnK/tPxf8= +golang.org/x/term v0.4.0/go.mod h1:9P2UbLfCdcvo3p/nzKvsmas4TnlujnuoV9hGgYzW1lQ= golang.org/x/text v0.3.0/go.mod h1:NqM8EUOU14njkJ3fqMW+pc6Ldnwhi/IjpwHt7yyuwOQ= golang.org/x/text v0.3.3/go.mod h1:5Zoc/QRtKVWzQhOtBMvqHzDpF6irO9z98xDceosuGiQ= golang.org/x/text v0.3.5 h1:i6eZZ+zk0SOf0xgBpEpPD18qWcJda6q1sxt3S0kzyUQ= golang.org/x/text v0.3.5/go.mod h1:5Zoc/QRtKVWzQhOtBMvqHzDpF6irO9z98xDceosuGiQ= +golang.org/x/text v0.3.7/go.mod h1:u+2+/6zg+i71rQMx5EYifcz6MCKuco9NR6JIITiCfzQ= +golang.org/x/text v0.6.0 h1:3XmdazWV+ubf7QgHSTWeykHOci5oeekaGJBLkrkaw4k= +golang.org/x/text v0.6.0/go.mod h1:mrYo+phRRbMaCq/xk9113O4dZlRixOauAjOtrjsXDZ8= golang.org/x/tools v0.0.0-20180525024113-a5b4c53f6e8b/go.mod h1:n7NCudcB/nEzxVGmLbDWY5pfWTLqBcC2KZ6jyYvM4mQ= golang.org/x/tools v0.0.0-20180917221912-90fa682c2a6e/go.mod h1:n7NCudcB/nEzxVGmLbDWY5pfWTLqBcC2KZ6jyYvM4mQ= golang.org/x/tools v0.0.0-20190206041539-40960b6deb8e/go.mod h1:n7NCudcB/nEzxVGmLbDWY5pfWTLqBcC2KZ6jyYvM4mQ= golang.org/x/tools v0.0.0-20190927191325-030b2cf1153e/go.mod h1:b+2E5dAYhXwXZwtnZ6UAqBI28+e2cm9otk0dWdXHAEo= +golang.org/x/tools v0.0.0-20191119224855-298f0cb1881e/go.mod h1:b+2E5dAYhXwXZwtnZ6UAqBI28+e2cm9otk0dWdXHAEo= +golang.org/x/tools v0.1.12/go.mod h1:hNGJHUnrk76NpqgfD5Aqm5Crs+Hm0VOH/i9J2+nxYbc= +golang.org/x/tools v0.5.0 h1:+bSpV5HIeWkuvgaMfI3UmKRThoTA5ODJTUd8T17NO+4= +golang.org/x/tools v0.5.0/go.mod h1:N+Kgy78s5I24c24dU8OfWNEotWjutIs8SnJvn5IDq+k= golang.org/x/xerrors v0.0.0-20190717185122-a985d3407aa7/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543 h1:E7g+9GITq07hpfrRu66IVDexMakfv52eLZ2CXBWiKr4= golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0=