mirror of https://github.com/mpastell/Weave.jl
Added a test for the script reader
parent
2b5e8ee968
commit
d2a298a1d6
|
@ -0,0 +1,79 @@
|
||||||
|
#' % 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.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 Gadfly, DSP
|
||||||
|
|
||||||
|
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, Geom.line,
|
||||||
|
Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Magnitude (db)"))
|
||||||
|
|
||||||
|
#' And again with default options
|
||||||
|
|
||||||
|
h_phase = unwrap(-atan2(imag(h),real(h)))
|
||||||
|
plot(y = h_phase, x = ws, Geom.line,
|
||||||
|
Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Phase (radians)"))
|
|
@ -0,0 +1,178 @@
|
||||||
|
% 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.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:
|
||||||
|
|
||||||
|
|
||||||
|
% 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.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:
|
||||||
|
|
||||||
|
|
||||||
|
~~~~{.julia}
|
||||||
|
using Gadfly, DSP
|
||||||
|
|
||||||
|
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$
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
## 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$
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
~~~~{.julia}
|
||||||
|
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.
|
||||||
|
|
||||||
|
|
||||||
|
## Plot the frequency and impulse response
|
||||||
|
|
||||||
|
The next code chunk is executed in term mode, see the [script](FIR_design.jl) for syntax.
|
||||||
|
|
||||||
|
~~~~{.julia}
|
||||||
|
julia>
|
||||||
|
h_db = log10(abs(h))
|
||||||
|
|
||||||
|
1024-element Array{Float32,1}:
|
||||||
|
0.0
|
||||||
|
-1.52727e-6
|
||||||
|
-6.03147e-6
|
||||||
|
-1.35386e-5
|
||||||
|
-2.39452e-5
|
||||||
|
-3.71738e-5
|
||||||
|
-5.31213e-5
|
||||||
|
-7.16583e-5
|
||||||
|
-9.25781e-5
|
||||||
|
-0.000115803
|
||||||
|
⋮
|
||||||
|
-3.21659
|
||||||
|
-3.18342
|
||||||
|
-3.15605
|
||||||
|
-3.13379
|
||||||
|
-3.11615
|
||||||
|
-3.10275
|
||||||
|
-3.09335
|
||||||
|
-3.08778
|
||||||
|
-3.08593
|
||||||
|
julia> ws = w/pi*(fs/2)
|
||||||
|
|
||||||
|
linspace(0.0,10.0,1024)
|
||||||
|
julia> plot(y = h_db, x = ws, Geom.line,
|
||||||
|
Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Magnitude (db)"))
|
||||||
|
~~~~~~~~~~~~~
|
||||||
|
|
||||||
|
|
||||||
|
![](figures/FIR_design_3_1.png)\
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
And again with default options
|
||||||
|
|
||||||
|
|
||||||
|
~~~~{.julia}
|
||||||
|
h_phase = unwrap(-atan2(imag(h),real(h)))
|
||||||
|
plot(y = h_phase, x = ws, Geom.line,
|
||||||
|
Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Phase (radians)"))
|
||||||
|
~~~~~~~~~~~~~
|
||||||
|
|
||||||
|
|
||||||
|
![](figures/FIR_design_4_1.png)\
|
||||||
|
|
|
@ -49,3 +49,8 @@ weave("documents/gadfly_markdown_test.jmd", doctype="github",plotlib="gadfly", i
|
||||||
result = @compat readstring(open("documents/gadfly_markdown_test.md"))
|
result = @compat readstring(open("documents/gadfly_markdown_test.md"))
|
||||||
ref = @compat readstring(open("documents/gadfly_markdown_test_ref.md"))
|
ref = @compat readstring(open("documents/gadfly_markdown_test_ref.md"))
|
||||||
@test result == ref
|
@test result == ref
|
||||||
|
|
||||||
|
weave("documents/FIR_design.jl", doctype="pandoc", plotlib="gadfly", informat="script")
|
||||||
|
result = @compat readstring(open("documents/FIR_design.md"))
|
||||||
|
ref = @compat readstring(open("documents/FIR_design_ref.md"))
|
||||||
|
@test result == ref
|
||||||
|
|
Loading…
Reference in New Issue