Reverse Iteration

< !-- begin{Schunk} !--> < !--\end{Schunk}!-->

Time to horrify some people.

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

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> source(“pretend.R”)

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

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:

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

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

< !-- ccc -->

> x = rnorm(length(t))

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!--> < !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

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

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

And filter it:

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> X = rfft(x)

< !-- ccc -->

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

< !-- ccc -->

> s = 2*pi*1i*f

< !-- ccc -->

> B = .2

< !-- ccc -->

> f0 = .5

< !-- ccc -->

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

< !-- ccc -->

> Y = H*X

< !-- ccc -->

> y = irfft(Y)

< !-- ccc -->

> min(f)

< !-- ccc --> < !-- end{Sinput} !-->

[1] 0

< !-- begin{Sinput} !-->

> max(f)

< !-- ccc --> < !-- end{Sinput} !-->

[1] 1.95

< !--\end{Schunk}!-->

So that it looks like:

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

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

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

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

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> B = .1

< !-- ccc -->

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

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

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

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> y.a = analytic(y)

< !-- ccc -->

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

< !-- ccc -->

> f0.est = qr.solve(

> cbind(t, 1),

> cbind(phase)

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

< !-- ccc -->

> f0.est

< !-- ccc --> < !-- end{Sinput} !-->

[1] 0.4818572

< !--\end{Schunk}!-->

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

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

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

< !-- ccc -->

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

> B < < - B.

> f0.est

> })

< !-- ccc -->

> f0.ests

< !-- ccc --> < !-- end{Sinput} !-->

[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

< !--\end{Schunk}!--> < !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> plot(Bs, f0.ests)

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

And make the same plot for a different center:

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> f0 = 1.0

< !-- ccc -->

> f0.ests

< !-- ccc --> < !-- end{Sinput} !-->

[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

< !-- begin{Sinput} !-->

> plot(Bs, f0.ests)

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

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

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> y.a = analytic(y)

< !-- ccc -->

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

< !-- ccc -->

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

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

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

< !-- begin{Schunk} !--> < !-- begin{Sinput} !-->

> analytic = function(u) {

> n = length(u)

> m = n/2 + 1

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

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

> fft(U, inv = T)

> }

< !-- ccc --> < !--\end{Sinput}!--> < !--\end{Schunk}!-->

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.