# Adjusting Point Sets in Frequency Space Using a Differentiable Fourier Transform

A couple years ago, Bart Wronski wrote an interesting blog post on how the Fourier transform was easily differentiable, which allowed using gradient descent to make blue noise textures. His post is here: https://bartwronski.com/2020/04/26/optimizing-blue-noise-dithering-backpropagation-through-fourier-transform-and-sorting/

In this post I’ll show the details of how to do the same for points sets instead of textures.

There are two web demos that go along with this post. They let you make point sets in 1D and 2D, and move them around in either normal space, or frequency space.
1D Point Set Frequencies: http://demofox.org/PointSetsFT1D.html
2D Point Set Frequencies: http://demofox.org/PointSetsFT2D.html

NOTE: We are going to adjust frequency AMPLITUDES of point sets, but not frequency PHASES.

## 1D Fourier Transform of Point Sets

In a previous blog post I wrote about the discrete Fourier transform of 1D data (https://blog.demofox.org/2016/08/11/understanding-the-discrete-fourier-transform/) which shows the formula:

$X_k = \sum_{n=0}^{N-1} x_n * e^{2 \pi i k n/N}$

As a quick recap, k is the frequency (1 for 1hz, 2 for 2hz, etc), $X_k$ is the complex frequency information for that frequency, n is the index of the current value, $x_n$ is the value at that index, and N is the total number of values.

If you have a 1D point set, one way to Fourier transform it is to make a 1D array of data, fill it with zeros, and then plot those points as ones into that 1D array. You then Fourier transform that array using the formula above.

It’s up to you what size to make that array. If you make the array smaller, the transform will be faster, but it will be less accurate in general. If the array is larger, the transform will be slower, but will be more accurate in general.

Here are some examples of different sizes arrays holding the points set (0.0, 0.25, 0.5, 0.75), where we assume the array holds the values [0,1):

• 4: [1111]
• 8: [10101010]
• 9: [101001010]
• 12: [100100100100]

See how the sized 9 array doesn’t have the points evenly spaced? If you have values that don’t evenly divide the array size, larger arrays make the plots more accurate.

Another way to transform these points though is to think of each of the points as a delta function that is zero everywhere, except at their value, where the function is one. We can then modify the Fourier transform to sum up the contribution of each point towards a specific frequency. Since we know the function of each point is zero everywhere except one location where it is one, we only have to evaluate it at one place. That gives us this:

$X_k = \sum_{p \in P} e^{2 \pi i k p}$

Where again, k is the frequency being analyzed (1 for 1hz, 2 for 2hz, etc), $X_k$ is the complex frequency information for that frequency, P is the total set of points in the point set, and p is the value of the current point in the point set.

## Adjusting 1D Point Sets in Frequency Space

To adjust frequency magnitudes of our point set, we are going to need to show the formula for calculating frequency magnitude and we are going to need to differentiate it to get the gradient of it, so we know which points to move in which direction to increase or decrease specific frequency magnitudes.

To calculate the magnitude of a specific frequency, we start with our equation from the last section.

$X_k = \sum_{p \in P} e^{2 \pi i k p}$

Remembering that $e^{ix} = \cos(x) + i*\sin(x)$, we can turn that to this:

$X_k = \sum_{p \in P} \cos(2 \pi k p) + i * \sin(2 \pi k p)$

The magnitude of the frequency is the length of the vector (real, imaginary), so let’s break up the real and imaginary parts:

$\Re(X_k) = \sum_{p \in P} \cos(2 \pi k p)$

$\Im(X_k) = \sum_{p \in P} \sin(2 \pi k p)$

The magnitude would then be:

$\lVert X_k \rVert= \sqrt{(\sum_{p \in P} \cos(2 \pi k p))^2 + (\sum_{p \in P} \sin(2 \pi k p))^2}$

We’ll remove the square root to make it easier to differentiate, and get this:

$\lVert X_k \rVert^2= (\sum_{p \in P} \cos(2 \pi k p))^2 + (\sum_{p \in P} \sin(2 \pi k p))^2$

Now we need to get the derivative of that function for each point p, which together is called the gradient. This tells us how far to move each point to make the frequency magnitude increase by 1. It also tells us the direction since a negative number would mean to move the point to the left and a positive number would mean to move the point to the right.

We can use the sum rule to break our function into 2 simpler derivatives and deal with them separately:

$(f(x)+g(x))' = f'(x) + g'(x)$

$(\lVert X_k \rVert^2)'= ((\sum_{p \in P} \cos(2 \pi k p))^2)' + ((\sum_{p \in P} \sin(2 \pi k p))^2)'$

Let’s start by differentiating the real (cosine) term:

$f' = ((\sum_{p \in P} \cos(2 \pi k p))^2)'$

We can use the product rule to deal with the squaring:

$(f(x)^2)' = (f(x)*f(x))' = f'(x)f(x) + f(x)f'(x) = 2f(x)f'(x)$

$f' = 2*(\sum_{p \in P} \cos(2 \pi k p))*(\sum_{p \in P} \cos(2 \pi k p))'$

That is all differentiated except for the last term:

$g' = (\sum_{p \in P} \cos(2 \pi k p))'$

If we are differentiating this function g for a specific p (we’ll call it $p_i$), all of the terms are constants and disappear, except for the $p_i$ term. This is easily differentiable because $\cos(ax)' = -a \sin(ax)$

$g' = (\cos(2 \pi k p_i))'$

$g' = -2 \pi k \sin(2 \pi k p_i)$

We can plug our g result back into f and get this as the derivative of our cosine (real) term:

$f' = 2*(\sum_{p \in P} \cos(2 \pi k p))*-2 \pi k \sin(2 \pi k p_i)$

$f' = -4 \pi k*(\sum_{p \in P} \cos(2 \pi k p))*\sin(2 \pi k p_i)$

If we do the same steps for the sine (imaginary) term we end up with this:

$h' = 4 \pi k*(\sum_{p \in P} \sin(2 \pi k p))*\cos(2 \pi k p_i)$

We can then add f’ and h’ together to get the full equation we were looking for!

$(\lVert X_k \rVert^2)' = -4 \pi k*(\sum_{p \in P} \cos(2 \pi k p))*\sin(2 \pi k p_i) + 4 \pi k*(\sum_{p \in P} \sin(2 \pi k p))*\cos(2 \pi k p_i)$

Here is that monster as javascript code from the demo:

If you calculate this value for each point and put it into an array, you’ll have the gradient of the frequency k’s magnitude (squared). If you want to increase the frequency k’s magnitude, you add the gradient to the points. If you want to decrease the frequency k’s magnitude, you subtract the gradient from the points.

The gradient is only guaranteed to be accurate for an infinitesimally small step, so what I do in the web app is normalize this gradient vector, and multiply it by a step size before I move the points. I also multiply the step size by how many points there are so that the distance that each point moves doesn’t decrease as more points are added (otherwise they would have to share a unit length movement among N points, which gets shorter as N gets bigger). I also have a step count which allows this operation to be done M times, re-evaluating the gradient each time. This results in an interactive user controlled gradient descent.

You can play with the demo here to see it all in action:

1D Point Set Frequencies: http://demofox.org/PointSetsFT1D.html

For more information about why you might want to normalize a gradient when doing gradient descent, give this a read: https://jermwatt.github.io/machine_learning_refined/notes/3_First_order_methods/3_9_Normalized.html

## 2D Fourier Transform of Point Sets

The 2D Fourier transform is just the 1D Fourier transform done on each axis. Because of this, you have a horizontal frequency j and a vertical frequency k, for an image that is MxN pixels:

$X_{j,k} = \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} x_{m,n} * e^{2 \pi i j m/M}* e^{2 \pi i k n/N}$

