RStudio AI Blog: Discrete Fourier Transform

Note: This put up is an excerpt from the forthcoming e book, Deep Learning and Scientific Computing with R torch. The chapter in query is on the Discrete Fourier Transform (DFT), and is situated partially three. Part three is devoted to scientific computation past deep studying.
There are two chapters on the Fourier Transform. The first strives to, in as “verbal” and lucid a method as was attainable to me, solid a lightweight on what’s behind the magic; it additionally exhibits how, surprisingly, you possibly can code the DFT in merely half a dozen strains. The second focuses on quick implementation (the Fast Fourier Transform, or FFT), once more with each conceptual/explanatory in addition to sensible, code-it-yourself components.
Together, these cowl much more materials than might sensibly match right into a weblog put up; subsequently, please think about what follows extra as a “teaser” than a completely fledged article.

In the sciences, the Fourier Transform is nearly in all places. Stated very typically, it converts knowledge from one illustration to a different, with none lack of info (if performed appropriately, that’s.) If you employ torch, it’s only a perform name away: torch_fft_fft() goes a method, torch_fft_ifft() the opposite. For the consumer, that’s handy – you “just” have to know interpret the outcomes. Here, I wish to assist with that. We begin with an instance perform name, taking part in round with its output, after which, attempt to get a grip on what’s going on behind the scenes.

Understanding the output of torch_fft_fft()

As we care about precise understanding, we begin from the best attainable instance sign, a pure cosine that performs one revolution over the whole sampling interval.

Starting level: A cosine of frequency 1

The method we set issues up, there might be sixty-four samples; the sampling interval thus equals N = 64. The content material of frequency(), the beneath helper perform used to assemble the sign, displays how we characterize the cosine. Namely:

[
f(x) = cos(frac{2 pi}{N} k x)
]

Here (x) values progress over time (or house), and (okay) is the frequency index. A cosine is periodic with interval (2 pi); so if we wish it to first return to its beginning state after sixty-four samples, and (x) runs between zero and sixty-three, we’ll need (okay) to be equal to (1). Like that, we’ll attain the preliminary state once more at place (x = frac{2 pi}{64} * 1 * 64).

Let’s rapidly verify this did what it was alleged to:

df <- data.frame(x = sample_positions, y = as.numeric(x))

ggplot(df, aes(x = x, y = y)) +
  geom_line() +
  xlab("time") +
  ylab("amplitude") +
  theme_minimal()

Now that we have now the enter sign, torch_fft_fft() computes for us the Fourier coefficients, that’s, the significance of the varied frequencies current within the sign. The variety of frequencies thought of will equal the variety of sampling factors: So (X) might be of size sixty-four as nicely.

(In our instance, you’ll discover that the second half of coefficients will equal the primary in magnitude. This is the case for each real-valued sign. In such circumstances, you would name torch_fft_rfft() as an alternative, which yields “nicer” (within the sense of shorter) vectors to work with. Here although, I wish to clarify the overall case, since that’s what you’ll discover performed in most expositions on the subject.)

Even with the sign being actual, the Fourier coefficients are complicated numbers. There are 4 methods to examine them. The first is to extract the actual half:

[1]  0 32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
[57] 0 0 0 0 0 0 0 32

Only a single coefficient is non-zero, the one at place 1. (We begin counting from zero, and will discard the second half, as defined above.)

Now wanting on the imaginary half, we discover it’s zero all through:

[1]  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[57] 0 0 0 0 0 0 0 0

At this level we all know that there’s only a single frequency current within the sign, specifically, that at (okay = 1). This matches (and it higher needed to) the best way we constructed the sign: specifically, as conducting a single revolution over the whole sampling interval.

Since, in concept, each coefficient might have non-zero actual and imaginary components, typically what you’d report is the magnitude (the sq. root of the sum of squared actual and imaginary components):

[1]  0 32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
[57] 0 0 0 0 0 0 0 32

Unsurprisingly, these values precisely replicate the respective actual components.

Finally, there’s the part, indicating a attainable shift of the sign (a pure cosine is unshifted). In torch, we have now torch_angle() complementing torch_abs(), however we have to keep in mind roundoff error right here. We know that in every however a single case, the actual and imaginary components are each precisely zero; however as a result of finite precision in how numbers are offered in a pc, the precise values will typically not be zero. Instead, they’ll be very small. If we take one in every of these “fake non-zeroes” and divide it by one other, as occurs within the angle calculation, massive values may end up. To stop this from taking place, our customized implementation rounds each inputs earlier than triggering the division.

