# Posit AI Weblog: Wavelet Remodel

Be aware: Like a number of prior ones, this put up is an excerpt from the forthcoming guide, Deep Studying and Scientific Computing with R torch. And like many excerpts, it’s a product of exhausting trade-offs. For added depth and extra examples, I’ve to ask you to please seek the advice of the guide.

## Wavelets and the Wavelet Remodel

What are wavelets? Just like the Fourier foundation, they’re capabilities; however they don’t lengthen infinitely. As a substitute, they’re localized in time: Away from the middle, they shortly decay to zero. Along with a location parameter, additionally they have a scale: At completely different scales, they seem squished or stretched. Squished, they may do higher at detecting excessive frequencies; the converse applies after they’re stretched out in time.

The essential operation concerned within the Wavelet Remodel is convolution – have the (flipped) wavelet slide over the info, computing a sequence of dot merchandise. This manner, the wavelet is mainly in search of similarity.

As to the wavelet capabilities themselves, there are numerous of them. In a sensible software, we’d wish to experiment and decide the one which works finest for the given knowledge. In comparison with the DFT and spectrograms, extra experimentation tends to be concerned in wavelet evaluation.

The subject of wavelets could be very completely different from that of Fourier transforms in different respects, as properly. Notably, there’s a lot much less standardization in terminology, use of symbols, and precise practices. On this introduction, I’m leaning closely on one particular exposition, the one in Arnt Vistnes’ very good guide on waves . In different phrases, each terminology and examples mirror the alternatives made in that guide.

## Introducing the Morlet wavelet

The Morlet, also called Gabor, wavelet is outlined like so:

[
Psi_{omega_{a},K,t_{k}}(t_n) = (e^{-i omega_{a} (t_n – t_k)} – e^{-K^2}) e^{- omega_a^2 (t_n – t_k )^2 /(2K )^2}
]

This formulation pertains to discretized knowledge, the sorts of information we work with in follow. Thus, (t_k) and (t_n) designate deadlines, or equivalently, particular person time-series samples.

This equation seems to be daunting at first, however we will “tame” it a bit by analyzing its construction, and pointing to the primary actors. For concreteness, although, we first have a look at an instance wavelet.

We begin by implementing the above equation:

Evaluating code and mathematical formulation, we discover a distinction. The operate itself takes one argument, (t_n); its realization, 4 (`omega`, `Okay`, `t_k`, and `t`). It’s because the `torch` code is vectorized: On the one hand, `omega`, `Okay`, and `t_k`, which, within the system, correspond to (omega_{a}), (Okay), and (t_k) , are scalars. (Within the equation, they’re assumed to be mounted.) `t`, alternatively, is a vector; it can maintain the measurement instances of the collection to be analyzed.

We decide instance values for `omega`, `Okay`, and `t_k`, in addition to a variety of instances to guage the wavelet on, and plot its values:

``````omega <- 6 * pi
Okay <- 6
t_k <- 5

sample_time <- torch_arange(3, 7, 0.0001)

create_wavelet_plot <- operate(omega, Okay, t_k, sample_time) {
morlet <- morlet(omega, Okay, t_k, sample_time)
df <- knowledge.body(
x = as.numeric(sample_time),
actual = as.numeric(morlet\$actual),
imag = as.numeric(morlet\$imag)
) %>%
pivot_longer(-x, names_to = "half", values_to = "worth")
ggplot(df, aes(x = x, y = worth, colour = half)) +
geom_line() +
scale_colour_grey(begin = 0.8, finish = 0.4) +
xlab("time") +
ylab("wavelet worth") +
ggtitle("Morlet wavelet",
subtitle = paste0("ω_a = ", omega / pi, "π , Okay = ", Okay)
) +
theme_minimal()
}

create_wavelet_plot(omega, Okay, t_k, sample_time)``````

What we see here’s a complicated sine curve – notice the actual and imaginary components, separated by a section shift of (pi/2) – that decays on either side of the middle. Wanting again on the equation, we will establish the components chargeable for each options. The primary time period within the equation, (e^{-i omega_{a} (t_n – t_k)}), generates the oscillation; the third, (e^{- omega_a^2 (t_n – t_k )^2 /(2K )^2}), causes the exponential decay away from the middle. (In case you’re questioning concerning the second time period, (e^{-Okay^2}): For given (Okay), it’s only a fixed.)

