Digital Signal Processing #1

The following example uses sequence and array operators to implement a simple digital filter. The example builds on the previous WAV I/O example. A digital filter frequently needs multiply-add operators, these may be implemented using Array.map2

> let fx = [|for x in {1..5} do yield sin((float)x)|];; 

val fx : float [] =
  [|0.8414709848; 0.9092974268; 0.1411200081; -0.7568024953; -0.9589242747|]

> let filter = [|1.0; -1.0; 0.0; 0.0; 0.0|];;

val filter : float [] = [|1.0; -1.0; 0.0; 0.0; 0.0|]

> Array.map2 (fun i j -> i*j) fx filter;;
val it : float [] = [|0.8414709848; -0.9092974268; 0.0; 0.0; 0.0|]

Full source:

open System
open System.IO
open System.Text

// https://ccrma.stanford.edu/courses/422/projects/WaveFormat/

type WavAudioFormat = WavPCM
exception WavUnknownFormat
exception WavUnknownNumChannels
exception WavUnknownBitsPerSample
exception WavUnknownFilter

let Sample16BitLE2Float i = float(BitConverter.ToInt16((Seq.toArray i), 0))/32768.0
let Float2Sample16BitLE s = BitConverter.GetBytes(int16(s*32768.0))

let breakBy windowSize (source: seq<_>) =
    seq { let arr = Array.create windowSize (Seq.nth 0 source)
          let i = ref 0
          use e = source.GetEnumerator()
          while e.MoveNext() do
              arr.[!i] <- e.Current
              i := !i + 1
              if !i = windowSize then
                yield Array.init windowSize (fun j -> arr.[j])
                i := 0
        }

