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 {