Time to horrify some people.

First let’s include the code we wrote last time,

> source(“pretend.R”)

and the dependency-tracking environment it creates will be used to run all the following examples.

Let’s look at, I don’t know, I’m just trying to demonstrate a language feature so uh… band-pass filtering Gaussian noise.

Here’s some noise:

> t = seq(0, 10, len=40)

> x = rnorm(length(t))

> plot(t, x, type=‘l’)

And filter it:

> X = rfft(x)

> f = 1/length(t)/(t[2]–t[1]) * (seq_along(X)–1)

> s = 2*pi*1i*f

> B = .2

> f0 = .5

> H = B/(s**2 + B*s + (2*pi*f0)**2)

> Y = H*X

> y = irfft(Y)

> min(f)

[1] 0

> max(f)

[1] 1.95

So that it looks like:

> plot(t, y, type=‘l’)

We’ve made it look wavier. We could decrease the bandwidth and make it cleaner:

> B = .1

> plot(t, y, type=‘l’)

We could unrwap the phase and find the (what we already know to be) the frequency of the wave we just made:

> y.a = analytic(y)

> phase = unwrap(Arg(y.a))

> f0.est = qr.solve(

> cbind(t, 1),

> cbind(phase)

> )[[1]] / (2*pi)

> f0.est

[1] 0.4818572

We can try various bandwidths and see what estimates we get for the center:

> Bs = seq(0.1, 0.3, len=12)

> f0.ests = sapply(Bs, function(B.) {

> B <<- B.

> f0.est

> })

> f0.ests

[1] 0.4818572 0.4814956 0.4811829 0.4809197 0.4807073 0.4805482 0.4804464
[8] 0.4804080 0.4804399 0.4805393 0.4806539 0.3708171

> plot(Bs, f0.ests)

And make the same plot for a different center:

> f0 = 1.0

> f0.ests

[1] 0.9406669 0.9400220 0.9393042 0.9385253 0.9376873 0.9367848 0.9358192
[8] 0.9348181 0.9338298 0.9328914 0.9320139 0.8671673

> plot(Bs, f0.ests)

Now the interesting thing about all that we’ve done here is that nowhere is the interesting code wrapped in either functions or loops: it’s not wrapped in anything at all. There are loops, but they occurr after. Code such as

> y.a = analytic(y)

> phase = unwrap(Arg(y.a))

> f0.est = qr.solve(cbind(t, 1), cbind(phase))[[1]]/(2 * pi)

was written entirely without repetition in mind.

I have to admit that I did cheat a little bit. If you go back and look at the code you’ll notice that each variable is assigned exactly once. Hence the hiding of logic like

> analytic = function(u) {

> n = length(u)

> m = n/2 + 1

> U = fft(u)/length(u)

> U[(m + 1):n] = 0

> fft(U, inv = T)

> }

into functions, because this is awkward to write without assigning U twice. Because in the main program each variable is only assigned once, when a modification is made, it is trivial to decide which other variables must be updated and how.