Reverse Iteration

July 3, 2011
By

(This article was first published on Struggling Through Problems, and kindly contributed to R-bloggers)

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.

To leave a comment for the author, please follow the link and comment on his blog: Struggling Through Problems.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.