Added firls

Designs FIR filter taps using the least squares method

src/Filters/Filters.jl

@@ -29,7 +29,8 @@ export  FilterType,
         digitalfilter,
         iirnotch,
         kaiserord,
-        FIRWindow
+        FIRWindow,
+        firls
 
 include("response.jl")
 export  freqs,

src/Filters/design.jl

@@ -472,6 +472,34 @@ function firprototype(n::Integer, ftype::Bandstop)
     out
 end
 
+
+#
+# http://cnx.org/contents/eb1ecb35-03a9-4610-ba87-41cd771c95f2@7/Linear-Phase_Fir_Filter_Design
+#
+function firls(n::Integer, bands::AbstractVector)
+    isodd(n) || throw(ArgumentError("Number of taps must be odd"))
+
+    m = int((n-1)/2)
+    Q = zeros(m+1,m+1)
+    b = zeros(m+1,m+1)
+
+    for (ƒl, ƒh, amplitude) in bands
+        for j in 0:m
+
+            b[j+1] += amplitude*2π*(ƒh*sinc(ƒh*j) - ƒl*sinc(ƒl*j))
+
+            for i in 0:m
+                Q[i+1,j+1] += π*( ƒh*sinc(ƒh*(i-j)) - ƒl*sinc(ƒl*(i-j)) ) # Q₁, Toeplitz
+                Q[i+1,j+1] += π*( ƒh*sinc(ƒh*(i+j)) - ƒl*sinc(ƒl*(i+j)) ) # Q₂, Hankel
+            end
+        end
+    end
+
+    a = Q \ b
+    h = [a[m+1:-1:2]./2, a[1], a[2:m+1]./2]
+end
+
+
 scalefactor(coefs::Vector, ::Union(Lowpass, Bandstop)) = sum(coefs)
 function scalefactor(coefs::Vector, ::Highpass)
     c = zero(coefs[1])

test/firls.jl

@@ -0,0 +1,67 @@
+using DSP
+
+# Lowpass
+N    = 41
+b    = [(0.0, 0.4, 1), (0.6, 1, 0)]
+h_lp = firls(N, b)
+
+# Highpass
+N    = 41
+b    = [(0.0, 0.4, 0), (0.6, 1, 1)]
+h_hp = firls(N, b)
+
+# Bandpass
+N    = 41
+b    = [(0, .25, -50dBa), (.4, .6, 0dBa), (.75, 1, -50dBa)]
+h_bp = firls(N, b)
+
+# Bandstop
+N    = 41
+b    = [(0, .2, 0dBa), (.4, .6, -10dBa), (.8, 1, 0dBa)]
+h_bs = firls(N, b)
+
+using PyPlot
+
+function Filters.freqz(h::Vector, w = linspace(0, π, 250))
+    pr = PolynomialRatio(h, [one(eltype(h))])
+    freqz(pr, w)
+end
+Util.amp2db(x::Complex) = amp2db(abs(x))
+Util.amp2db(x::AbstractVector) = [amp2db(x) for x in x]
+
+clf()
+subplot(2,2,1)
+title("Lowpass")
+ylabel("dB")
+xlabel("Normalized ƒ (π rad/s)")
+H = amp2db(freqz(h_lp))
+f = linspace(0,1, length(H))
+plot(f, H)
+grid()
+
+subplot(2,2,2)
+title("Highpass")
+ylabel("dB")
+xlabel("Normalized ƒ (π rad/s)")
+H = amp2db(freqz(h_hp))
+f = linspace(0,1, length(H))
+plot(f, H)
+grid()
+
+subplot(2,2,3)
+title("Bandpass")
+ylabel("dB")
+xlabel("Normalized ƒ (π rad/s)")
+H = amp2db(freqz(h_bp))
+f = linspace(0,1, length(H))
+plot(f, H)
+grid()
+
+subplot(2,2,4)
+title("Bandstop")
+ylabel("dB")
+xlabel("Normalized ƒ (π rad/s)")
+H = amp2db(freqz(h_bs))
+f = linspace(0,1, length(H))
+plot(f, H)
+grid()