<commit_before>#' % FIR filter design with Julia
#' % <NAME>
#' % 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.

#' ## Functions for frequency, phase, impulse and step response

#' Let's first define functions to plot filter
#' properties.

using Gadfly, DSP

plot(x=1:10)

#' ## Lowpass FIR filter

#' Designing a lowpass FIR filter is very simple to do with SciPy, all you
#' need to do is to define the window length, cut off frequency and the
#' window.

#' The Hamming window is defined as:
#' $w(n) = \alpha - \beta\cos\frac{2\pi n}{N-1}$, where $\alpha=0.54$ and $\beta=0.46$

#' The next code chunk is executed in term mode, see the [Python script](FIR_design.py) for syntax.
#' Notice also that Pweave can now catch multiple figures/code chunk.

#+ term=true
n = 61



#' ## Highpass FIR Filter

#' Let's define a highpass FIR filter, if you compare to original blog
#' post you'll notice that it has become easier since 2009. You don't
#' need to do ' spectral inversion "manually" anymore!

n = 101

#' ## Bandpass FIR filter

#' Notice that the plot has a caption defined in code chunk options.

#+ caption = "Bandpass FIR filter."

n = 1001<commit_msg>Added example script for publishing<commit_after>#' % FIR filter design with Julia
#' % <NAME>
#' % 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.

#' The Hamming window is defined as:
#' $w(n) = \alpha - \beta\cos\frac{2\pi n}{N-1}$, where $\alpha=0.54$ and $\beta=0.46$

#' The next code chunk is executed in term mode, see the [script](FIR_design.jl) for syntax.
#+ term=true


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

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)"))

h_phase = unwrap(-atan2(imag(h),real(h)))
plot(y = h_phase, x = ws, Geom.line,
    Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Phase (radians)"))
{'fraction_non_alphanumeric': 0.12408272181454302, 'fraction_numerical': 0.016010673782521682, 'mean_word_length': 3.7851063829787233, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '883191', 'n_tokens_mistral': 1585, 'n_tokens_neox': 1386, 'n_words': 601}