The third time period truly is a Gaussian, with location parameter (t_k) and scale (Okay). We’ll speak about (Okay) in nice element quickly, however what’s with (t_k)? (t_k) is the middle of the wavelet; for the Morlet wavelet, that is additionally the situation of most amplitude. As distance from the middle will increase, values shortly method zero. That is what is supposed by wavelets being localized: They’re “energetic” solely on a brief vary of time.

## The roles of (Okay) and (omega_a)

Now, we already stated that (Okay) is the dimensions of the Gaussian; it thus determines how far the curve spreads out in time. However there’s additionally (omega_a). Wanting again on the Gaussian time period, it, too, will impression the unfold.

First although, what’s (omega_a)? The subscript (a) stands for “evaluation”; thus, (omega_a) denotes a single frequency being probed.

Now, let’s first examine visually the respective impacts of (omega_a) and (Okay).

``````p1 <- create_wavelet_plot(6 * pi, 4, 5, sample_time)
p2 <- create_wavelet_plot(6 * pi, 6, 5, sample_time)
p3 <- create_wavelet_plot(6 * pi, 8, 5, sample_time)
p4 <- create_wavelet_plot(4 * pi, 6, 5, sample_time)
p5 <- create_wavelet_plot(6 * pi, 6, 5, sample_time)
p6 <- create_wavelet_plot(8 * pi, 6, 5, sample_time)

(p1 | p4) /
(p2 | p5) /
(p3 | p6)``````

Within the left column, we preserve (omega_a) fixed, and differ (Okay). On the best, (omega_a) modifications, and (Okay) stays the identical.

Firstly, we observe that the upper (Okay), the extra the curve will get unfold out. In a wavelet evaluation, which means that extra deadlines will contribute to the remodel’s output, leading to excessive precision as to frequency content material, however lack of decision in time. (We’ll return to this – central – trade-off quickly.)

As to (omega_a), its impression is twofold. On the one hand, within the Gaussian time period, it counteracts – precisely, even – the dimensions parameter, (Okay). On the opposite, it determines the frequency, or equivalently, the interval, of the wave. To see this, check out the best column. Akin to the completely different frequencies, we have now, within the interval between 4 and 6, 4, six, or eight peaks, respectively.

This double function of (omega_a) is the rationale why, all-in-all, it does make a distinction whether or not we shrink (Okay), protecting (omega_a) fixed, or improve (omega_a), holding (Okay) mounted.

This state of issues sounds difficult, however is much less problematic than it might sound. In follow, understanding the function of (Okay) is essential, since we have to decide smart (Okay) values to attempt. As to the (omega_a), alternatively, there will probably be a mess of them, similar to the vary of frequencies we analyze.

So we will perceive the impression of (Okay) in additional element, we have to take a primary have a look at the Wavelet Remodel.

## Wavelet Remodel: An easy implementation

Whereas total, the subject of wavelets is extra multifaceted, and thus, could appear extra enigmatic than Fourier evaluation, the remodel itself is less complicated to understand. It’s a sequence of native convolutions between wavelet and sign. Right here is the system for particular scale parameter (Okay), evaluation frequency (omega_a), and wavelet location (t_k):

[
W_{K, omega_a, t_k} = sum_n x_n Psi_{omega_{a},K,t_{k}}^*(t_n)
]

That is only a dot product, computed between sign and complex-conjugated wavelet. (Right here complicated conjugation flips the wavelet in time, making this convolution, not correlation – a proven fact that issues quite a bit, as you’ll see quickly.)

Correspondingly, easy implementation ends in a sequence of dot merchandise, every similar to a distinct alignment of wavelet and sign. Under, in `wavelet_transform()`, arguments `omega` and `Okay` are scalars, whereas `x`, the sign, is a vector. The result’s the wavelet-transformed sign, for some particular `Okay` and `omega` of curiosity.

