Reverse Iteration

July 3, 2011
By

[This article was first published on Struggling Through Problems, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 their blog: Struggling Through Problems.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.



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.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)