type WavFile() =
  let mutable AudioFormat:WavAudioFormat = WavPCM
  let mutable NumChannels:int = 0
  let mutable SampleRate:int = 0
  let mutable BitsPerSample:int = 0
  let mutable Samples:seq<float> = seq [||]

  member this.Read(filename:string) =
    let a = File.ReadAllBytes(filename)
    if a.Length < 40 then
      raise(WavUnknownFormat)
    if a.[ 0.. 3] <> Encoding.ASCII.GetBytes("RIFF")
    || a.[ 8..11] <> Encoding.ASCII.GetBytes("WAVE")
    || a.[12..15] <> Encoding.ASCII.GetBytes("fmt ") then
      raise(WavUnknownFormat)

    let ChunkSize    = BitConverter.ToInt32(a.[4..7], 0)
    let SubChunk1Size = BitConverter.ToInt32(a.[16..19], 0)
    let AudioFormatVal = BitConverter.ToInt16(a.[20..21], 0)
    AudioFormat <- match AudioFormatVal with
                   | 1s -> WavPCM
                   | _ -> raise(WavUnknownFormat)

    NumChannels <- int(BitConverter.ToInt16(a.[22..23], 0))
    if NumChannels < 1 || NumChannels > 2 then raise(WavUnknownNumChannels)

    SampleRate <- int(BitConverter.ToInt32(a.[24..27], 0))

    BitsPerSample <- int(BitConverter.ToInt16(a.[34..35], 0))
    if BitsPerSample % 8 <> 0 then raise(WavUnknownBitsPerSample)

    let SubChunk2Start = 12 + 8 + SubChunk1Size
    if a.[SubChunk2Start..(SubChunk2Start+3)] <>
       Encoding.ASCII.GetBytes("data") then
      raise(WavUnknownFormat)

    let SubChunk2Size = BitConverter.ToInt32(
                          a.[(SubChunk2Start+4)..(SubChunk2Start+7)], 0)
    let data = a.[(SubChunk2Start+8)..(SubChunk2Start+8+SubChunk2Size-1)]
    Samples <- match BitsPerSample with
               | 16 -> data |> breakBy 2 |> Seq.map Sample16BitLE2Float
               | _ -> raise(WavUnknownBitsPerSample)

  member this.Write(filename:string) =
    let NumSamples = Seq.length Samples

    let mutable bytes = Encoding.ASCII.GetBytes "RIFF"
    bytes <- Array.append bytes (BitConverter.GetBytes
      (44 + (NumSamples * BitsPerSample/8)))
    bytes <- Array.append bytes (Encoding.ASCII.GetBytes "WAVE")

    bytes <- Array.append bytes (Encoding.ASCII.GetBytes "fmt ")
    match AudioFormat with
    | WavPCM ->
      bytes <- Array.append bytes (BitConverter.GetBytes 16)
      bytes <- Array.append bytes (BitConverter.GetBytes (int16 1))
      bytes <- Array.append bytes (BitConverter.GetBytes (int16 NumChannels))
      bytes <- Array.append bytes (BitConverter.GetBytes SampleRate)
      bytes <- Array.append bytes (BitConverter.GetBytes
        (SampleRate*NumChannels*BitsPerSample/8))
      bytes <- Array.append bytes (BitConverter.GetBytes
        (int16 (NumChannels*BitsPerSample/8)))
      bytes <- Array.append bytes (BitConverter.GetBytes (int16 BitsPerSample))

      bytes <- Array.append bytes (Encoding.ASCII.GetBytes "data")
      bytes <- Array.append bytes (BitConverter.GetBytes
        (NumSamples * NumChannels * BitsPerSample/8))

      let data = Samples |> Seq.map Float2Sample16BitLE
                 |> Seq.collect (fun i -> i) |> Seq.toArray
      bytes <- Array.append bytes data

    File.WriteAllBytes(filename, bytes)

  member this.Print():unit =
    printfn "RIFF ChunkSize= ";

    printf "WAVE AudioFormat=";
    if AudioFormat = WavPCM then printf "PCM "
    else printf "unknown "

    printfn "NumChannels=%d SampleRate=%d BitsPerSample=%d "
      NumChannels SampleRate BitsPerSample

    let NumSamples = Seq.length Samples
    printfn "%d samples, %dmin %ds %dms playtime" NumSamples
      (NumSamples/SampleRate/60)
      ((NumSamples/SampleRate)%60)
      (((NumSamples%SampleRate)*1000)/SampleRate)

  member this.GetSamples() = Samples

  member this.LowPass(samples) =
    // Simple sliding average window filter
    let window_size = 32
    let xcoeff = [| for i in 1..window_size do yield 1.0/float(window_size) |]
    let buffer = ref [|for i in 1..(Array.length xcoeff) do yield 0.0|]
    Seq.map (fun s ->
      buffer := Array.append (!buffer).[1..((Array.length !buffer)-1)] [|s|]
      Array.map2 (fun x b -> x*b) xcoeff !buffer |> Array.sum
      ) samples

  member this.Filter(filter):unit =
    match filter with
    | "none" -> ()
    | "lowpass" ->  Samples <- this.LowPass(Samples)
    | _ -> raise WavUnknownFilter

let args = System.Environment.GetCommandLineArgs()
if (Array.length args) = 4 then
  let filter = args.[1]
  let infilename = args.[2]
  let outfilename = args.[3]
  let w = new WavFile()
  w.Read(infilename)
  w.Print()
  w.Filter(filter)
  w.Write(outfilename)
else
  printf "wav2.exe <filter:none,sinc> <input filename.wav> <output filename.wav>\n"

Invocation:

~> fsc.exe wav2.fs
Microsoft (R) F# 2.0 Compiler build 2.0.0.0
Copyright (c) Microsoft Corporation. All Rights Reserved.
~> time ./wav2.exe lowpass /usr/share/sounds/alsa/test.wav out.wav
RIFF ChunkSize= 
WAVE AudioFormat=PCM NumChannels=2 SampleRate=44100 BitsPerSample=16 
670320 samples, 0min 15s 200ms playtime

real    0m4.253s
user    0m4.221s
sys     0m0.095s

Further reading:

Advertisements
This entry was posted in Algorithms, IO. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s