``````wavelet_transform <- operate(x, omega, Okay) {
n_samples <- dim(x)
W <- torch_complex(
torch_zeros(n_samples), torch_zeros(n_samples)
)
for (i in 1:n_samples) {
# transfer heart of wavelet
t_k <- x[i, 1]
m <- morlet(omega, Okay, t_k, x[, 1])
# compute native dot product
# notice wavelet is conjugated
dot <- torch_matmul(
m\$conj()\$unsqueeze(1),
x[, 2]\$to(dtype = torch_cfloat())
)
W[i] <- dot
}
W
}``````

To check this, we generate a easy sine wave that has a frequency of 100 Hertz in its first half, and double that within the second.

``````gencos <- operate(amp, freq, section, fs, period) {
x <- torch_arange(0, period, 1 / fs)[1:-2]\$unsqueeze(2)
y <- amp * torch_cos(2 * pi * freq * x + section)
torch_cat(record(x, y), dim = 2)
}

# sampling frequency
fs <- 8000

f1 <- 100
f2 <- 200
section <- 0
period <- 0.25

s1 <- gencos(1, f1, section, fs, period)
s2 <- gencos(1, f2, section, fs, period)

s3 <- torch_cat(record(s1, s2), dim = 1)
s3[(dim(s1) + 1):(dim(s1) * 2), 1] <-
s3[(dim(s1) + 1):(dim(s1) * 2), 1] + period

df <- knowledge.body(
x = as.numeric(s3[, 1]),
y = as.numeric(s3[, 2])
)
ggplot(df, aes(x = x, y = y)) +
geom_line() +
xlab("time") +
ylab("amplitude") +
theme_minimal()`````` An instance sign, consisting of a low-frequency and a high-frequency half.

Now, we run the Wavelet Remodel on this sign, for an evaluation frequency of 100 Hertz, and with a `Okay` parameter of two, discovered by way of fast experimentation:

``````Okay <- 2
omega <- 2 * pi * f1

res <- wavelet_transform(x = s3, omega, Okay)
df <- knowledge.body(
x = as.numeric(s3[, 1]),
y = as.numeric(res\$abs())
)

ggplot(df, aes(x = x, y = y)) +
geom_line() +
xlab("time") +
ylab("Wavelet Remodel") +
theme_minimal()`````` Wavelet Remodel of the above two-part sign. Evaluation frequency is 100 Hertz.

The remodel appropriately picks out the a part of the sign that matches the evaluation frequency. In the event you really feel like, you would possibly wish to double-check what occurs for an evaluation frequency of 200 Hertz.

Now, in actuality we’ll wish to run this evaluation not for a single frequency, however a variety of frequencies we’re all for. And we’ll wish to attempt completely different scales `Okay`. Now, when you executed the code above, you is likely to be anxious that this might take a lot of time.

Effectively, it by necessity takes longer to compute than its Fourier analogue, the spectrogram. For one, that’s as a result of with spectrograms, the evaluation is “simply” two-dimensional, the axes being time and frequency. With wavelets there are, as well as, completely different scales to be explored. And secondly, spectrograms function on entire home windows (with configurable overlap); a wavelet, alternatively, slides over the sign in unit steps.

Nonetheless, the state of affairs isn’t as grave because it sounds. The Wavelet Remodel being a convolution, we will implement it within the Fourier area as an alternative. We’ll try this very quickly, however first, as promised, let’s revisit the subject of various `Okay`.

## Decision in time versus in frequency

We already noticed that the upper `Okay`, the extra spread-out the wavelet. We are able to use our first, maximally easy, instance, to research one fast consequence. What, for instance, occurs for `Okay` set to twenty?

``````Okay <- 20

res <- wavelet_transform(x = s3, omega, Okay)
df <- knowledge.body(
x = as.numeric(s3[, 1]),
y = as.numeric(res\$abs())
)

ggplot(df, aes(x = x, y = y)) +
geom_line() +
xlab("time") +
ylab("Wavelet Remodel") +
theme_minimal()`````` Wavelet Remodel of the above two-part sign, with Okay set to twenty as an alternative of two.

The Wavelet Remodel nonetheless picks out the right area of the sign – however now, as an alternative of a rectangle-like outcome, we get a considerably smoothed model that doesn’t sharply separate the 2 areas.

