Merge pull request #18 from rprechelt/biquads

Initial Implementation of Biquads

src/DSP.jl

@@ -11,6 +11,18 @@ type DFTSetup{T}
     end
 end
 
+type Biquad{T}
+    setup::Ptr{Void}
+    sections::Int
+
+    function Biquad(setup::Ptr{Void}, sections::Int)
+        biquadsetup = new(setup, sections)
+        finalizer(biquadsetup, biquaddestroy)
+        biquadsetup
+    end
+end
+
+## === FUNCTIONS == ##
 
 ## === FUNCTIONS == ##
 for (T, suff) in ((Float64, "D"), (Float32, ""))
@@ -62,7 +74,6 @@ for (T, suff) in ((Float64, "D"), (Float32, ""))
         end
     end
 
-
     """
     Cross-correlation of two Vector{T}'s 'X' and 'Y'.
     Result vector should have at least length(X) + length(Y) - 1 elements
@@ -121,9 +132,70 @@ for (T, suff) in ((Float64, "D"), (Float32, ""))
             xcorr(X, X)
         end
     end
-
 end
 
+## == Biquadratic/IIR filtering
+for (T, suff) in ((Float64, "D"), )
+
+    """
+    Initializes a vDSP_biquad_setup for use with vDSP_biquad. A multi-section filter
+    can be initialized with a single call to biquad_create_setup. coefficients must
+    contain 5 coefficients for each section. The three feed-forward coefficients are
+    specified first, followed by the two feedback coefficients. Returns a Biquad object.
+
+    Returns: Biquad{T}
+    """
+    @eval begin
+        function biquadcreate(coefficients::Vector{$T},  sections::Int)
+            if length(coefficients) < 5*sections
+                error("Incomplete biquad specification provided - coefficients must
+                            contain 5 elements for each filter section")
+            end
+            setup = ccall(($(string("vDSP_biquad_CreateSetup", suff), libacc)),  Ptr{Void},
+                          (Ptr{$T}, UInt64),
+                          coefficients, sections)
+            return Biquad{$T}(setup, sections)
+        end
+    end
+
+
+    """
+    Filters an input array X with the specified Biquad filter and filter delay values provided
+    in delays; only numelem elements are filtered. After execution, delays contains the final
+    state data of the filter.
+
+    Returns: Vector{T}
+    """
+    @eval begin
+        function biquad(X::Vector{$T}, delays::Vector{$T}, numelem::Int, biquad::Biquad{$T})
+            if length(delays) < (2*(biquad.sections)+2)
+                error("Incomplete delay specification provided - delays must contain 2*M + 2
+                                values where M is the number of sections in the biquad")
+            end
+            result::Vector{$T} = similar(X)
+            ccall(($(string("vDSP_biquad", suff), libacc)),  Void,
+                  (Ptr{Void},  Ptr{$T},  Ptr{$T},  Int64,  Ptr{$T},  Int64, UInt64),
+                  biquad.setup, delays, X,  1, result,  1,  numelem)
+            return result
+        end
+    end
+
+
+    """
+    Frees all resources associated with a particular Biquad previously
+    created through a call to biquad_create_setup. This is called automatically
+    when the setup object is no longer visible to the garbage collector.
+
+    Returns: Void
+    """
+    @eval begin
+        function biquaddestroy(biquad::Biquad{$T})
+            ccall(($(string("vDSP_biquad_DestroySetup", suff), libacc)),  Void,
+                  (Ptr{Void}, ),
+                  biquad.setup)
+        end
+    end
+end
 
 ## == WINDOW GENERATION == ##
 for (T, suff) in ((Float32, ""), (Float64, "D"))

test/REQUIRE

@@ -1 +1,2 @@
 BaseTestNext
+DSP

test/runtests.jl

@@ -1,4 +1,5 @@
 using AppleAccelerate
+using DSP
 
 if VERSION >= v"0.5-"
     using Base.Test
@@ -121,8 +122,8 @@ end
 
 for T in (Float32,  Float64)
     @testset "Convolution & Correlation::$T" begin
-        X::Array{T} = randn(N)
-        Y::Array{T} = randn(N)
+        X::Vector{T} = randn(N)
+        Y::Vector{T} = randn(N)
         @testset "Testing $f::$T" for f in [:conv, :xcorr]
             @eval fb = $f
             @eval fa = AppleAccelerate.$f
@@ -137,6 +138,19 @@ for T in (Float32,  Float64)
     end
 end
 
+for T in (Float64, )
+    @testset "Biquadratic Flitering::$T" begin
+        @testset "Single Section::$T" begin
+            X::Vector{T} = randn(10)
+            d::Vector{T} = zeros(4)
+            c::Vector{T} = [x%0.5 for x in randn(5)]
+            fdsp = DSP.Biquad(c[1], c[2], c[3], c[4], c[5])
+            fa = AppleAccelerate.biquadcreate(c, 1)
+            @test DSP.filt(fdsp, X) ≈ AppleAccelerate.biquad(X, d, length(X), fa)
+        end
+    end
+end
+
 
 for T in (Float32, Float64)
     @testset "Window Functions::$T" begin