Audio-Filter - please try, comment, optimise, notify bugs

[jamirokwai]

Hi there,

as I am working on an audio-tool for retro-needs, I turned to ChatGPT for getting help with audio-filters. Although, I am an engineer, my knowledge in audio-algorithms is... neglectable...

But, ChatGPT had a full day to set things right. And it somehow delivered. After about 20 iterations, the following code emerged.

I'd appreciate, if some of you would have a look, try things out and give me some feedback about improving the code, enhancing the experience and creating better code for all of us. I'd like to donate this part to the community - but it should be flawless!

Thanks!

Features so far:
- cross platform
- works with 8, 16, 24, 32 bits
- mono and stereo
- should also work within 8 to 192 khz
- uses a hamming-filter, you can set the window-size (the lesser, the more noise, but also less computing efforts)
- functions for bandpass, lowpass, highpass, notch, resonance (how does this work?), 3-band-EQ (highly questionable implementation?)
- create a sound, e.g. with "GenerateWhiteNoise", apply a filter, save as .wav
- or, load a .wav and put the data into a buffer.

- optimisations in mind:
- use a lookup table for sinus
- add more filters (like a BassBoost (tm), a TrebleEnlighter (tm), frequencies sweeps, more EQs, what you need?
- all routines may be optimised enough to use them in a realtime environment, but that needs to be tested

Code: Select all

Global DefaultSampleRate.i = 44100
Global DefaultBitDepth.i = 16
Global DefaultChannels.i = 2

#TwoPI = 2 * #PI

;- Filter

; globale variablen
Global Dim coefficients.i(0) ; filterkoeffizienten
Global Dim hamming.i(0)      ; hamming-fenster
Global windowsize

; hamming-fenster vorberechnen
Procedure PrecomputeHamming(newwindowsize.i)
  Protected i
  
  windowsize = newwindowsize
  ReDim hamming.i(windowsize - 1)
  ReDim coefficients.i(windowsize - 1)
  
  ; Hamming-Werte berechnen
  For i = 0 To windowsize - 1
    ; Verwende die Hamming-Formel
    Protected scaledValue.d = 0.54 - 0.46 * Cos((#TwoPI * i) / (windowsize - 1))
    hamming(i) = Int(scaledValue * 100000) ; Skalierung auf Integer
  Next
  
EndProcedure

; FIR-Smooth-Filter für verschiedene Bit-Tiefen
Procedure SmoothFilter(*tempBuffer)
  Protected samples.i = MemorySize(*tempBuffer) / (DefaultBitDepth / 8)
  Protected Dim filtered.i(samples - 1)
  Protected i.i, j.i
  
  ; Normierungsfaktor berechnen
  Protected normalizationFactor.i = 0
  For i = 0 To windowsize - 1
    normalizationFactor + coefficients(i)
  Next
  
  If normalizationFactor = 0
    Debug "Normalization factor is zero in SmoothFilter. Check coefficients."
    ProcedureReturn
  EndIf
  
  ; Normierung der Koeffizienten
  For i = 0 To windowsize - 1
    coefficients(i) = coefficients(i) * 100000 / normalizationFactor
  Next
  
  ; Debug: Normierte Summe prüfen
  Protected totalCoefficients.i = 0
  For i = 0 To windowsize - 1
    totalCoefficients + coefficients(i)
  Next
  Debug "Sum of coefficients (should be 100000): " + Str(totalCoefficients)
  
  ; Filter anwenden
  For i = 0 To samples - 1
    For j = 0 To windowsize - 1
      If (i - j) >= 0 And (i - j) < samples
        Select DefaultBitDepth
          Case 8
            filtered(i) + ((PeekA(*tempBuffer + (i - j)) - 128) * coefficients(j)) / 100000
          Case 16
            filtered(i) + (PeekW(*tempBuffer + (i - j) * 2) * coefficients(j)) / 100000
          Case 24
            ; 24-Bit Wert manuell berechnen
            Protected sample24.i
            sample24 = PeekB(*tempBuffer + (i - j) * 3)
            sample24 | (PeekB(*tempBuffer + (i - j) * 3 + 1) << 8)
            sample24 | (PeekB(*tempBuffer + (i - j) * 3 + 2) << 16)
            If sample24 & $800000
              sample24 - 16777216 ; Negativbereich für 24 Bit
            EndIf
            filtered(i) + (sample24 * coefficients(j)) / 100000
          Case 32
            filtered(i) + (PeekL(*tempBuffer + (i - j) * 4) * coefficients(j)) / 100000
        EndSelect
      EndIf
    Next
  Next
  
  ; Ergebnis zurückschreiben in den Originalbuffer
  For i = 0 To samples - 1
    Protected value.i = filtered(i)
    Select DefaultBitDepth
      Case 8
        If value > 127 : value = 127 : ElseIf value < -128 : value = -128 : EndIf
        PokeA(*tempBuffer + i, value + 128)
      Case 16
        If value > 32767 : value = 32767 : ElseIf value < -32768 : value = -32768 : EndIf
        PokeW(*tempBuffer + i * 2, value)
      Case 24
        If value > 8388607 : value = 8388607 : ElseIf value < -8388608 : value = -8388608 : EndIf
        PokeB(*tempBuffer + i * 3, value & $FF)
        PokeB(*tempBuffer + i * 3 + 1, (value >> 8) & $FF)
        PokeB(*tempBuffer + i * 3 + 2, (value >> 16) & $FF)
      Case 32
        If value > 2147483647 : value = 2147483647 : ElseIf value < -2147483648 : value = -2147483648 : EndIf
        PokeL(*tempBuffer + i * 4, value)
    EndSelect
  Next
EndProcedure

; FIR-Bandpass-Filter
Procedure ApplyBandpassFilter(*tempBuffer, lowFreq.i, highFreq.i)
  Protected i
  
  If lowFreq < 0 Or highFreq < 0 Or lowFreq >= DefaultSampleRate / 2 Or highFreq >= DefaultSampleRate / 2 Or lowFreq >= highFreq
    Debug "Invalid frequency range: lowFreq=" + Str(lowFreq) + ", highFreq=" + Str(highFreq)
    ProcedureReturn
  EndIf
  
  Protected scaledLowFreq.d = lowFreq / DefaultSampleRate
  Protected scaledHighFreq.d = highFreq / DefaultSampleRate
  Protected freqdiff.d = scaledHighFreq - scaledLowFreq
  Protected freqdiff_double.d = freqdiff * 2
  Protected pi_freqdiff = #PI * freqdiff
  Protected half_windowsize = windowsize / 2
  Protected i_minus_windowsize
  
  ; Filterkoeffizienten berechnen
  For i = 0 To windowsize - 1
    i_minus_windowsize = i - half_windowsize
    If i_minus_windowsize <> 0
      coefficients(i) = (2 * Sin(pi_freqdiff * i_minus_windowsize) / i_minus_windowsize) * hamming(i)
    Else
      coefficients(i) = freqdiff_double * hamming(i)
    EndIf
  Next
  
  ; Buffer filtern
  SmoothFilter(*tempBuffer)
EndProcedure

; FIR-Tiefpass-Filter
Procedure ApplyLowpassFilter(*tempBuffer, cutoffFreq.i)
  Protected i.i
  
  ; Filterkoeffizienten berechnen
  Protected scaledCutoff.d = cutoffFreq / DefaultSampleRate ; Grenzfrequenz skalieren
  Protected half_windowsize = windowsize / 2
  Protected double_scaledCutoff.d = scaledCutoff * 2
  Protected pi_scaledCutoff.d = #TwoPI * scaledCutoff
  Protected i_minus_windowsize
  
  For i = 0 To windowsize - 1
    i_minus_windowsize = i - half_windowsize
    If i_minus_windowsize <> 0
      coefficients(i) = (Sin(pi_scaledCutoff * i_minus_windowsize) / i_minus_windowsize) * hamming(i)
    Else
      coefficients(i) = double_scaledCutoff * hamming(i) ; Division durch Null vermeiden
    EndIf
  Next
  
  ; Buffer direkt filtern
  SmoothFilter(*tempBuffer)
EndProcedure

; FIR-Hochpass-Filter
Procedure ApplyHighpassFilter(*tempBuffer, cutoffFreq.i)
  Protected i.i
  
  ; Skalierung der Grenzfrequenz
  Protected scaledCutoff.d = cutoffFreq / DefaultSampleRate
  Protected half_windowsize = windowsize / 2
  Protected double_scaledCutoff.d = 2 * scaledCutoff
  Protected pi_scaledCutoff.d = #TwoPI * scaledCutoff
  Protected i_minus_windowsize
  
  For i = 0 To windowsize - 1
    i_minus_windowsize = i - half_windowsize
    If i_minus_windowsize <> 0
      coefficients(i) = (-Sin(pi_scaledCutoff * i_minus_windowsize) / i_minus_windowsize) * hamming(i)
    Else
      coefficients(i) = (1 - double_scaledCutoff) * hamming(i) ; Division durch Null vermeiden
    EndIf
  Next
  
  ; Buffer direkt filtern
  SmoothFilter(*tempBuffer)
EndProcedure

; FIR-Notch-Filter
Procedure ApplyNotchFilter(*tempBuffer, centerFreq.i, bandwidth.i)
  Protected i.i
  
  ; Skalierung der Frequenzen
  Protected scaledCenterFreq.d = centerFreq / DefaultSampleRate
  Protected scaledBandwidth.d = bandwidth / DefaultSampleRate
  Protected half_windowsize = windowsize / 2
  Protected double_scaledbandwith.d = #TwoPI * scaledBandwidth
  Protected pi_scaledbandwith.d = #PI * scaledBandwidth
  Protected i_minus_windowsize
  
  ; Filterkoeffizienten berechnen
  For i = 0 To windowsize - 1
    i_minus_windowsize = i - half_windowsize
    If i_minus_windowsize <> 0
      coefficients(i) = ((Sin(pi_scaledbandwith * i_minus_windowsize) - Sin(double_scaledbandwith * i_minus_windowsize)) / i_minus_windowsize) * hamming(i)
    Else
      coefficients(i) = (1 - scaledBandwidth) * hamming(i) ; Division durch Null vermeiden
    EndIf
  Next
  
  ; Buffer direkt filtern
  SmoothFilter(*tempBuffer)
EndProcedure

; FIR-Resonanz-Filter
Procedure ApplyResonanceFilter(*tempBuffer, centerFreq.i, resonanceFactor.f)
  Protected i.i
  
  ; Frequenzskalierung
  Protected scaledCenterFreq.d = centerFreq / DefaultSampleRate
  Protected half_windowsize = windowsize / 2
  Protected double_scaledCenterFreq.d = #TwoPI * scaledCenterFreq
  Protected i_minus_windowsize
  
  ; Filterkoeffizienten berechnen
  For i = 0 To windowsize - 1
    i_minus_windowsize = i - half_windowsize
    If i_minus_windowsize <> 0
      coefficients(i) = (Sin(double_scaledCenterFreq * i_minus_windowsize) / i_minus_windowsize) * resonanceFactor * hamming(i)
    Else
      coefficients(i) = resonanceFactor * hamming(i) ; Division durch Null vermeiden
    EndIf
  Next
  
  ; Buffer direkt filtern
  SmoothFilter(*tempBuffer)
EndProcedure

; 3 Band-Equalizer !!
Procedure Apply3BandEqualizer(*tempBuffer, lowFreq.i, lowGain.i, midFreq.i, midGain.i, highFreq.i, highGain.i, lowWidth.i = 100, midWidth.i = 200, highWidth.i = 1000)
  Protected totalSamples.i = MemorySize(*tempBuffer) / (DefaultBitDepth / 8)
  Protected i
  Dim originalBuffer.i(totalSamples - 1) ; Puffer für das Original
  
  ; Originaldaten speichern
  CopyMemory(*tempBuffer, @originalBuffer(), MemorySize(*tempBuffer))
  
  ; Band 1: Low
  Protected scaledLowGain.d = lowGain / 100.0 ; Verstärkungsfaktor skalieren
  ApplyBandpassFilter(*tempBuffer, lowFreq - lowWidth / 2, lowFreq + lowWidth / 2)
  For i = 0 To totalSamples - 1
    Select DefaultBitDepth
      Case 8
        PokeA(*tempBuffer + i, (PeekA(*tempBuffer + i) - 128) * scaledLowGain + 128)
      Case 16
        PokeW(*tempBuffer + i * 2, PeekW(*tempBuffer + i * 2) * scaledLowGain)
      Case 24
        Protected sample24.i = PeekB(*tempBuffer + i * 3) | (PeekB(*tempBuffer + i * 3 + 1) << 8) | (PeekB(*tempBuffer + i * 3 + 2) << 16)
        If sample24 & $800000
          sample24 | -16777216
        EndIf
        sample24 * scaledLowGain
        PokeB(*tempBuffer + i * 3, sample24 & $FF)
        PokeB(*tempBuffer + i * 3 + 1, (sample24 >> 8) & $FF)
        PokeB(*tempBuffer + i * 3 + 2, (sample24 >> 16) & $FF)
      Case 32
        PokeL(*tempBuffer + i * 4, PeekL(*tempBuffer + i * 4) * scaledLowGain)
    EndSelect
  Next
  
  ; Band 2: Mid
  CopyMemory(@originalBuffer(), *tempBuffer, MemorySize(*tempBuffer))
  Protected scaledMidGain.d = midGain / 100.0
  ApplyBandpassFilter(*tempBuffer, midFreq - midWidth / 2, midFreq + midWidth / 2)
  For i = 0 To totalSamples - 1
    Select DefaultBitDepth
      Case 8
        PokeA(*tempBuffer + i, (PeekA(*tempBuffer + i) - 128) * scaledMidGain + 128)
      Case 16
        PokeW(*tempBuffer + i * 2, PeekW(*tempBuffer + i * 2) * scaledMidGain)
      Case 24
        sample24.i = PeekB(*tempBuffer + i * 3) | (PeekB(*tempBuffer + i * 3 + 1) << 8) | (PeekB(*tempBuffer + i * 3 + 2) << 16)
        If sample24 & $800000
          sample24 | -16777216
        EndIf
        sample24 * scaledMidGain
        PokeB(*tempBuffer + i * 3, sample24 & $FF)
        PokeB(*tempBuffer + i * 3 + 1, (sample24 >> 8) & $FF)
        PokeB(*tempBuffer + i * 3 + 2, (sample24 >> 16) & $FF)
      Case 32
        PokeL(*tempBuffer + i * 4, PeekL(*tempBuffer + i * 4) * scaledMidGain)
    EndSelect
  Next
  
  ; Band 3: High
  CopyMemory(@originalBuffer(), *tempBuffer, MemorySize(*tempBuffer))
  Protected scaledHighGain.d = highGain / 100.0
  ApplyBandpassFilter(*tempBuffer, highFreq - highWidth / 2, highFreq + highWidth / 2)
  For i = 0 To totalSamples - 1
    Select DefaultBitDepth
      Case 8
        PokeA(*tempBuffer + i, (PeekA(*tempBuffer + i) - 128) * scaledHighGain + 128)
      Case 16
        PokeW(*tempBuffer + i * 2, PeekW(*tempBuffer + i * 2) * scaledHighGain)
      Case 24
        sample24.i = PeekB(*tempBuffer + i * 3) | (PeekB(*tempBuffer + i * 3 + 1) << 8) | (PeekB(*tempBuffer + i * 3 + 2) << 16)
        If sample24 & $800000
          sample24 | -16777216
        EndIf
        sample24 * scaledHighGain
        PokeB(*tempBuffer + i * 3, sample24 & $FF)
        PokeB(*tempBuffer + i * 3 + 1, (sample24 >> 8) & $FF)
        PokeB(*tempBuffer + i * 3 + 2, (sample24 >> 16) & $FF)
      Case 32
        PokeL(*tempBuffer + i * 4, PeekL(*tempBuffer + i * 4) * scaledHighGain)
    EndSelect
  Next
EndProcedure


Procedure SaveWaveFile(filename.s, *tempBuffer)
  If *tempBuffer = 0
    ProcedureReturn -1
  EndIf
  
  Protected file, datasize, headersize = 44
  datasize = MemorySize(*tempBuffer)
  
  file = CreateFile(#PB_Any, filename)
  If file
    ; RIFF-Header
    WriteString(file, "RIFF", #PB_Ascii)
    WriteLong(file, datasize + headersize - 8)
    WriteString(file, "WAVE", #PB_Ascii)
    
    ; fmt-Chunk
    WriteString(file, "fmt ", #PB_Ascii)
    WriteLong(file, 16)               ; Größe des fmt-Chunks
    WriteWord(file, 1)                ; PCM-Format
    WriteWord(file, DefaultChannels)  ; Kanäle (Mono oder Stereo)
    WriteLong(file, DefaultSampleRate); Abtastrate
    WriteLong(file, DefaultSampleRate * DefaultChannels * (DefaultBitDepth / 8)) ; Byte-Rate
    WriteWord(file, DefaultChannels * (DefaultBitDepth / 8))                     ; Blockgröße
    WriteWord(file, DefaultBitDepth)                                             ; Bits pro Sample
    
    ; data-Chunk
    WriteString(file, "data", #PB_Ascii)
    WriteLong(file, datasize)         ; Größe des Datenbereichs
    WriteData(file, *tempBuffer, datasize)
    
    CloseFile(file)
  EndIf
EndProcedure

Procedure GenerateWhiteNoise(durationMs.i)
  Protected totalSamples = DefaultSampleRate * durationMs / 1000 * DefaultChannels
  Protected *whiteNoise = AllocateMemory(totalSamples * (DefaultBitDepth / 8))
  
  If *whiteNoise
    RandomData(*whitenoise, MemorySize(*whiteNoise))
  EndIf
  
  ProcedureReturn *whiteNoise
EndProcedure

PrecomputeHamming(101) ; the higher, the "better sounding" the filter should work
Define *noise = GenerateWhiteNoise(2000) ; make 2 seconds
ApplyLowpassFilter(*noise, 500)
SaveWaveFile("noise.wav", *noise)

[infratec]

If you need speed (realtime), forget Peek() and Poke().
These are procedures which needs to be called -> needs additional time.

Use something like:

Code: Select all

Structure AsciiArray_Structure
  a.a[0]
EndStructure

Procedure SmoothFilter(*tempBuffer.AsciiArray_Structure)
...
*tempBuffer\a[i - j]

[infratec]

Example:

Code: Select all

Structure ByteArray_Structure
  b.b[0]
EndStructure

Structure WordArray_Structure
  w.w[0]
EndStructure

Structure LongArray_Structure
  l.l[0]
EndStructure


Procedure SmoothFilter1(*tempBuffer)
  Protected samples.i = MemorySize(*tempBuffer) / (DefaultBitDepth / 8)
  Protected Dim filtered.i(samples - 1)
  Protected i.i, j.i
  Protected *b.ByteArray_Structure, *w.WordArray_Structure, *l.LongArray_Structure
  Protected Index.i
  
  *b = *tempBuffer
  *w = *tempBuffer
  *l = *tempBuffer
  
  ; Normierungsfaktor berechnen
  Protected normalizationFactor.i = 0
  For i = 0 To windowsize
    normalizationFactor + coefficients(i)
  Next
  
  If normalizationFactor = 0
    Debug "Normalization factor is zero in SmoothFilter. Check coefficients."
    ProcedureReturn
  EndIf
  
  ; Normierung der Koeffizienten
  For i = 0 To windowsize
    coefficients(i) = coefficients(i) * 100000 / normalizationFactor
  Next
  
  ; Debug: Normierte Summe prüfen
  Protected totalCoefficients.i = 0
  For i = 0 To windowsize
    totalCoefficients + coefficients(i)
  Next
  Debug "Sum of coefficients (should be 100000): " + Str(totalCoefficients)
  
  ; Filter anwenden
  For i = 0 To samples - 1
    For j = 0 To windowsize
      Index = i - j
      If Index >= 0 And Index < samples
        Select DefaultBitDepth
          Case 8
            ;filtered(i) + ((PeekA(*tempBuffer + (i - j)) - 128) * coefficients(j)) / 100000
            filtered(i) + ((*b\b[Index] - 128) * coefficients(j)) / 100000
          Case 16
            ;filtered(i) + (PeekW(*tempBuffer + (i - j) * 2) * coefficients(j)) / 100000
            filtered(i) + (*w\w[Index] * coefficients(j)) / 100000
          Case 24
            ; 24-Bit Wert manuell berechnen
            Protected sample24.i
;             sample24 = PeekB(*tempBuffer + (i - j) * 3)
;             sample24 | (PeekB(*tempBuffer + (i - j) * 3 + 1) << 8)
;             sample24 | (PeekB(*tempBuffer + (i - j) * 3 + 2) << 16)
            Index = Index * 3
            sample24 = *b\b[Index]
            sample24 | (*b\b[Index + 1] << 8)
            sample24 | (*b\b[Index + 2] << 16)
            If sample24 & $800000
              sample24 - 16777216 ; Negativbereich für 24 Bit
            EndIf
            filtered(i) + (sample24 * coefficients(j)) / 100000
          Case 32
            filtered(i) + (*l\l[Index] * coefficients(j)) / 100000
        EndSelect
      EndIf
    Next
  Next
  
  ; Ergebnis zurückschreiben in den Originalbuffer
  For i = 0 To samples - 1
    Protected value.i = filtered(i)
    Select DefaultBitDepth
      Case 8
        If value > 127 : value = 127 : ElseIf value < -128 : value = -128 : EndIf
        ;PokeA(*tempBuffer + i, value + 128)
        *b\b[i] = value + 128
      Case 16
        If value > 32767 : value = 32767 : ElseIf value < -32768 : value = -32768 : EndIf
        ;PokeW(*tempBuffer + i * 2, value)
        *w\w[i] = value
      Case 24
        If value > 8388607 : value = 8388607 : ElseIf value < -8388608 : value = -8388608 : EndIf
;         PokeB(*tempBuffer + i * 3, value & $FF)
;         PokeB(*tempBuffer + i * 3 + 1, (value >> 8) & $FF)
;         PokeB(*tempBuffer + i * 3 + 2, (value >> 16) & $FF)
        Index = i * 3
        *b\b[Index] = value & $FF
        *b\b[Index + 1] = (value >> 8) & $FF
        *b\b[Index + 2] = (value >> 16) & $FF
      Case 32
        If value > 2147483647 : value = 2147483647 : ElseIf value < -2147483648 : value = -2147483648 : EndIf
        ;PokeL(*tempBuffer + i * 4, value)
        *l\l[i] = value
    EndSelect
  Next
EndProcedure

I also reduced window in general by 1 to avoid recalculations in every for loop.


oops...

I just meassured the time ... my version is slower

[jamirokwai]

Hi Infratec,

thanks for the array-hint!

But, how can it be slower than Peek & Poke, I wonder? The overhead should be less.
Maybe the array-positions need to be calculated, while Peek & Poke accesses memory directly?

Example:

Code: Select all

Structure ByteArray_Structure
  b.b[0]
EndStructure

Structure WordArray_Structure
  w.w[0]
EndStructure

Structure LongArray_Structure
  l.l[0]
EndStructure


Procedure SmoothFilter1(*tempBuffer)
  Protected samples.i = MemorySize(*tempBuffer) / (DefaultBitDepth / 8)
  Protected Dim filtered.i(samples - 1)
  Protected i.i, j.i
  Protected *b.ByteArray_Structure, *w.WordArray_Structure, *l.LongArray_Structure
  Protected Index.i
  
  *b = *tempBuffer
  *w = *tempBuffer
  *l = *tempBuffer
  
  ; Normierungsfaktor berechnen
  Protected normalizationFactor.i = 0
  For i = 0 To windowsize
    normalizationFactor + coefficients(i)
  Next
  
  If normalizationFactor = 0
    Debug "Normalization factor is zero in SmoothFilter. Check coefficients."
    ProcedureReturn
  EndIf
  
  ; Normierung der Koeffizienten
  For i = 0 To windowsize
    coefficients(i) = coefficients(i) * 100000 / normalizationFactor
  Next
  
  ; Debug: Normierte Summe prüfen
  Protected totalCoefficients.i = 0
  For i = 0 To windowsize
    totalCoefficients + coefficients(i)
  Next
  Debug "Sum of coefficients (should be 100000): " + Str(totalCoefficients)
  
  ; Filter anwenden
  For i = 0 To samples - 1
    For j = 0 To windowsize
      Index = i - j
      If Index >= 0 And Index < samples
        Select DefaultBitDepth
          Case 8
            ;filtered(i) + ((PeekA(*tempBuffer + (i - j)) - 128) * coefficients(j)) / 100000
            filtered(i) + ((*b\b[Index] - 128) * coefficients(j)) / 100000
          Case 16
            ;filtered(i) + (PeekW(*tempBuffer + (i - j) * 2) * coefficients(j)) / 100000
            filtered(i) + (*w\w[Index] * coefficients(j)) / 100000
          Case 24
            ; 24-Bit Wert manuell berechnen
            Protected sample24.i
;             sample24 = PeekB(*tempBuffer + (i - j) * 3)
;             sample24 | (PeekB(*tempBuffer + (i - j) * 3 + 1) << 8)
;             sample24 | (PeekB(*tempBuffer + (i - j) * 3 + 2) << 16)
            Index = Index * 3
            sample24 = *b\b[Index]
            sample24 | (*b\b[Index + 1] << 8)
            sample24 | (*b\b[Index + 2] << 16)
            If sample24 & $800000
              sample24 - 16777216 ; Negativbereich für 24 Bit
            EndIf
            filtered(i) + (sample24 * coefficients(j)) / 100000
          Case 32
            filtered(i) + (*l\l[Index] * coefficients(j)) / 100000
        EndSelect
      EndIf
    Next
  Next
  
  ; Ergebnis zurückschreiben in den Originalbuffer
  For i = 0 To samples - 1
    Protected value.i = filtered(i)
    Select DefaultBitDepth
      Case 8
        If value > 127 : value = 127 : ElseIf value < -128 : value = -128 : EndIf
        ;PokeA(*tempBuffer + i, value + 128)
        *b\b[i] = value + 128
      Case 16
        If value > 32767 : value = 32767 : ElseIf value < -32768 : value = -32768 : EndIf
        ;PokeW(*tempBuffer + i * 2, value)
        *w\w[i] = value
      Case 24
        If value > 8388607 : value = 8388607 : ElseIf value < -8388608 : value = -8388608 : EndIf
;         PokeB(*tempBuffer + i * 3, value & $FF)
;         PokeB(*tempBuffer + i * 3 + 1, (value >> 8) & $FF)
;         PokeB(*tempBuffer + i * 3 + 2, (value >> 16) & $FF)
        Index = i * 3
        *b\b[Index] = value & $FF
        *b\b[Index + 1] = (value >> 8) & $FF
        *b\b[Index + 2] = (value >> 16) & $FF
      Case 32
        If value > 2147483647 : value = 2147483647 : ElseIf value < -2147483648 : value = -2147483648 : EndIf
        ;PokeL(*tempBuffer + i * 4, value)
        *l\l[i] = value
    EndSelect
  Next
EndProcedure

I also reduced window in general by 1 to avoid recalculations in every for loop.


oops...

I just meassured the time ... my version is slower

[benubi]

This is a glitch of the Assembler backend IMO.

In an other thread I noticed that Lookup tables were SLOWER that Sqr() function itself. But this only happened on assembler backend. With C backend and optimizer the lookup table would be 10-15x faster than Sqr().

I like to use an "Any" structure often to save stack load and variable names

Code: Select all

Structure Any
  StructureUnion
    a.a
    b.b
    c.c
    d.d
    f.f
    i.i
    l.l
    q.q
    u.u
    w.w
;(...)
  EndStructureUnion
EndStructure

I have never used this but the compiler doesn't complain:

Code: Select all

Structure AnyArray
  StructureUnion
    a.a[0]
    b.b[0]
    c.c[0]
    d.d[0]
    f.f[0]
    i.i[0]
    l.l[0]
    q.q[0]
    u.u[0]
    w.w[0]
;(...)
  EndStructureUnion
EndStructure

[jamirokwai]

Hi Benubi,

thanks for the link and the explanation. On MacOS, there is no Asm backend available, AFAIK, only C. I tried with as, words and longs. Seems both are reasonable fast and calculations are off by about 3 to 4 %. That may be inside measurement inaccuracy.

However, if there is a real difference of about 3 %, that *could* impact realtime operations. But people in need of every channel or effect in realtime won't code their software by themselves, would they?

With "Optimize generated code" checked, times are peek, poke 1654, array 1489. That is a difference of about 9 % and it's about 15% faster!

Code: Select all

Structure AudioBuffer_Structure
  a.a[0] ; Für 8-Bit-Audio				; peek, poke 3887, array 3796 ms
  w.w[0] ; Für 16-Bit-Audio				; peek, poke 1928, array 1885 ms
  l.l[0] ; Für 32-Bit-Audio				; peek, poke 974, array 938 ms
EndStructure

*buffer.AudioBuffer_Structure = AllocateMemory(5 * 1024 * 1024)
*bufferpoke = AllocateMemory(5 * 1024 * 1024)

von = ElapsedMilliseconds()
For i = 0 To 49
  For j = 0 To MemorySize(*bufferpoke) - 1 Step 2
    a.w = PeekW(*bufferpoke + j)
    PokeW(*bufferpoke + j, a)
  Next j
Next i
Debug ElapsedMilliseconds() - von

von = ElapsedMilliseconds()
For i = 0 To 49
  For j = 0 To MemorySize(*buffer) - 1 Step 2
    a.w = *buffer\w[j]
    *buffer\w[j] = a
  Next j
Next i
Debug ElapsedMilliseconds() - von

Back to your comment: Yes, pre-calculation of the sinus-values is on the list for my next improvements

This is a glitch of the Assembler backend IMO.

In an other thread I noticed that Lookup tables were SLOWER that Sqr() function itself. But this only happened on assembler backend. With C backend and optimizer the lookup table would be 10-15x faster than Sqr().

I like to use an "Any" structure often to save stack load and variable names

Code: Select all

Structure Any
  StructureUnion
    a.a
    b.b
    c.c
    d.d
    f.f
    i.i
    l.l
    q.q
    u.u
    w.w
;(...)
  EndStructureUnion
EndStructure

I have never used this but the compiler doesn't complain:

Code: Select all

Structure AnyArray
  StructureUnion
    a.a[0]
    b.b[0]
    c.c[0]
    d.d[0]
    f.f[0]
    i.i[0]
    l.l[0]
    q.q[0]
    u.u[0]
    w.w[0]
;(...)
  EndStructureUnion
EndStructure

[Shardik]

On MacOS, there is no Asm backend available, AFAIK, only C.

That's incorrect. Until Ventura the Asm backend is working like a charm. But you have to use pbcompilerasm (in the Compilers folder) instead of pbcompiler. Even in Sonoma and Sequoia the Asm backend is working if you put

Code: Select all

Import "-fno-pie" : EndImport

at the top of your code.

[mk-soft]

You should replace all ASM code in macOS or stop using it. There are (unfortunately) no more new Macs with Intel.

[jamirokwai]

On MacOS, there is no Asm backend available, AFAIK, only C.

That's incorrect. Until Ventura the Asm backend is working like a charm. But you have to use pbcompilerasm (in the Compilers folder) instead of pbcompiler. Even in Sonoma and Sequoia the Asm backend is working if you put

Code: Select all

Import "-fno-pie" : EndImport

at the top of your code.

Thanks for making this clear! My fault.

[jamirokwai]

Hi there,

working on the filter, I got many things working. But this BassBooster is not working as intended. It should enhance low frequencies below the given frequency. Here is a working example, which seems to do something, but not the intended amplification below 400 Hertz.

Code: Select all

Global DefaultSampleRate.i = 44100 ; 8000 to 192000
Global DefaultBitDepth.i = 16 ; or 8
Global DefaultChannels.i = 1  ; 2 = Stereo 

#TwoPI = 2 * #PI

; globale variablen
Global useHammingWindow.b = #False
Global Dim hamming.d(0)      ; Hamming-Fenster (Double)
Global windowsize

Structure EQBand
  centerFreq.i    ; Zentrale Frequenz
  width.i         ; Bandbreite
  gain.f          ; Verstärkung
EndStructure

#filter_06db = 1.0
#filter_12db = 0.5
#filter_24db = 0.25

Procedure MixBuffers(*buffer1, *buffer2)
  ; Berechne die Anzahl der Samples aus der Puffergröße
  Protected size1 = MemorySize(*buffer1)
  Protected size2 = MemorySize(*buffer2)
  Protected samples1 = size1 / (DefaultBitDepth / 8) / DefaultChannels
  Protected samples2 = size2 / (DefaultBitDepth / 8) / DefaultChannels
  
  ; Bestimme die größere Anzahl an Samples
  Protected totalSamples = samples1
  If samples2 > totalSamples
    totalSamples = samples2
  EndIf
  
  ; Speicher für den gemischten Puffer anlegen
  Protected *mixedBuffer = AllocateMemory(totalSamples * DefaultChannels * (DefaultBitDepth / 8))
  If *mixedBuffer = 0
    ProcedureReturn 0 ; Fehler: Kein Speicher verfügbar
  EndIf
  
  ; Mischen der beiden Puffer
  Protected i, sample1.i, sample2.i, mixedSample.i, channel.i
  For i = 0 To totalSamples - 1
    For channel = 0 To DefaultChannels - 1
      If i < samples1
        Select DefaultBitDepth
          Case 8
            sample1 = PeekA(*buffer1 + (i * DefaultChannels + channel)) - 128
          Case 16
            sample1 = PeekW(*buffer1 + (i * DefaultChannels + channel) * 2)
        EndSelect
      Else
        sample1 = 0
      EndIf
      
      If i < samples2
        Select DefaultBitDepth
          Case 8
            sample2 = PeekA(*buffer2 + (i * DefaultChannels + channel)) - 128
          Case 16
            sample2 = PeekW(*buffer2 + (i * DefaultChannels + channel) * 2)
        EndSelect
      Else
        sample2 = 0
      EndIf
      
      ; Addiere die Samples und normalisiere
      mixedSample = (sample1 + sample2) / 2
      Select DefaultBitDepth
        Case 8
          PokeA(*mixedBuffer + (i * DefaultChannels + channel), mixedSample + 128)
        Case 16
          PokeW(*mixedBuffer + (i * DefaultChannels + channel) * 2, mixedSample)
      EndSelect
    Next
  Next
  
  ProcedureReturn *mixedBuffer
EndProcedure

; hamming-fenster vorberechnen
Procedure PrecomputeHamming(newwindowsize.i)
  Protected i
  
  windowsize = newwindowsize
  ReDim hamming.d(windowsize - 1)
  ;   ReDim coefficients.d(windowsize - 1)
  
  ; Hamming-Werte berechnen
  For i = 0 To windowsize - 1
    ; Verwende die Hamming-Formel
    hamming(i) = 0.54 - 0.46 * Cos((#TwoPI * i) / (windowsize - 1))
  Next
  
EndProcedure

; FIR-Smooth-Filter für verschiedene Bit-Tiefen
Procedure SmoothFilter(*tempBuffer, *coefficients)
  Protected samples.i = MemorySize(*tempBuffer) / (DefaultBitDepth / 8) / DefaultChannels
  Protected Dim filtered.d(samples - 1)
  Protected i.i, j.i, channel.i
  Protected coefficient.d
  
  Select DefaultBitDepth
    Case 8
      ; Filter anwenden für 8-Bit
      For channel = 0 To DefaultChannels - 1
        For j = 0 To windowsize - 1
          coefficient = PeekD(*coefficients + j * SizeOf(Double))
          For i = j To samples - 1
            filtered(i) + ((PeekA(*tempBuffer + ((i - j) * DefaultChannels + channel)) - 128) * coefficient)
          Next
        Next
        
        ; Ergebnis zurückschreiben in den Originalbuffer
        For i = 0 To samples - 1
          Protected value.d = filtered(i)
          If value > 127 : value = 127 : ElseIf value < -128 : value = -128 : EndIf
          PokeA(*tempBuffer + (i * DefaultChannels + channel), value + 128)
        Next
      Next
      
    Case 16
      ; Filter anwenden für 16-Bit
      For channel = 0 To DefaultChannels - 1
        For j = 0 To windowsize - 1
          coefficient = PeekD(*coefficients + j * SizeOf(Double))
          For i = j To samples - 1
            filtered(i) + (PeekW(*tempBuffer + ((i - j) * DefaultChannels + channel) * 2) * coefficient)
          Next
        Next
        
        ; Ergebnis zurückschreiben in den Originalbuffer
        For i = 0 To samples - 1
          value.d = filtered(i)
          If value > 32767 : value = 32767 : ElseIf value < -32768 : value = -32768 : EndIf
          PokeW(*tempBuffer + (i * DefaultChannels + channel) * 2, value)
        Next
      Next
  EndSelect
EndProcedure

Procedure ApplyLowpassFilter(*tempBuffer, cutoffFreq.d, attenuation.d = #filter_24db)
  Protected i.i
  
  ; Filterkoeffizienten berechnen
  Protected scaledCutoff.d = cutoffFreq / DefaultSampleRate ; Grenzfrequenz skalieren
  Protected half_windowsize = windowsize / 2
  Protected double_scaledCutoff.d = scaledCutoff * 2.0
  Protected pi_scaledCutoff.d = #TwoPI * scaledCutoff
  Protected i_minus_windowsize
  Protected Dim coefficients.d(windowsize) ; Filterkoeffizienten (Double)
  
  ; Filterkoeffizienten berechnen
  For i = 0 To windowsize - 1
    i_minus_windowsize = i - half_windowsize
    If i_minus_windowsize <> 0
      coefficients(i) = ((Sin(pi_scaledCutoff * i_minus_windowsize) / i_minus_windowsize) * hamming(i)) * attenuation
    Else
      coefficients(i) = (double_scaledCutoff * hamming(i)) * attenuation ; Division durch Null vermeiden
    EndIf
  Next
  
  ; Filter anwenden
  SmoothFilter(*tempBuffer, @coefficients())
EndProcedure

; Bass-Boost !!
Procedure ApplyBassBoost(*tempBuffer, boost.i, cutoffFreq.i = 200)
  Protected scaledBoost.d = boost / 100.0
  
  ; Kopie des Original-Buffers erstellen
  Define *filteredBuffer = AllocateMemory(MemorySize(*tempBuffer))
  If *filteredBuffer = 0
    Debug "Memory allocation failed for filteredBuffer."
    ProcedureReturn
  EndIf
  CopyMemory(*tempBuffer, *filteredBuffer, MemorySize(*tempBuffer))
  
  ; Tiefpass-Filter auf die Kopie anwenden
  ApplyLowpassFilter(*filteredBuffer, cutoffFreq)
  
  ; Gefilterten Puffer mit dem Original mischen
  Define *mixedBuffer = MixBuffers(*tempBuffer, *filteredBuffer)
  
  If *mixedBuffer = 0
    Debug "MixBuffers failed."
    FreeMemory(*filteredBuffer)
    ProcedureReturn
  EndIf
  
  ; Ergebnis zurückkopieren und temporäre Puffer freigeben
  CopyMemory(*mixedBuffer, *tempBuffer, MemorySize(*tempBuffer))
  FreeMemory(*filteredBuffer)
  FreeMemory(*mixedBuffer)
EndProcedure

;- init
PrecomputeHamming(101) ; Standard ist 101, 63 wäre ausreichend für realtime, 31 wäre gut für realtime, aber schlechtere Qualität


Procedure SaveWaveFile(filename.s, *tempBuffer)
  If *tempBuffer = 0
    ProcedureReturn -1
  EndIf
  
  Protected file, datasize, headersize = 44
  datasize = MemorySize(*tempBuffer)
  
  file = CreateFile(#PB_Any, filename)
  If file
    ; RIFF-Header
    WriteString(file, "RIFF", #PB_Ascii)
    WriteLong(file, datasize + headersize - 8)
    WriteString(file, "WAVE", #PB_Ascii)
    
    ; fmt-Chunk
    WriteString(file, "fmt ", #PB_Ascii)
    WriteLong(file, 16)               ; Größe des fmt-Chunks
    WriteWord(file, 1)                ; PCM-Format
    WriteWord(file, DefaultChannels)  ; Kanäle (Mono oder Stereo)
    WriteLong(file, DefaultSampleRate); Abtastrate
    WriteLong(file, DefaultSampleRate * DefaultChannels * (DefaultBitDepth / 8)) ; Byte-Rate
    WriteWord(file, DefaultChannels * (DefaultBitDepth / 8))                     ; Blockgröße
    WriteWord(file, DefaultBitDepth)                                             ; Bits pro Sample
    
    ; data-Chunk
    WriteString(file, "data", #PB_Ascii)
    WriteLong(file, datasize)         ; Größe des Datenbereichs
    WriteData(file, *tempBuffer, datasize)
    
    CloseFile(file)
  EndIf
EndProcedure

Procedure GenerateWhiteNoise(durationMs.i)
  Protected totalSamples = DefaultSampleRate * durationMs / 1000 * DefaultChannels
  Protected *whiteNoise = AllocateMemory(totalSamples * (DefaultBitDepth / 8))
  
  If *whiteNoise
    RandomData(*whitenoise, MemorySize(*whiteNoise))
  EndIf
  
  ProcedureReturn *whiteNoise
EndProcedure

PrecomputeHamming(101) ; the higher, the "better sounding" the filter should work
Define *noise = GenerateWhiteNoise(2000) ; make 2 seconds
ApplyBassBoost(*noise, 4, 400)
SaveWaveFile("noise.wav", *noise)

[SMaag]

For Speed up:

1. Calculations in integer do not make sense if it's the basic function is a float.
Integer base calculation will speed up your code if you use SSE Registers, which are much faster in Integer!

coefficients(i) = coefficients(i) * 100000 / normalizationFactor
Better do this all as float, in this case you can eliminate the division what will speed up your code a lot.

2. The selection of the BitDeepth
better do the BitDeepth selection outside the For Loops and use a seperate Loop Code for each Depth.

Code: Select all

Select BitDepth
  Case 8
     For I = 1 to windowsize
     	For J
     	Next
     Next

  Case   16
     For I = 1 to windowsize
     	For J
     	Next
     Next
  
  Case 24
     For I = 1 to windowsize
     	For J
     	Next
     Next