$X_{j,k}$ is the complex frequency information for the frequency, j is the horizontal frequency, k is the vertical frequency, $x_{m,n}$ is the pixel value at location (m, n).

We can change this to be a Fourier transform of points p again like this:

$X_{j,k} = \sum_{p \in P} e^{2 \pi i j p_x}* e^{2 \pi i k p_y}$

Where $p_x$ is the x component of the point and $p_y$ is the y component of the point.

We can use the identity $e^{a} * e^{b} = e^{a+b}$ to simplify the equation a bit:

$X_{j,k} = \sum_{p \in P} e^{2 \pi i j p_x + 2 \pi i k p_y}$

$X_{j,k} = \sum_{p \in P} e^{2 \pi i (j p_x + k p_y)}$

… and we are done 🙂

## Adjusting 2D Point Sets in Frequency Space

Let’s jump to breaking the function into real and imaginary parts.

$\Re(X_{j,k}) = \sum_{p \in P} \cos(2 \pi (j p_x + k p_y))$

$\Im(X_{j,k}) = \sum_{p \in P} \sin(2 \pi (j p_x + k p_y))$

The magnitude would then be:

$\lVert X_k \rVert= \sqrt{(\sum_{p \in P} \cos(2 \pi (j p_x + k p_y)))^2 + (\sum_{p \in P} \sin(2 \pi (j p_x + k p_y)))^2}$