part <- perform(Ft, threshold = 1e5) {
  torch_atan2(
    torch_abs(torch_round(Ft$imag * threshold)),
    torch_abs(torch_round(Ft$actual * threshold))
  )
}

as.numeric(part(Ft)) %>% spherical(5)

[1]  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[57] 0 0 0 0 0 0 0 0

As anticipated, there isn’t a part shift within the sign.

Let’s visualize what we discovered.

create_plot <- perform(x, y, amount) {
  df <- data.frame(
    x_ = x,
    y_ = as.numeric(y) %>% spherical(5)
  )
  ggplot(df, aes(x = x_, y = y_)) +
    geom_col() +
    xlab("frequency") +
    ylab(amount) +
    theme_minimal()
}

p_real <- create_plot(
  sample_positions,
  real_part,
  "actual half"
)
p_imag <- create_plot(
  sample_positions,
  imag_part,
  "imaginary half"
)
p_magnitude <- create_plot(
  sample_positions,
  magnitude,
  "magnitude"
)
p_phase <- create_plot(
  sample_positions,
  part(Ft),
  "part"
)

p_real + p_imag + p_magnitude + p_phase

It’s honest to say that we have now no cause to doubt what torch_fft_fft() has performed. But with a pure sinusoid like this, we are able to perceive precisely what’s occurring by computing the DFT ourselves, by hand. Doing this now will considerably assist us later, after we’re writing the code.

Reconstructing the magic

One caveat about this part. With a subject as wealthy because the Fourier Transform, and an viewers who I think about to range broadly on a dimension of math and sciences schooling, my possibilities to fulfill your expectations, expensive reader, have to be very near zero. Still, I wish to take the chance. If you’re an knowledgeable on these items, you’ll anyway be simply scanning the textual content, searching for items of torch code. If you’re reasonably aware of the DFT, you should still like being reminded of its interior workings. And – most significantly – for those who’re fairly new, and even fully new, to this subject, you’ll hopefully take away (a minimum of) one factor: that what looks like one of many biggest wonders of the universe (assuming there’s a actuality someway similar to what goes on in our minds) could be a marvel, however neither “magic” nor a factor reserved to the initiated.

In a nutshell, the Fourier Transform is a foundation transformation. In the case of the DFT – the Discrete Fourier Transform, the place time and frequency representations each are finite vectors, not features – the brand new foundation seems like this:

[
begin{aligned}
&mathbf{w}^{0n}_N = e^{ifrac{2 pi}{N}* 0 * n} = 1
&mathbf{w}^{1n}_N = e^{ifrac{2 pi}{N}* 1 * n} = e^{ifrac{2 pi}{N} n}
&mathbf{w}^{2n}_N = e^{ifrac{2 pi}{N}* 2 * n} = e^{ifrac{2 pi}{N}2n}& …
&mathbf{w}^{(N-1)n}_N = e^{ifrac{2 pi}{N}* (N-1) * n} = e^{ifrac{2 pi}{N}(N-1)n}
end{aligned}
]

Here (N), as earlier than, is the variety of samples (64, in our case); thus, there are (N) foundation vectors. With (okay) working via the idea vectors, they are often written:

