News: Sign up for "The Monthly Mean," the newsletter that dares to call itself average, www.pmean.com/news.
| P.Mean: Simple examples of wavelets (created 2013-12-07).
News: Sign up for "The Monthly Mean," the newsletter that dares to call itself average, www.pmean.com/news. |
|
Wavelets represent a method to decompose a signal or an image. It represents an alternative to the Fourier decomposition. It is useful in image compression, edge analysis and flexible de-noising of even complicated discontinuous functions. Although there are adaptations to more complex settings, wavelets work best when analyzing a signal that is sampled at evenly spaced time points.
To understand wavelets well, you need to refresh yourself on how the Fourier decomposition really works. Think about Fourier as a way to travel to Trigland, the land of sines and cosines.
If you are curious, I wrote a few simple functions to display information nicely.
comp <- function(by.hand, using.r, lab) {
out <- as.matrix(c(by.hand, using.r))
dimnames(out) <- list(c("By hand", "Using R"), lab)
out
}
dilute <- function(x) {
n <- length(x)
y <- rep(0, 2 * n)
y[2 * (1:n) - 1] <- x
return(y)
}
display.formula <- function(form, cx = 2) {
par(mar = rep(0, 4))
plot(0:1, 0:1, axes = FALSE, type = "n")
text(0.1, 0.5, form, adj = 0, cex = cx)
}
dot.prod <- function(signal, filter, i = 1) {
padded.signal <- c(signal, signal)
m <- length(filter)
j <- i:(i + m - 1)
sum(padded.signal[j] * filter)
}
show.axis <- function() {
par(mar = c(0.1, 0.1, 0.1, 0.1))
plot(c(0, 2400), c(0, 1), axes = FALSE, type = "n")
axis(side = 1, at = c(0, 512, 1024, 1536, 2048), pos = 1)
}
show.single <- function(squiggle, lab, y.range = range(squiggle)) {
par(mar = c(0.1, 0.1, 0.1, 0.1), yaxs = "i")
plot(c(0, 2400), y.range, axes = FALSE, type = "n")
lines(1:2048, squiggle)
text(2100, squiggle[2048], lab, adj = 0)
}
There is data set with the wavelets package in R called “ecg” that represents 2048 evenly spaced values from a typical electrocardiogram. You can use this file to apply both a Fourier transform and a wavelet transform. Here's the code to get, display, and graph that signal.
library("wavelets")
data(ecg)
head(ecg) # print the first few data points
## [1] -0.10477 -0.09314 -0.08150 -0.11641 -0.08150 -0.11641
tail(ecg) # print the last few data points
## [1] -0.4306 -0.4422 -0.4771 -0.4422 -0.4771 -0.4771
show.single(ecg, "ecg\nsignal")
You can ignore the units of the signal itself, as they are arbitrary (at least for this application), but the minimum value is -1.48 and the maximum is 1.44. You compute the Fourier decomposition using the fft command, which is part of the base R package. The letters stand for “fast fourier transform”.
ft.ecg <- fft(ecg)/2048
head(ft.ecg)
## [1] -0.040860+0.000000i -0.144005-0.128781i -0.004312+0.007707i
## [4] 0.023643+0.046245i 0.004964-0.026598i -0.010242+0.017355i
tail(ft.ecg)
## [1] 0.010344+0.000403i -0.010242-0.017355i 0.004964+0.026598i
## [4] 0.023643-0.046245i -0.004312-0.007707i -0.144005+0.128781i
ift.ecg <- fft(ft.ecg, inverse = TRUE)
head(ift.ecg)
## [1] -0.10477-0i -0.09314+0i -0.08150-0i -0.11641+0i -0.08150-0i -0.11641+0i
tail(ift.ecg)
## [1] -0.4306+0i -0.4422+0i -0.4771+0i -0.4422-0i -0.4771-0i -0.4771-0i
The object ft.ecg contains the Fourier decomposition, and ift.ecg is the inverse of the decomposition. You need to scale the Fourier decomposition by the number of data points (2,048).
The inverse gets the original signal (-0.1048, -0.09314, -0.0815 …) except that there is a +0i or -0i tacked on to the end.
The complex numbers in the Fourier decomposition are a short-hand for the mix of sines and cosines. The first number has only a real component and no imaginary component. The real component represent the mean of the entire signal.
comp(mean(ecg), Re(ft.ecg[1]), "ft[1]")
## ft[1]
## By hand -0.04086
## Using R -0.04086
For the remaining numbers, the real part represents the cosine function and the imaginary part represents the sine function. So
Here's a plot of the cosine portion,
t.unit <- (1:2048)/2048
cos1 <- +Re(ft.ecg[2]) * cos(2 * pi * t.unit)
show.single(cos1, "cosine\nfunction", c(-0.35, 0.35))
the sine portion,
sin1 <- -Im(ft.ecg[2]) * sin(2 * pi * t.unit)
show.single(sin1, "sine\nfunction", c(-0.35, 0.35))
and the sum.
show.single(sin1 + cos1, "Sum of\nsine and\ncosine", c(-0.35, 0.35))
The graph of the cosine function shows a peak halfway through (1024) and a trough at the end (or equivalently at the beginning). This is actually flipped around from the normal cosine function, but recall that the multiplier front of the cosine function is negative.
The graph of the sine function shows a peak at ¼ of the the way through (512) and a trough ¾ of the way through (1536). Notice that the coefficient of the sine function is actually the negative of the imaginary portion. This is one of those scaling issues that will drive you nuts.
The sum shows the peak for the sum that is somewhere between the peaks of the sine and cosine functions. It's actually slightly closer to the peak of the cosine function, but not by much. The trough is also between the two troughs. The relative weights given to the sine and cosine functions and the positive or negative values associated with these weights will determine where the peak and trough occur for the sum.
If you are familiar with the arctangent function, you can compute the location of peak/trough. The Mod function will tell you the amplitude.
ft.ratio <- -Im(ft.ecg[2])/Re(ft.ecg[2])
comp(round(2048 * (atan(ft.ratio) + pi)/(2 * pi)), which(sum1 == max(sum1)),
"Peak")
## Peak
## By hand 786
## Using R 786
comp(Mod(ft.ecg[2]), max(sum1), "Amplitude")
## Amplitude
## By hand 0.1932
## Using R 0.1932
Now this was a lot of work to decode the second component of the Fourier decomposition. You might forget to scale by the length of the signal. You could reverse the sines and cosines. You might compute the sines and cosines on a range from -pi to pi instead of from 0 to 2pi. You might leave out an important minus sign. Trust me, I did all these things.
There's a faster and more convenient way. Let the computer do the work for you. Take the Fourier decomposition, zero out everything except the second component, and then invert the decomposition. This allows you to graph the second component directly.
ft.comp2 <- rep(0, length(ft.ecg))
ft.comp2[2] <- ft.ecg[2]
comp2 <- Re(fft(ft.comp2, inverse = TRUE))
show.single(comp2, "Second\nFourier\ncomponent", c(-0.35, 0.35))
You can use the same method to show the shape of the third Fourier component.
ft.comp3 <- rep(0, length(ft.ecg))
ft.comp3[3] <- ft.ecg[3]
comp3 <- Re(fft(ft.comp3, inverse = TRUE))
mx.3 <- max(comp3)
show.single(comp3, "Third\nFourier\ncomponent", 1.8 * mx.3 * c(-1, 1))
This component has a higher frequency, with two peaks rather than one. Do the same for the fourth Fourier component
ft.comp4 <- rep(0, length(ft.ecg))
ft.comp4[4] <- ft.ecg[4]
comp4 <- Re(fft(ft.comp4, inverse = TRUE))
mx.4 <- max(comp4)
show.single(comp4, "Fourth\nFourier\ncomponent", 1.8 * mx.4 * c(-1, 1))
The fourth Fourier component has three peaks. The interesting thing is that the second and third and fourth Fourier components are orthogonal to one another.
dot.prod(comp2, comp3)
## [1] -5.59e-20
dot.prod(comp2, comp4)
## [1] 3.725e-16
dot.prod(comp3, comp4)
## [1] 3.515e-20
Actually, every Fourier component is orthogonal to every other Fourier component. Mathematical folks like “orthogonal” because it allows you to derive the mathematical properties of the decomposition with greater ease. For a practical person, what “orthogonal” means is that the information contained in one Fourier component does not “overlap” with that in any other Fourier component. It also means you can run the decomposition in any order (low frequency to high frequency, high frequency to low frequency, or even a haphazard frequency). Finally, orthogonality also preserves independence properties in certain settings.
The Fourier decomposition continues with components of higher and higher frequency. The following code computes the fiftieth component of the Fourier transform.
ft.comp50 <- rep(0, length(ft.ecg))
ft.comp50[50] <- ft.ecg[50]
comp50 <- Re(fft(ft.comp50, inverse = TRUE))
mx.50 <- max(comp50)
mx <- 0.35 * mx.50/mx.2
show.single(comp50, "Fiftieth\nFourier\ncomponent", c(-mx, mx))
Here is the graph.
Here's the fiftieth Fourier component, and if you count carefully, you will get 49 peaks. Count them carefully. No, don't bother, I'll do it for you. Because that are packed so closely together, I will onlyshow you every third peak.
show.single(comp50, "Fiftieth\nFourier\ncomponent", c(-mx, mx))
text(9 + 3 * (0:16) * 2048/49, 0.9 * mx, 3 * (0:16) + 1)
If you continue down the list, you get a higher and higher frequency or equivalently more and more peaks. At about the midpoint, something strange happens.
ft.comp1025 <- rep(0, length(ft.ecg))
ft.comp1025[1025] <- ft.ecg[1025]
comp1025 <- Re(fft(ft.comp1025, inverse = TRUE))
mx.1025 <- max(comp1025)
mx <- 0.35 * mx.1025/mx.2
show.single(comp1025, "1,025th\nFourier\ncomponent", c(-mx, mx))
The 1,025th component, which has 1,024 peaks, looks like a black blob. But if you plot just the points, it looks different.
par(mar = rep(0, 4))
plot(c(0, 2400), c(-mx, mx), type = "n", axes = FALSE)
points(1:2048, comp1025)
Notice that there are just two values. This sort of makes sense. There are 1,024 peaks and 1,024 troughs and that doesn't leave any room for anything in the middle. It gets even stranger when you look at the 1,026th component.
ft.comp1026 <- rep(0, length(ft.ecg))
ft.comp1026[1026] <- ft.ecg[1026]
comp1026 <- Re(fft(ft.comp1026, inverse = TRUE))
mx.1026 <- max(comp1026)
mx <- 0.35 * mx.1026/mx.2
show.single(comp1026, "1,026th\nFourier\ncomponent", c(-mx, mx))
When you plot the points, you see two competing waves.
par(mar = rep(0, 4))
plot(c(0, 2400), c(-mx, mx), type = "n", axes = FALSE)
points(1:2048, comp1026)
Also, the coefficients start pairing together with earlier coefficients.
ft.ecg[c(1024, 1026)]
## [1] 0.000379+0.001086i 0.000379-0.001086i
ft.ecg[c(1023, 1027)]
## [1] 0.001134+0.000438i 0.001134-0.000438i
ft.ecg[c(1022, 1028)]
## [1] 0.000819-0.001452i 0.000819+0.001452i
The very last coefficients in the Fourier transform match up with the very first coefficients.
ft.ecg[c(2, 2048)]
## [1] -0.144-0.1288i -0.144+0.1288i
ft.ecg[c(3, 2047)]
## [1] -0.004312+0.007707i -0.004312-0.007707i
ft.ecg[c(4, 2046)]
## [1] 0.02364+0.04625i 0.02364-0.04625i
This relates to the concept of aliasing. Wikipedia has a nice explanation about aliasing that is worth reading.
This page was written by
Steve Simon and is licensed under the Creative
Commons Attribution 3.0 United States License. Need more
information? I have a page with general help
resources. You can also browse for pages similar to this one at .R Software