From fd1f13eb33215396e4ebcbdaaf4c068831264655 Mon Sep 17 00:00:00 2001 From: Matti Pastell Date: Sun, 30 Oct 2016 18:50:02 +0200 Subject: [PATCH] Update examples for Julia 0.5, #54 --- examples/FIR_design_plots.jl | 80 ++++++++++++++++++++++++++++++++++++ examples/julia_sample.mdw | 2 +- 2 files changed, 81 insertions(+), 1 deletion(-) create mode 100644 examples/FIR_design_plots.jl diff --git a/examples/FIR_design_plots.jl b/examples/FIR_design_plots.jl new file mode 100644 index 0000000..6bfa39b --- /dev/null +++ b/examples/FIR_design_plots.jl @@ -0,0 +1,80 @@ +#' % FIR filter design with Julia +#' % Matti Pastell +#' % 21th April 2016 + +#' # Introduction + +#' This an example of a julia script that can be published using +#' [Weave](http://mpastell.github.io/Weave.jl/latest/usage/). +#' The script can be executed normally using Julia +#' or published to HTML or pdf with Weave. +#' Text is written in markdown in lines starting with "`#'` " and code +#' is executed and results are included in the published document. + +#' Notice that you don't need to define chunk options, but you can using +#' `#+`. just before code e.g. `#+ term=True, caption='Fancy plots.'`. +#' If you're viewing the published version have a look at the +#' [source](FIR_design_plots.jl) to see the markup. + + +#' # FIR Filter Design + +#' We'll implement lowpass, highpass and ' bandpass FIR filters. If +#' you want to read more about DSP I highly recommend [The Scientist +#' and Engineer's Guide to Digital Signal +#' Processing](http://www.dspguide.com/) which is freely available +#' online. + +#' ## Calculating frequency response + +#' DSP.jl package doesn't (yet) have a method to calculate the +#' the frequency response of a FIR filter so we define it: + +using Plots, DSP +plotlyjs() + +function FIRfreqz(b::Array, w = linspace(0, π, 1024)) + n = length(w) + h = Array{Complex64}(n) + sw = 0 + for i = 1:n + for j = 1:length(b) + sw += b[j]*exp(-im*w[i])^-j + end + h[i] = sw + sw = 0 + end + return h +end + + +#' ## Design Lowpass FIR filter + +#' Designing a lowpass FIR filter is very simple to do with DSP.jl, all you +#' need to do is to define the window length, cut off frequency and the +#' window. We will define a lowpass filter with cut off frequency at 5Hz for a signal +#' sampled at 20 Hz. +#' We will use the Hamming window, which is defined as: +#' $w(n) = \alpha - \beta\cos\frac{2\pi n}{N-1}$, where $\alpha=0.54$ and $\beta=0.46$ + + +fs = 20 +f = digitalfilter(Lowpass(5, fs = fs), FIRWindow(hamming(61))) +w = linspace(0, pi, 1024) +h = FIRfreqz(f, w) + +#' ## Plot the frequency and impulse response + +#' The next code chunk is executed in term mode, see the [script](FIR_design.jl) for syntax. +#+ term=true + +h_db = log10(abs(h)) +ws = w/pi*(fs/2) +plot(y = h_db, x = ws, + xlabel = "Frequency (Hz)", ylabel = "Magnitude (db)") + +#' And again with default options + +h_phase = unwrap(-atan2(imag(h),real(h))) +plot(y = h_phase, x = ws, + xlabel = "Frequency (Hz)", ylabel = "Phase (radians)") diff --git a/examples/julia_sample.mdw b/examples/julia_sample.mdw index dbee48d..4ac2163 100644 --- a/examples/julia_sample.mdw +++ b/examples/julia_sample.mdw @@ -17,7 +17,7 @@ weave("examples/julia_sample.mdw") <>= x = 1:10 -d = {"Weave" => "testing"} +d = Dict("Weave" => "testing") y = [2, 4 ,8] @