[
mathbf{w}^{kn}_N = e^{ifrac{2 pi}{N}k n}
] {#eq-dft-1}

Like (okay), (n) runs from (0) to (N-1). To perceive what these foundation vectors are doing, it’s useful to briefly change to a shorter sampling interval, (N = 4), say. If we accomplish that, we have now 4 foundation vectors: (mathbf{w}^{0n}_N), (mathbf{w}^{1n}_N), (mathbf{w}^{2n}_N), and (mathbf{w}^{3n}_N). The first one seems like this:

[
mathbf{w}^{0n}_N
=
begin{bmatrix}
e^{ifrac{2 pi}{4}* 0 * 0}
e^{ifrac{2 pi}{4}* 0 * 1}
e^{ifrac{2 pi}{4}* 0 * 2}
e^{ifrac{2 pi}{4}* 0 * 3}
end{bmatrix}
=
begin{bmatrix}
1
1
1
1
end{bmatrix}
]

The second, like so:

[
mathbf{w}^{1n}_N
=
begin{bmatrix}
e^{ifrac{2 pi}{4}* 1 * 0}
e^{ifrac{2 pi}{4}* 1 * 1}
e^{ifrac{2 pi}{4}* 1 * 2}
e^{ifrac{2 pi}{4}* 1 * 3}
end{bmatrix}
=
begin{bmatrix}
1
e^{ifrac{pi}{2}}
e^{i pi}
e^{ifrac{3 pi}{4}}
end{bmatrix}
=
begin{bmatrix}
1
i
-1
-i
end{bmatrix}
]

This is the third:

[
mathbf{w}^{2n}_N
=
begin{bmatrix}
e^{ifrac{2 pi}{4}* 2 * 0}
e^{ifrac{2 pi}{4}* 2 * 1}
e^{ifrac{2 pi}{4}* 2 * 2}
e^{ifrac{2 pi}{4}* 2 * 3}
end{bmatrix}
=
begin{bmatrix}
1
e^{ipi}
e^{i 2 pi}
e^{ifrac{3 pi}{2}}
end{bmatrix}
=
begin{bmatrix}
1
-1
1
-1
end{bmatrix}
]

And lastly, the fourth:

[
mathbf{w}^{3n}_N
=
begin{bmatrix}
e^{ifrac{2 pi}{4}* 3 * 0}
e^{ifrac{2 pi}{4}* 3 * 1}
e^{ifrac{2 pi}{4}* 3 * 2}
e^{ifrac{2 pi}{4}* 3 * 3}
end{bmatrix}
=
begin{bmatrix}
1
e^{ifrac{3 pi}{2}}
e^{i 3 pi}
e^{ifrac{9 pi}{2}}
end{bmatrix}
=
begin{bmatrix}
1
-i
-1
i
end{bmatrix}
]

We can characterize these 4 foundation vectors by way of their “speed”: how briskly they transfer across the unit circle. To do that, we merely have a look at the rightmost column vectors, the place the ultimate calculation outcomes seem. The values in that column correspond to positions pointed to by the revolving foundation vector at completely different time limits. This implies that taking a look at a single “update of position”, we are able to see how briskly the vector is transferring in a single time step.

Looking first at (mathbf{w}^{0n}_N), we see that it doesn’t transfer in any respect. (mathbf{w}^{1n}_N) goes from (1) to (i) to (-1) to (-i); yet one more step, and it might be again the place it began. That’s one revolution in 4 steps, or a step measurement of (frac{pi}{2}). Then (mathbf{w}^{2n}_N) goes at double that tempo, transferring a distance of (pi) alongside the circle. That method, it finally ends up finishing two revolutions total. Finally, (mathbf{w}^{3n}_N) achieves three full loops, for a step measurement of (frac{3 pi}{2}).

The factor that makes these foundation vectors so helpful is that they’re mutually orthogonal. That is, their dot product is zero:

[
langle mathbf{w}^{kn}_N, mathbf{w}^{ln}_N rangle = sum_{n=0}^{N-1} ({e^{ifrac{2 pi}{N}k n}})^* e^{ifrac{2 pi}{N}l n} = sum_{n=0}^{N-1} ({e^{-ifrac{2 pi}{N}k n}})e^{ifrac{2 pi}{N}l n} = 0
] {#eq-dft-2}

Let’s take, for instance, (mathbf{w}^{2n}_N) and (mathbf{w}^{3n}_N). Indeed, their dot product evaluates to zero.

[
begin{bmatrix}
1 & -1 & 1 & -1
end{bmatrix}
begin{bmatrix}
1
-i
-1
i
end{bmatrix}
=
1 + i + (-1) + (-i) = 0
]

Now, we’re about to see how the orthogonality of the Fourier foundation considerably simplifies the calculation of the DFT. Did you discover the similarity between these foundation vectors and the best way we wrote the instance sign? Here it’s once more:

[
f(x) = cos(frac{2 pi}{N} k x)
]

If we handle to characterize this perform by way of the idea vectors (mathbf{w}^{kn}_N = e^{ifrac{2 pi}{N}okay n}), the interior product between the perform and every foundation vector might be both zero (the “default”) or a a number of of 1 (in case the perform has a part matching the idea vector in query). Luckily, sines and cosines can simply be transformed into complicated exponentials. In our instance, that is how that goes:

[
begin{aligned}
mathbf{x}_n &= cos(frac{2 pi}{64} n)
&= frac{1}{2} (e^{ifrac{2 pi}{64} n} + e^{-ifrac{2 pi}{64} n})
&= frac{1}{2} (e^{ifrac{2 pi}{64} n} + e^{ifrac{2 pi}{64} 63n})
&= frac{1}{2} (mathbf{w}^{1n}_N + mathbf{w}^{63n}_N)
end{aligned}
]

Here step one straight outcomes from Euler’s components, and the second displays the truth that the Fourier coefficients are periodic, with frequency -1 being the identical as 63, -2 equaling 62, and so forth.

Now, the (okay)th Fourier coefficient is obtained by projecting the sign onto foundation vector (okay).

Due to the orthogonality of the idea vectors, solely two coefficients won’t be zero: these for (mathbf{w}^{1n}_N) and (mathbf{w}^{63n}_N). They are obtained by computing the interior product between the perform and the idea vector in query, that’s, by summing over (n). For every (n) ranging between (0) and (N-1), we have now a contribution of (frac{1}{2}), leaving us with a remaining sum of (32) for each coefficients. For instance, for (mathbf{w}^{1n}_N):

[
begin{aligned}
X_1 &= langle mathbf{w}^{1n}_N, mathbf{x}_n rangle
&= langle mathbf{w}^{1n}_N, frac{1}{2} (mathbf{w}^{1n}_N + mathbf{w}^{63n}_N) rangle
&= frac{1}{2} * 64
&= 32
end{aligned}
]

And analogously for (X_{63}).

Now, wanting again at what torch_fft_fft() gave us, we see we have been in a position to arrive on the similar end result. And we’ve discovered one thing alongside the best way.

As lengthy as we stick with alerts composed of a number of foundation vectors, we are able to compute the DFT on this method. At the tip of the chapter, we’ll develop code that may work for all alerts, however first, let’s see if we are able to dive even deeper into the workings of the DFT. Three issues we’ll wish to discover:

• What would occur if frequencies modified – say, a melody have been sung at a better pitch?

• What about amplitude adjustments – say, the music have been performed twice as loud?

• What about part – e.g., there have been an offset earlier than the piece began?

In all circumstances, we’ll name torch_fft_fft() solely as soon as we’ve decided the end result ourselves.

And lastly, we’ll see how complicated sinusoids, made up of various elements, can nonetheless be analyzed on this method, supplied they are often expressed by way of the frequencies that make up the idea.

Varying frequency

Assume we quadrupled the frequency, giving us a sign that appeared like this:

[
mathbf{x}_n = cos(frac{2 pi}{N}*4*n)
]

Following the identical logic as above, we are able to specific it like so:

[
mathbf{x}_n = frac{1}{2} (mathbf{w}^{4n}_N + mathbf{w}^{60n}_N)
]

We already see that non-zero coefficients might be obtained just for frequency indices (4) and (60). Picking the previous, we get hold of

[
begin{aligned}
X_4 &= langle mathbf{w}^{4n}_N, mathbf{x}_n rangle
&= langle mathbf{w}^{4n}_N, frac{1}{2} (mathbf{w}^{4n}_N + mathbf{w}^{60n}_N) rangle
&= 32
end{aligned}
]

For the latter, we’d arrive on the similar end result.

Now, let’s make certain our evaluation is right. The following code snippet incorporates nothing new; it generates the sign, calculates the DFT, and plots them each.

x <- torch_cos(frequency(4, N) * sample_positions)

plot_ft <- perform(x)  plot_spacer()) /
    (p_real 

plot_ft(x)

This does certainly verify our calculations.

A particular case arises when sign frequency rises to the best one “allowed”, within the sense of being detectable with out aliasing. That would be the case at one half of the variety of sampling factors. Then, the sign will seem like so:

[
mathbf{x}_n = frac{1}{2} (mathbf{w}^{32n}_N + mathbf{w}^{32n}_N)
]

Consequently, we find yourself with a single coefficient, similar to a frequency of 32 revolutions per pattern interval, of double the magnitude (64, thus). Here are the sign and its DFT:

x <- torch_cos(frequency(32, N) * sample_positions)
plot_ft(x)

Varying amplitude

Now, let’s take into consideration what occurs after we range amplitude. For instance, say the sign will get twice as loud. Now, there might be a multiplier of two that may be taken outdoors the interior product. In consequence, the one factor that adjustments is the magnitude of the coefficients.

Let’s confirm this. The modification is predicated on the instance we had earlier than the final one, with 4 revolutions over the sampling interval:

x <- 2 * torch_cos(frequency(4, N) * sample_positions)
plot_ft(x)

So far, we have now not as soon as seen a coefficient with non-zero imaginary half. To change this, we add in part.

Adding part

Changing the part of a sign means shifting it in time. Our instance sign is a cosine, a perform whose worth is 1 at (t=0). (That additionally was the – arbitrarily chosen – place to begin of the sign.)

Now assume we shift the sign ahead by (frac{pi}{2}). Then the height we have been seeing at zero strikes over to (frac{pi}{2}); and if we nonetheless begin “recording” at zero, we should discover a worth of zero there. An equation describing that is the next. For comfort, we assume a sampling interval of (2 pi) and (okay=1), in order that the instance is a straightforward cosine:

[
f(x) = cos(x – phi)
]

The minus signal might look unintuitive at first. But it does make sense: We now wish to get hold of a price of 1 at (x=frac{pi}{2}), so (x – phi) ought to consider to zero. (Or to any a number of of (pi).) Summing up, a delay in time will seem as a unfavorable part shift.

Now, we’re going to calculate the DFT for a shifted model of our instance sign. But for those who like, take a peek on the phase-shifted model of the time-domain image now already. You’ll see {that a} cosine, delayed by (frac{pi}{2}), is nothing else than a sine beginning at 0.

To compute the DFT, we comply with our familiar-by-now technique. The sign now seems like this:

[
mathbf{x}_n = cos(frac{2 pi}{N}*4*x – frac{pi}{2})
]

First, we specific it by way of foundation vectors:

[
begin{aligned}
mathbf{x}_n &= cos(frac{2 pi}{64} 4 n – frac{pi}{2})
&= frac{1}{2} (e^{ifrac{2 pi}{64} 4n – frac{pi}{2}} + e^{ifrac{2 pi}{64} 60n – frac{pi}{2}})
&= frac{1}{2} (e^{ifrac{2 pi}{64} 4n} e^{-i frac{pi}{2}} + e^{ifrac{2 pi}{64} 60n} e^{ifrac{pi}{2}})
&= frac{1}{2} (e^{-i frac{pi}{2}} mathbf{w}^{4n}_N + e^{i frac{pi}{2}} mathbf{w}^{60n}_N)
end{aligned}
]

Again, we have now non-zero coefficients just for frequencies (4) and (60). But they’re complicated now, and each coefficients are now not an identical. Instead, one is the complicated conjugate of the opposite. First, (X_4):

[
begin{aligned}
X_4 &= langle mathbf{w}^{4n}_N, mathbf{x}_n rangle
&=langle mathbf{w}^{4n}_N, frac{1}{2} (e^{-i frac{pi}{2}} mathbf{w}^{4n}_N + e^{i frac{pi}{2}} mathbf{w}^{60n}_N) rangle
&= 32 *e^{-i frac{pi}{2}}
&= -32i
end{aligned}
]

And right here, (X_{60}):

[
begin{aligned}
X_{60} &= langle mathbf{w}^{60n}_N, mathbf{x}_N rangle
&= 32 *e^{i frac{pi}{2}}
&= 32i
end{aligned}
]

As standard, we examine our calculation utilizing torch_fft_fft().

x <- torch_cos(frequency(4, N) * sample_positions - pi / 2)

plot_ft(x)

For a pure sine wave, the non-zero Fourier coefficients are imaginary. The part shift within the coefficients, reported as (frac{pi}{2}), displays the time delay we utilized to the sign.

Finally – earlier than we write some code – let’s put all of it collectively, and have a look at a wave that has greater than a single sinusoidal part.

Superposition of sinusoids

The sign we assemble should be expressed by way of the idea vectors, however it’s now not a pure sinusoid. Instead, it’s a linear mixture of such:

[
begin{aligned}
mathbf{x}_n &= 3 sin(frac{2 pi}{64} 4n) + 6 cos(frac{2 pi}{64} 2n) +2cos(frac{2 pi}{64} 8n)
end{aligned}
]

I gained’t undergo the calculation intimately, however it’s no completely different from the earlier ones. You compute the DFT for every of the three elements, and assemble the outcomes. Without any calculation, nevertheless, there’s fairly a couple of issues we are able to say:

• Since the sign consists of two pure cosines and one pure sine, there might be 4 coefficients with non-zero actual components, and two with non-zero imaginary components. The latter might be complicated conjugates of one another.
• From the best way the sign is written, it’s straightforward to find the respective frequencies, as nicely: The all-real coefficients will correspond to frequency indices 2, 8, 56, and 62; the all-imaginary ones to indices 4 and 60.
• Finally, amplitudes will end result from multiplying with (frac{64}{2}) the scaling components obtained for the person sinusoids.

Let’s examine:

x <- 3 * torch_sin(frequency(4, N) * sample_positions) +
  6 * torch_cos(frequency(2, N) * sample_positions) +
  2 * torch_cos(frequency(8, N) * sample_positions)

plot_ft(x)

Now, how can we calculate the DFT for much less handy alerts?

Coding the DFT

Fortunately, we already know what needs to be performed. We wish to mission the sign onto every of the idea vectors. In different phrases, we’ll be computing a bunch of interior merchandise. Logic-wise, nothing adjustments: The solely distinction is that basically, it won’t be attainable to characterize the sign by way of only a few foundation vectors, like we did earlier than. Thus, all projections will really need to be calculated. But isn’t automation of tedious duties one factor we have now computer systems for?

Let’s begin by stating enter, output, and central logic of the algorithm to be carried out. As all through this chapter, we keep in a single dimension. The enter, thus, is a one-dimensional tensor, encoding a sign. The output is a one-dimensional vector of Fourier coefficients, of the identical size because the enter, every holding details about a frequency. The central concept is: To get hold of a coefficient, mission the sign onto the corresponding foundation vector.

To implement that concept, we have to create the idea vectors, and for each, compute its interior product with the sign. This might be performed in a loop. Surprisingly little code is required to perform the objective:

dft <- perform(x) {
  n_samples <- size(x)

  n <- torch_arange(0, n_samples - 1)$unsqueeze(1)

  Ft <- torch_complex(
    torch_zeros(n_samples), torch_zeros(n_samples)
  )

  for (okay in 0:(n_samples - 1)) {
    w_k <- torch_exp(-1i * 2 * pi / n_samples * okay * n)
    dot <- torch_matmul(w_k, x$to(dtype = torch_cfloat()))
    Ft[k + 1] <- dot
  }
  Ft
}

To take a look at the implementation, we are able to take the final sign we analysed, and evaluate with the output of torch_fft_fft().

[1]  0 0 192 0 0 0 0 0 64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
[57] 64 0 0 0 0 0 192 0

[1]  0 0 0 0 -96 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
[57] 0 0 0 0 96 0 0 0

Reassuringly – for those who look again – the outcomes are the identical.

Above, did I say “little code”? In reality, a loop just isn’t even wanted. Instead of working with the idea vectors one-by-one, we are able to stack them in a matrix. Then every row will maintain the conjugate of a foundation vector, and there might be (N) of them. The columns correspond to positions (0) to (N-1); there might be (N) of them as nicely. For instance, that is how the matrix would search for (N=4):

[
mathbf{W}_4
=
begin{bmatrix}
e^{-ifrac{2 pi}{4}* 0 * 0} & e^{-ifrac{2 pi}{4}* 0 * 1} & e^{-ifrac{2 pi}{4}* 0 * 2} & e^{-ifrac{2 pi}{4}* 0 * 3}
e^{-ifrac{2 pi}{4}* 1 * 0} & e^{-ifrac{2 pi}{4}* 1 * 1} & e^{-ifrac{2 pi}{4}* 1 * 2} & e^{-ifrac{2 pi}{4}* 1 * 3}
e^{-ifrac{2 pi}{4}* 2 * 0} & e^{-ifrac{2 pi}{4}* 2 * 1} & e^{-ifrac{2 pi}{4}* 2 * 2} & e^{-ifrac{2 pi}{4}* 2 * 3}
e^{-ifrac{2 pi}{4}* 3 * 0} & e^{-ifrac{2 pi}{4}* 3 * 1} & e^{-ifrac{2 pi}{4}* 3 * 2} & e^{-ifrac{2 pi}{4}* 3 * 3}
end{bmatrix}
] {#eq-dft-3}

Or, evaluating the expressions:

[
mathbf{W}_4
=
begin{bmatrix}
1 & 1 & 1 & 1
1 & -i & -1 & i
1 & -1 & 1 & -1
1 & i & -1 & -i
end{bmatrix}
]

With that modification, the code seems much more elegant:

dft_vec <- perform(x) {
  n_samples <- size(x)

  n <- torch_arange(0, n_samples - 1)$unsqueeze(1)
  okay <- torch_arange(0, n_samples - 1)$unsqueeze(2)

  mat_k_m <- torch_exp(-1i * 2 * pi / n_samples * okay * n)

  torch_matmul(mat_k_m, x$to(dtype = torch_cfloat()))
}

As you possibly can simply confirm, the end result is identical.

Thanks for studying!

Photo by Trac Vu on Unsplash