From 80ab5e476877e5f673e25906c97727913cfeb2e9 Mon Sep 17 00:00:00 2001
From: ausocean-david <davidsutton@ausocean.org>
Date: Wed, 14 Dec 2022 01:53:10 +1030
Subject: [PATCH] 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<n; i++ {
+	// 	audio[i] = generate(float64((i)*(22000/n)))
+	// }
+	audio[0] = generate(2000)
+	audio[1] = generate(10000)
+
+
+	// combine audio
+	combinedAudio := make([]float64, length)
+	for i := range audio[0] {
+		combinedAudio[i] = 0
+		for j:=0; j<n; j++ {
+			combinedAudio[i] += audio[j][i]
+		}
+	}
+	fmt.Println("audio generation:", time.Since(start))
+	start = time.Now()
+
+	filterLen := 500
+	filter := LowPass(filterLen, 20000)
+	fmt.Println("audio filter generation:", time.Since(start))
+	start = time.Now()
+
+	// convolve sinc with audio
+	filteredAudio := Convolve(combinedAudio, filter)
+	fmt.Println("convolution:", time.Since(start))
+	start = time.Now()
+	SaveAudioData(filteredAudio, "lowpass")
+	SaveAudioData(combinedAudio, "unfiltered")
+	fmt.Println("audio saving:", time.Since(start))
+	start = time.Now()
+
+	// SaveAudioData(filter, "test")
+
+	fmt.Println(time.Since(start))
+
+}
+
+func generate(Frequency float64) []float64 {
+	
+	// deteremine number of samples based off duration and sample rate
+	nsamps := Duration * SampleRate
+	
+	// generate x-values
+	var angle float64 = tau / float64(nsamps)
+
+	// create sample array
+	samp := make([]float64, int(nsamps))
+
+	// generate samples and write to file
+	for i := 0; i < int(nsamps); i++ {
+		samp[i] = math.Sin(angle * Frequency * float64(i))
+	}
+
+	return samp
+
+}
+
+func TopFive (a []float64) (topVal []float64, topI []int) {
+
+	length := 5
+	runMax := make([]float64, length, length)
+	indices := make([]int, length, length)
+	
+	for i:=range a {
+		switch {
+			case a[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<convLen; n++ {
+		go SubConvolve(n, x, h, y)
+	}
+
+	return y
+
+}
+
+func SubConvolve (n int, x, h, y []float64) {
+	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]
+		}
+	}
+	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<length; i++ {
+		binary.LittleEndian.PutUint32(buf[:], math.Float32bits(float32(signal[i])))
+		_, _ = f.Write(buf[:])
+	}
+
+}
+
+func Sinc (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
+
+}
\ No newline at end of file