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 href="http://strugglingthroughproblems.blogspot.com/2011/06/robjecttables-are-awesome.html">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.