**John Myles White: Die Sudelbücher » Statistics**, and kindly contributed to R-bloggers)

This post is a continuation of my series dealing with matrix operations for image processing. My next goal is to demonstrate the construction of simple low-pass and high-pass spatial frequency filters in R. It’s easy enough to construct simple versions of these filters in R using the Fast Fourier Transform (also known as the FFT), but, because the FFT is a slightly complicated tool, I’m going to build up to using it progressively over a few posts.

For starters, I want to review the use of complex numbers in R. As always, if you’re interested in reviewing the mathematics of complex numbers, I’d start by browsing references online. The tutorial here seemed good to me at first glance, though I can’t claim to have read it through.

Because complex numbers are implemented in the “base” package, it’s very easy to start working with them. To construct the complex number *x + iy*, you use `complex`

:

1 2 3 4 5 6 7 |
x <- 1 y <- 1 z <- complex(real = x, imaginary = y) z # [1] 1+1i |

It’s conventional in mathematics to use *z* to refer to a complex number, so I’ll continue on with that tradition.

As always occurs with mathematical data types in R, you can convert other objects to class “complex” using `as.complex`

:

1 2 |
as.complex(-1) # [1] -1+0i |

And you can test that an object is complex using `is.complex`

:

1 2 |
is.complex(as.complex(-1)) # [1] TRUE |

Beyond those standard operations, there are five essential mathematical operations you’d want to use on complex numbers.

First off, you want to be able to extract the real and imaginary components of a complex number. You can do this using `Re`

and `Im`

respectively:

1 2 3 4 5 6 7 |
z <- complex(real = 2, imaginary = 1) Re(z) # [1] 2 Im(z) # [1] 1 |

The *(x, y)* representation of numbers is easier to understand at first, but a polar coordinates representation is often more practical. You can get the relevant components of this representation by finding the modulus and complex argument of a complex number. In R, you would use `Mod`

and `Arg`

:

1 2 3 4 5 6 7 8 9 10 |
z <- complex(real = 0, imaginary = 1) Mod(z) # [1] 1 Arg(z) # [1] 1.570796 pi / 2 # [1] 1.570796 |

This corresponds to the intuition that *i* should be at a distance 1 from the origin and an angle of *pi / 2*.

Finally, you’ll want to be able to take the complex conjugate of a complex number; to do that in R, you can use `Conj`

:

1 2 3 4 5 |
Conj(z) # [1] 0-1i Mod(z) == z * Conj(z) # [1] TRUE |

As you can see, the modulus of *z* equals z times the conjugate of z, which is exactly what you expect.

Now, historically, complex numbers were invented so that you could find the square root of negative numbers. By default, `sqrt`

does not return a complex number when you ask for the square root of a negative number. Instead, it produces a `NaN`

error, as you can see below:

1 2 3 4 |
sqrt(-1) # [1] NaN # Warning message: # In sqrt(-1) : NaNs produced |

To get the complex square root, you need to cast your negative number as a complex number using `as.complex`

before applying `sqrt`

:

1 2 3 |
sqrt(as.complex(-1)) # [1] 0+1i |

In practice, complex numbers can be used to solve a huge range of problems, not least of which is the efficient representation of the FFT of an input vector. For fun, I thought I’d show a simple application: building a fractal using the Mandelbrot set.

To quote Wikipedia,

Mathematically the Mandelbrot set can be defined as the set of complex values of

cfor which the orbit of 0 under iteration of the complex quadratic polynomialz_{n+1}=z_{n}^{2}+cremains bounded. That is, a complex number,c, is in the Mandelbrot set if, when starting withz_{0}= 0 and applying the iteration repeatedly, the absolute value ofz_{n}never exceeds a certain number (that number depends onc) however largengets.

We can translate this definition into R pretty easily by making certain assumptions about the exactness we want in our results. For my purposes, I’ll say that a number is in the Mandelbrot set if, after 100 iterations of the Mandelbrot algorithm, it’s never once had a modulus greater than 1,000,000. As you’ll see from the image I produced, this assumption works out pretty well, though it takes some real number crunching to get results out of it for reasonable sized data sets.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
in.mandelbrot.set <- function(c, iterations = 100, bound = 1000000) { z <- 0 for (i in 1:iterations) { z <- z ** 2 + c if (Mod(z) > bound) { return(FALSE) } } return(TRUE) } |

Once we have this algorithm, we can easily generate an image of the Mandelbrot set by producing a matrix of complex numbers and depicting their inclusion in the set using `image`

:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
resolution <- 0.001 sequence <- seq(-1, 1, by = resolution) m <- matrix(nrow = length(sequence), ncol = length(sequence)) for (x in sequence) { for (y in sequence) { mandelbrot <- in.mandelbrot.set(complex(real = x, imaginary = y)) m[round((x + resolution + 1) / resolution), round((y + resolution + 1) / resolution)] <- mandelbrot } } png('mandelbrot.png') image(m) dev.off() |

And here’s the resulting image:

UPDATE: Thanks to Giuseppe for pointing out the rounding error in the original code, producing an image with lines in it.

**leave a comment**for the author, please follow the link and comment on their blog:

**John Myles White: Die Sudelbücher » Statistics**.

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...