Notably, the primary 0.05 seconds, too, present appreciable smoothing. The bigger a wavelet, the extra element-wise merchandise will probably be misplaced on the finish and the start. It’s because transforms are computed aligning the wavelet in any respect sign positions, from the very first to the final. Concretely, once we compute the dot product at location `t_k = 1`, only a single pattern of the sign is taken into account.

Other than presumably introducing unreliability on the boundaries, how does wavelet scale have an effect on the evaluation? Effectively, since we’re correlating (convolving, technically; however on this case, the impact, ultimately, is similar) the wavelet with the sign, point-wise similarity is what issues. Concretely, assume the sign is a pure sine wave, the wavelet we’re utilizing is a windowed sinusoid just like the Morlet, and that we’ve discovered an optimum `Okay` that properly captures the sign’s frequency. Then every other `Okay`, be it bigger or smaller, will lead to much less point-wise overlap.

## Performing the Wavelet Remodel within the Fourier area

Quickly, we’ll run the Wavelet Remodel on an extended sign. Thus, it’s time to velocity up computation. We already stated that right here, we profit from time-domain convolution being equal to multiplication within the Fourier area. The general course of then is that this: First, compute the DFT of each sign and wavelet; second, multiply the outcomes; third, inverse-transform again to the time area.

The DFT of the sign is shortly computed:

``F <- torch_fft_fft(s3[ , 2])``

With the Morlet wavelet, we don’t even need to run the FFT: Its Fourier-domain illustration could be said in closed kind. We’ll simply make use of that formulation from the outset. Right here it’s:

``````morlet_fourier <- operate(Okay, omega_a, omega) {
2 * (torch_exp(-torch_square(
Okay * (omega - omega_a) / omega_a
)) -
torch_exp(-torch_square(Okay)) *
torch_exp(-torch_square(Okay * omega / omega_a)))
}``````

Evaluating this assertion of the wavelet to the time-domain one, we see that – as anticipated – as an alternative of parameters `t` and `t_k` it now takes `omega` and `omega_a`. The latter, `omega_a`, is the evaluation frequency, the one we’re probing for, a scalar; the previous, `omega`, the vary of frequencies that seem within the DFT of the sign.

In instantiating the wavelet, there’s one factor we have to pay particular consideration to. In FFT-think, the frequencies are bins; their quantity is decided by the size of the sign (a size that, for its half, instantly is dependent upon sampling frequency). Our wavelet, alternatively, works with frequencies in Hertz (properly, from a consumer’s perspective; since this unit is significant to us). What this implies is that to `morlet_fourier`, as `omega_a` we have to cross not the worth in Hertz, however the corresponding FFT bin. Conversion is completed relating the variety of bins, `dim(x)`, to the sampling frequency of the sign, `fs`:

``````# once more search for 100Hz components
omega <- 2 * pi * f1

# want the bin similar to some frequency in Hz
omega_bin <- f1/fs * dim(s3)``````

We instantiate the wavelet, carry out the Fourier-domain multiplication, and inverse-transform the outcome:

``````Okay <- 3

m <- morlet_fourier(Okay, omega_bin, 1:dim(s3))
prod <- F * m
reworked <- torch_fft_ifft(prod)``````

Placing collectively wavelet instantiation and the steps concerned within the evaluation, we have now the next. (Be aware `wavelet_transform_fourier`, we now, conveniently, cross within the frequency worth in Hertz.)

``````wavelet_transform_fourier <- operate(x, omega_a, Okay, fs) {
N <- dim(x)
omega_bin <- omega_a / fs * N
m <- morlet_fourier(Okay, omega_bin, 1:N)
x_fft <- torch_fft_fft(x)
prod <- x_fft * m
w <- torch_fft_ifft(prod)
w
}``````

We’ve already made vital progress. We’re prepared for the ultimate step: automating evaluation over a variety of frequencies of curiosity. This can lead to a three-dimensional illustration, the wavelet diagram.

## Creating the wavelet diagram