We can square both sides again to make it easier to differentiate:

$\lVert X_k \rVert ^2= (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y)))^2 + (\sum_{p \in P} \sin(2 \pi (j p_x + k p_y)))^2$

Remembering that we can differentiate each term independently, let’s start with the cosine (real) part again and start differentiation using the product rule:

$f' = 2 * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y))) * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y)))'$

We now have two variables to get derivatives for though, the x and y components of a specific point. We’ll call them $p_{ix}$ and $p_{iy}$. We’ll start with $p_{ix}$. Remember that all terms of the sum except the ones with $p_i$ in them are constants and will become zero.

$\frac{df}{dp_{ix}} = 2 * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y))) * -2 \pi j \sin(2 \pi (j p_{ix} + k p_{iy}))$

$\frac{df}{dp_{ix}} = -4 \pi j * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y))) * \sin(2 \pi (j p_{ix} + k p_{iy}))$

We can do the say with the y component and get:

$\frac{df}{dp_{iy}} = -4 \pi k * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y))) * \sin(2 \pi (j p_{ix} + k p_{iy}))$

We can do the same process with the sine (imaginary) part and get these two equations:

$\frac{dh}{dp_{ix}} = 4 \pi j * (\sum_{p \in P} \sin(2 \pi (j p_x + k p_y))) * \cos(2 \pi (j p_{ix} + k p_{iy}))$

$\frac{dh}{dp_{iy}} = 4 \pi k * (\sum_{p \in P} \sin(2 \pi (j p_x + k p_y))) * \cos(2 \pi (j p_{ix} + k p_{iy}))$

We then add the real and imaginary x functions together to get a value for x, and we add the real and imaginary y functions together to get a value for y.

$\lVert X_k \rVert ^2 \frac{d}{dp_{ix}} = -4 \pi j * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y))) * \sin(2 \pi (j p_{ix} + k p_{iy})) + 4 \pi j * (\sum_{p \in P} \sin(2 \pi (j p_x + k p_y))) * \cos(2 \pi (j p_{ix} + k p_{iy}))$

$\lVert X_k \rVert ^2 \frac{d}{dp_{iy}} = -4 \pi k * (\sum_{p \in P} \cos(2 \pi (j p_x + k p_y))) * \sin(2 \pi (j p_{ix} + k p_{iy})) + 4 \pi k * (\sum_{p \in P} \sin(2 \pi (j p_x + k p_y))) * \cos(2 \pi (j p_{ix} + k p_{iy}))$

That sure is a mouthful, but it ends up looking a lot nicer in code. Javascript version from the demo below:

So the gradient of the function now has a 2D vector for every entry instead of a scalar. I’m not sure if this is still called a gradient. To normalize it, I sum up the length of each of the 2D vectors and then divide them all by that value. I then multiply by the number of points for the same reason as before, and I also have a step size multiplier, and a step count, just as in the 1D case.

You can play with the demo here:

2D Point Set Frequencies: http://demofox.org/PointSetsFT2D.html

## Challenge

If you are up for a challenge, I have one for you!

Blue noise points are “randomized but roughly evenly spaced” and yet, I once stumbled on how to make a point set which shows a blue noise spectrum but has clumps. It’s a preset in the 1D demo and it looks like this:

Red noise on the other hand is defined as randomized but clumping points.

The challenge is this: Can you find points which are reasonably spread out on the numberline, or otherwise not clumping, but show a red noise spectrum?

I’d be real interested in seeing it if so. You can leave a comment here or find me on twitter at https://twitter.com/Atrix256

Thanks for reading and hope you found this useful and/or interesting!