Merge pull request #111 from JuliaDSP/jsk/resamplers

resample function

Author: simonster

doc/filters.rst

@@ -132,6 +132,20 @@ Filter application
     Apply FIR filter ``b`` along the first dimension of array ``x``
     using an FFT-based overlap-save algorithm.
 
+.. function:: resample(x, rate[, coef])
+
+    Resample ``x`` at rational or arbitrary ``rate``.
+    ``coef`` is an optional vector of FIR filter taps. If ``coef``
+    is not provided, the taps will be computed using a Kaiser window.
+
+    Internally, ``resample`` uses a polyphase ``FIRFilter`` object,
+    but performs additional operations to make resampling a signal easier.
+    It compensates for for the ``FIRFilter``'s delay (ramp-up), and appends
+    zeros to ``x``. The result is that when the input and output signals
+    are plotted on top of each other, they correlate very well, but one
+    signal will have more samples that the other.
+
+
 
 Filter design
 -------------

src/Filters/Filters.jl

@@ -29,7 +29,8 @@ export  FilterType,
         digitalfilter,
         iirnotch,
         kaiserord,
-        FIRWindow
+        FIRWindow,
+        resample_filter
 
 include("response.jl")
 export  freqs,
@@ -42,6 +43,9 @@ include("stream_filt.jl")
 export  FIRFilter,
         outputlength,
         inputlength,
-        reset!
+        reset!,
+        resample,
+        setphase!,
+        timedelay
 
 end

src/Filters/design.jl

@@ -423,9 +423,7 @@ function kaiserord(transitionwidth::Real, attenuation::Real=60)
         β = 0.0
     end
 
-    α = π*β
-
-    return n, α
+    return n, β/π
 end
 
 immutable FIRWindow{T}
@@ -496,3 +494,42 @@ function digitalfilter(ftype::FilterType, proto::FIRWindow)
     out = coefs .* proto.window
     proto.scale ? scale!(out, 1/scalefactor(out, ftype)) : out
 end
+
+
+# Compute FIR coefficients necessary for arbitrary rate resampling
+function resample_filter(rate::FloatingPoint, Nϕ = 32, rel_bw = 1.0, attenuation = 60)
+    f_nyq       = rate >= 1.0 ? 1.0/Nϕ : rate/Nϕ
+    cutoff      = f_nyq * rel_bw
+    trans_width = cutoff * 0.2
+
+    # Determine resampling filter order
+    hLen, α = kaiserord(trans_width, attenuation)
+
+    # Round the number of taps up to a multiple of Nϕ.
+    # Otherwise the missing taps will be filled with 0.
+    hLen = Nϕ * ceil(Int, hLen/Nϕ)
+
+    # Design filter
+    h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(hLen, α)))
+    scale!(h, Nϕ)
+end
+
+# Compute FIR coefficients necessary for rational rate resampling
+function resample_filter(rate::Rational, rel_bw = 1.0, attenuation = 60)
+    Nϕ          = num(rate)
+    decimation  = den(rate)
+    f_nyq       = min(1/Nϕ, 1/decimation)
+    cutoff      = f_nyq * rel_bw
+    trans_width = cutoff * 0.2
+
+    # Determine resampling filter order
+    hLen, α = kaiserord(trans_width, attenuation)
+
+    # Round the number of taps up to a multiple of Nϕ (same as interpolation factor).
+    # Otherwise the missing taps will be filled with 0.
+    hLen = Nϕ * ceil(Int, hLen/Nϕ)
+
+    # Design filter
+    h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(hLen, α)))
+    scale!(h, Nϕ)
+end