Within the Fourier Remodel, the variety of coefficients we receive is dependent upon sign size, and successfully reduces to half the sampling frequency. With its wavelet analogue, since anyway we’re doing a loop over frequencies, we would as properly resolve which frequencies to research.

Firstly, the vary of frequencies of curiosity could be decided working the DFT. The following query, then, is about granularity. Right here, I’ll be following the advice given in Vistnes’ guide, which is predicated on the relation between present frequency worth and wavelet scale, `Okay`.

Iteration over frequencies is then carried out as a loop:

``````wavelet_grid <- operate(x, Okay, f_start, f_end, fs) {
# downsample evaluation frequency vary
# as per Vistnes, eq. 14.17
num_freqs <- 1 + log(f_end / f_start)/ log(1 + 1/(8 * Okay))
freqs <- seq(f_start, f_end, size.out = flooring(num_freqs))

reworked <- torch_zeros(
num_freqs, dim(x),
dtype = torch_cfloat()
)
for(i in 1:num_freqs) {
w <- wavelet_transform_fourier(x, freqs[i], Okay, fs)
reworked[i, ] <- w
}
record(reworked, freqs)
}``````

Calling `wavelet_grid()` will give us the evaluation frequencies used, along with the respective outputs from the Wavelet Remodel.

Subsequent, we create a utility operate that visualizes the outcome. By default, `plot_wavelet_diagram()` shows the magnitude of the wavelet-transformed collection; it could actually, nevertheless, plot the squared magnitudes, too, in addition to their sq. root, a way a lot beneficial by Vistnes whose effectiveness we’ll quickly have alternative to witness.

The operate deserves a number of additional feedback.

Firstly, similar as we did with the evaluation frequencies, we down-sample the sign itself, avoiding to counsel a decision that’s not truly current. The system, once more, is taken from Vistnes’ guide.

Then, we use interpolation to acquire a brand new time-frequency grid. This step might even be essential if we preserve the unique grid, since when distances between grid factors are very small, R’s `picture()` might refuse to simply accept axes as evenly spaced.

Lastly, notice how frequencies are organized on a log scale. This results in far more helpful visualizations.

``````plot_wavelet_diagram <- operate(x,
freqs,
grid,
Okay,
fs,
f_end,
sort = "magnitude") {
grid <- change(sort,
magnitude = grid\$abs(),
magnitude_squared = torch_square(grid\$abs()),
magnitude_sqrt = torch_sqrt(grid\$abs())
)

# downsample time collection
# as per Vistnes, eq. 14.9
new_x_take_every <- max(Okay / 24 * fs / f_end, 1)
new_x_length <- flooring(dim(grid) / new_x_take_every)
new_x <- torch_arange(
x,
x[dim(x)],
step = x[dim(x)] / new_x_length
)

# interpolate grid
new_grid <- nnf_interpolate(
grid\$view(c(1, 1, dim(grid), dim(grid))),
c(dim(grid), new_x_length)
)\$squeeze()
out <- as.matrix(new_grid)

# plot log frequencies
freqs <- log10(freqs)

picture(
x = as.numeric(new_x),
y = freqs,
z = t(out),
ylab = "log frequency [Hz]",
xlab = "time [s]",
col = hcl.colours(12, palette = "Gentle grays")
)
major <- paste0("Wavelet Remodel, Okay = ", Okay)
sub <- change(sort,
magnitude = "Magnitude",
magnitude_squared = "Magnitude squared",
magnitude_sqrt = "Magnitude (sq. root)"
)

mtext(facet = 3, line = 2, at = 0, adj = 0, cex = 1.3, major)
mtext(facet = 3, line = 1, at = 0, adj = 0, cex = 1, sub)
}``````

Let’s use this on a real-world instance.

## An actual-world instance: Chaffinch’s track

For the case research, I’ve chosen what, to me, was essentially the most spectacular wavelet evaluation proven in Vistnes’ guide. It’s a pattern of a chaffinch’s singing, and it’s accessible on Vistnes’ web site.

``````url <- "http://www.physics.uio.no/pow/wavbirds/chaffinch.wav"

obtain.file(
file.path(url),
destfile = "/tmp/chaffinch.wav"
)``````

