IIR equalizer biquad filters

README.md

@@ -3,13 +3,12 @@ DSP.jl
 
 [![Build Status](https://travis-ci.org/JuliaDSP/DSP.jl.svg?branch=master)](https://travis-ci.org/JuliaDSP/DSP.jl)
 [![Coverage Status](https://coveralls.io/repos/JuliaDSP/DSP.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaDSP/DSP.jl?branch=master)
-[![Package Evaluator](http://iainnz.github.io/packages.julialang.org/badges/DSP_release.svg)](http://iainnz.github.io/packages.julialang.org/?pkg=DSP&ver=release)
-[![Package Evaluator](http://iainnz.github.io/packages.julialang.org/badges/DSP_nightly.svg)](http://iainnz.github.io/packages.julialang.org/?pkg=DSP&ver=nightly)
+[![Package Evaluator](http://pkg.julialang.org/badges/DSP_0.3.svg)](http://pkg.julialang.org/?pkg=DSP)
+[![Package Evaluator](http://pkg.julialang.org/badges/DSP_0.4.svg)](http://pkg.julialang.org/?pkg=DSP)
 
 DSP.jl provides a number of common DSP routines in Julia.  So far, the following are implemented:
 
 - [Periodogram estimation](http://dspjl.readthedocs.org/en/latest/periodograms.html)
 - [Window functions](http://dspjl.readthedocs.org/en/latest/windows.html)
 - [Filter design and filtering](http://dspjl.readthedocs.org/en/latest/filters.html)
 - [Utility functions](http://dspjl.readthedocs.org/en/latest/util.html)
-

doc/filters.rst

@@ -174,6 +174,15 @@ Such filters are documented below.
     bandwidth ``bandwidth``. If ``fs`` is not specified, ``Wn`` is
     interpreted as a normalized frequency in half-cycles/sample.
 
+.. function:: iirEQ(response, Wn, gain[, Q, fs])
+
+    Design typical IIR equalizer biquad filters. Available response types
+    are 'lowpass', 'highpass, 'bandpass', 'bandstop', 'peak', 'highshelf',
+    'lowshelf'. Optional parameter ``Q`` does not have any effect on shelving
+    filter designs and causes a resonance on low/highpass filters. If ``fs``
+    is not specified, ``Wn`` is interpreted as a normalized frequency in
+    half-cycles/sample.
+
 .. _response-types:
 
 Filter response types

src/Filters/Filters.jl

@@ -30,6 +30,7 @@ export  FilterType,
         analogfilter,
         digitalfilter,
         iirnotch,
+        iirEQ,
         kaiserord,
         FIRWindow,
         resample_filter

src/Filters/design.jl

@@ -406,6 +406,103 @@ function iirnotch(w::Real, bandwidth::Real; fs=2)
     Biquad(b, -2b*cosw0, b, -2b*cosw0, 2b-1)
 end
 
+# Design of biquad EQ filters based on formula from the book
+# DAFx, Udo Zoelzer, Wiley 2011
+# Available filter types are: 'lowpass', 'highpass, 'bandpass', 'bandstop',
+#                             'peak', 'highshelf', 'lowshelf'
+function iirEQ(filttype::AbstractString, w::Real, gain::Real, Q::Real=1.0/sqrt(2.0), fs::Real=2)
+    w = normalize_freq(w, fs)
+
+    K = tan(pi * w);
+    N = K*K*Q + K + Q;
+    V0 = 10^(gain/20)
+    a1 = (2*Q*(K*K-1)) / N;
+    a2 = (K*K*Q-K+Q) / N;
+
+    if filttype == "lowpass"
+        b0 = K*K*Q / N * V0
+        b1 = 2*K*K*Q / N * V0
+        b2 = b0
+    elseif filttype == "highpass"
+        b0 = Q / N * V0
+        b1 = -2*Q / N * V0
+        b2 = b0
+    elseif filttype == "bandpass"
+        b0 = K / N * V0
+        b1 = 0
+        b2 = -b0
+    elseif filttype == "bandstop"
+        b0 = param*(1+K*K) / N
+        b1 = 2*param*(K*K-1) / N
+        b2 = b0
+    elseif filttype == "peak"
+        K2 = K*K
+
+        if gain<=0
+            denom = 1 + 1/(V0*Q)*K + K2
+            b0 = (1 + 1/Q*K + K2) / denom
+            b1 = 2*(K2-1) / denom
+            b2 = (1 - 1/Q*K + K2) / denom
+
+            a1 = 2*(K2 - 1) / denom
+            a2 = (1 - 1/(V0*Q)*K + K2) / denom
+        else
+            denom = (1 + 1/Q*K + K2)
+            b0 = (1 + V0/Q*K + K2) / denom
+            b1 = 2*(K2-1) / denom
+            b2 = (1 - V0/Q*K + K2) / denom
+
+            a1 = 2*(K2 - 1) / denom
+            a2 = (1 - 1/Q*K + K2) / denom
+        end
+    elseif filttype == "highshelf"
+        K2 = K*K
+
+        if gain<=0
+            denom = (1 + sqrt(2*V0)*K + V0*K2);
+            b0 = V0*(1 + sqrt(2)*K + K2) / denom;
+            b1 = 2*V0*(K2-1) / denom;
+            b2 = V0*(1 - sqrt(2)*K + K2) / denom;
+
+            a1 = 2*(V0*K2 - 1) / denom;
+            a2 = (1 - sqrt(2*V0)*K + V0*K2) / denom;
+        else
+            denom = (1 + sqrt(2)*K + K2);
+            b0 = (V0 + sqrt(2*V0)*K + K2) / denom;
+            b1 = 2*(K2-V0) / denom;
+            b2 = (V0 - sqrt(2*V0)*K + K2) / denom;
+
+            a1 = 2*(K2 - 1) / denom;
+            a2 = (1 - sqrt(2)*K + K2) / denom;
+        end
+    elseif filttype == "lowshelf"
+        K2 = K*K;
+        V0 = 10^(gain/20);
+
+        if gain<=0
+            denom = (V0 + sqrt(2*V0)*K + K2);
+            b0 = V0*(1 + sqrt(2)*K + K2) / denom;
+            b1 = 2*V0*(K2-1) / denom;
+            b2 = V0*(1 - sqrt(2)*K + K2) / denom;
+
+            a1 = 2*(K2 - V0) / denom;
+            a2 = (V0 - sqrt(2*V0)*K + K2) / denom;
+        else
+            denom = (1 + sqrt(2)*K + K2);
+            b0 = (1 + sqrt(2*V0)*K + V0*K2) / denom;
+            b1 = 2*(V0*K2-1) / denom;
+            b2 = (1 - sqrt(2*V0)*K + V0*K2) / denom;
+
+            a1 = 2*(K2 - 1) / denom;
+            a2 = (1 - sqrt(2)*K + K2) / denom;
+        end
+    else
+        throw(ArgumentError("Unknown filter type: $filttype"))
+    end
+
+    Biquad(b0, b1, b2, a1, a2)
+end
+
 #
 # FIR filter design
 #