From 70afcdb8164e29f9882616466cf43cdeb6178ed1 Mon Sep 17 00:00:00 2001 From: ausocean-david Date: Wed, 28 Dec 2022 20:37:39 +1030 Subject: [PATCH] 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=