We use `torchaudio` to load the file, and convert from stereo to mono utilizing `tuneR`’s appropriately named `mono()`. (For the form of evaluation we’re doing, there isn’t a level in protecting two channels round.)

``````library(torchaudio)
library(tuneR)

wav <- mono(wav, "each")
wav``````
``````Wave Object
Variety of Samples:      1864548
Length (seconds):     42.28
Samplingrate (Hertz):   44100
Channels (Mono/Stereo): Mono
PCM (integer format):   TRUE
Bit (8/16/24/32/64):    16 ``````

For evaluation, we don’t want the entire sequence. Helpfully, Vistnes additionally revealed a advice as to which vary of samples to research.

``````waveform_and_sample_rate <- transform_to_tensor(wav)
x <- waveform_and_sample_rate[]\$squeeze()
fs <- waveform_and_sample_rate[]

# http://www.physics.uio.no/pow/wavbirds/chaffinchInfo.txt
begin <- 34000
N <- 1024 * 128
finish <- begin + N - 1
x <- x[start:end]

dim(x)``````
`` 131072``

How does this look within the time area? (Don’t miss out on the event to truly pay attention to it, in your laptop computer.)

``````df <- knowledge.body(x = 1:dim(x), y = as.numeric(x))
ggplot(df, aes(x = x, y = y)) +
geom_line() +
xlab("pattern") +
ylab("amplitude") +
theme_minimal()``````

Now, we have to decide an inexpensive vary of study frequencies. To that finish, we run the FFT:

On the x-axis, we plot frequencies, not pattern numbers, and for higher visibility, we zoom in a bit.

``````bins <- 1:dim(F)
freqs <- bins / N * fs

# the bin, not the frequency
cutoff <- N/4

df <- knowledge.body(
x = freqs[1:cutoff],
y = as.numeric(F\$abs())[1:cutoff]
)
ggplot(df, aes(x = x, y = y)) +
geom_col() +
xlab("frequency (Hz)") +
ylab("magnitude") +
theme_minimal()``````

Based mostly on this distribution, we will safely prohibit the vary of study frequencies to between, roughly, 1800 and 8500 Hertz. (That is additionally the vary beneficial by Vistnes.)

First, although, let’s anchor expectations by making a spectrogram for this sign. Appropriate values for FFT measurement and window measurement have been discovered experimentally. And although, in spectrograms, you don’t see this achieved typically, I discovered that displaying sq. roots of coefficient magnitudes yielded essentially the most informative output.

``````fft_size <- 1024
window_size <- 1024
energy <- 0.5

spectrogram <- transform_spectrogram(
n_fft = fft_size,
win_length = window_size,
normalized = TRUE,
energy = energy
)

spec <- spectrogram(x)
dim(spec)``````
`` 513 257``

Like we do with wavelet diagrams, we plot frequencies on a log scale.

``````bins <- 1:dim(spec)
freqs <- bins * fs / fft_size
log_freqs <- log10(freqs)

frames <- 1:(dim(spec))
seconds <- (frames / dim(spec))  * (dim(x) / fs)

picture(x = seconds,
y = log_freqs,
z = t(as.matrix(spec)),
ylab = 'log frequency [Hz]',
xlab = 'time [s]',
col = hcl.colours(12, palette = "Gentle grays")
)
major <- paste0("Spectrogram, window measurement = ", window_size)
sub <- "Magnitude (sq. root)"
mtext(facet = 3, line = 2, at = 0, adj = 0, cex = 1.3, major)
mtext(facet = 3, line = 1, at = 0, adj = 0, cex = 1, sub)``````

The spectrogram already exhibits a particular sample. Let’s see what could be achieved with wavelet evaluation. Having experimented with a number of completely different `Okay`, I agree with Vistnes that `Okay = 48` makes for a superb alternative:

``````f_start <- 1800
f_end <- 8500

Okay <- 48
c(grid, freqs) %<-% wavelet_grid(x, Okay, f_start, f_end, fs)
plot_wavelet_diagram(
torch_tensor(1:dim(grid)),
freqs, grid, Okay, fs, f_end,
sort = "magnitude_sqrt"
)``````