Yannik Nelson - Harmonic Oscillators

Overview

For a project I'm working on, I needed to use harmonic oscillators (or rather a series of harmonic oscillators added together, aka a super-position of harmonic oscillators) to make the motion of some actors feel more organic.

I looked into Unreal's built in wave sampling, however I discovered the exisiting classes didn't provide the interface/functionality I was looking for. As such, I built my own, as well as a way to visualise the resulting wave form of the superposition.

Harmonic Oscillator Desired Interface

My desired interface for the super-position was to have functions to add/remove sine waves, as well as a function to sample the super-position at a given timestamp.
Similarly my desired interface for sine waves was to have getters/setters for the wave's amplitude, frequency, phase and offset and a function to sample the wave at a given timestamp.

Exisiting Unreal Oscillators

From what I could find, Unreal's exisiting utilities for harmonic oscillators and wave sampling are mostly used by the camera (WaveOscillatorCameraShakePattern.h's FWaveOscillator or LegacyCamerShake.h's FFOscillator) and audio (SoundModulationLFO.cpp's FLFOGenerator or WaveTableOsc.h's FWaveTableOsc) systems for camera shaking and andio synthesis/signal processing respectively.

The issue I encountered with these oscillators was their sampling functions (in actuallity called update/updateOffset) did not take in a time value from the start of the wave, but rather a delta time parameter from the last sample time. For some of the oscillators (e.g. FWaveTableOsc) the sample function (Generate) didn't even take a delta time, just two float pointers for outputting the results.

Due to this I decided to build my own implementation.

ISamplableWave Interface

I started implementation with an interface for my base wave generators. I wanted to use this interface for both the sine wave generators and the superposition wave generator, so while both would have sample functions, I had to be conservative with the getters and setters.

The sine wave generators would have getters and setters for all the wave attributes, however the super-position doesn't share such attributes, e.g. a super position of sine waves doesn't itself have a frequency. As such the interface couldn't have any attribute setters.

I also wanted to make sure that I could pass any implementation of this interface to a visualiser and not have to tailor the code to each implementation. As such I had to consider what information such a visualiser would need. The visualiser would need to know the maximum and minimum possible values the wave could produce as well as the timespan overwhich it should display the waveform. The minimum and maximum possible values are the amplitude of the wave and the timespan overwhich to display is the wavelength of the wave.

A super position of sine waves doesn't really have an amplitude in the strict sense, however its maximum/minimum possible value is the sum the maximum/minimum possible values of all the wave it's comprised of. Super-positions of sine waves do have wavelengths, however the wavelength is the lowest common multiple of the wavelengths of all the wave's it's comprised of, which could be quite large, as such I was satisfied with simply using the longest wavelength from all the waves the super-position is made of.

As such the Interface consisted of three functions: a Sample function; a GetAmplitude function; and a GetWavelength function.

USamplableWaveTable

Before I could build the super-position of sine waves, I needed to build a class to handle sampling sine waves themselves. For this it took inspiration from FWaveTableOsc in the Signal Processing part of the engine.

I started by making the USamplableWaveTable class act as it's own factory, making the constructor private and exposing a static CreateWaveTable function returning a TStrongObjectPtr to the created USamplableWaveTable. This allows me control the creation of wave tables more tightly, including allowing me to be smarter about memory management later down the line, e.g. I intend on having all wave tables share an array of pointers to their actual data as a memory optimization (this works as we always generate an unscaled/shifted wave for the table samples, though this means all tables use the same number of samples).

CreateWaveTable takes in a parent object (in order to create the new USamplableWaveTable UObject); a EWaveTableType enum value; the desired size of the wave table (how many samples the table will consist of); and the attributes of the wave. I used the code from FWaveTableOsc's CreateWaveTable for generating the wave samples as it already supports 7 types of wave and I didn't see much point in re-inventing this particular wheel.
USamplableWaveTable's CreateWaveTable also calls the setters for all of the attributes as well. The setters for the attributes mostly just set the value of a variable holding that attribute, however when setting the frequency or the phase I do a little bit of extra computation up front, to save doing it at every cycle.

As mentioned earlier, the wave generated and stored in the wave table always has and amplitude of 1; a frequency of 1; a phase of 0; and a offset of 0. We generate the wave this way because if the sample rate we use to build the wave table is a clean factor or multiple of the wave's wavelength (which is proportional to the frequency) then the values we store in the table to represent the wave can do so very poorly, e.g. if the sample rate is a multiple of half the wavelength, then every sample will read 0 and our wavetable will represent a straight line at 0.
As such we apply the wave parameters when we sample from our table.
Applying the amplitude to our sample result is simple, we multiply the value by the amplitude, as is applying the offset, we add the offset after scaling the sample by the amplitude, however applying the frequency and phase parameters is little more complex.

Rather than thinking of the frequency and phase of the wave, it's more useful for us to think of the wavelength and the offset in time the phase causes relative to a phase of 0.
Our phase parameter runs from 0 to 1 with 0 being no phase offset and 1 being a phase of $2 π$ , due to this, and the periodic nature of the wave, we can treat the phase as the percentage of the wavelength we've moved the start of the wave backwards in time. This means to apply the phase, we simply subtract the time offset $phase \times wavelength$ from the time we wish to sample.
To apply our frequency, we also use the wavelength of the wave and the fact that the wave is periodic. If our wave has a wavelength of 2 seconds and we want to sample it at 5 seconds, we don't really care about the 2 full periods of the wave that occurred before our sample time, we only care about how far through the current period of the wavelength our sample time is.
So if we mod (%) the sample time by the wavelength (ensuring the result is +ve, -3 % 4 = 1; 3 % 4 = 3) and then divide the result by the wavelength, we will end up with the percentage progress through the wave our sample is at.
We can then use this, and the number of samples we took when building the table, to find the indices of the table samples on either side of the desired sample by multiplying the number of samples by this progress. The floor of this value (rounded down) is the previous sample and the ceil (rounded up) is the next sample. We can then use the fractional part of the value to lerp the value of those samples for a more accurate approximation.

As the phase, and frequency don't change frequently, when the phase changes, we can calculate and cache the new value of the TimeOffset and when the frequency updates we calculate and cache the new value of the wavelength and then re-compute and cache the TimeOffset. This is the extra work we perform in the setters for the frequency and the phase.

UHarmonicOscillator

We now have everything needed to build our Super-Position class UHarmonicOscillator.

The UHarmonicOscillator class has two private members. An array of FSineWave UStructs which old only the frequency, amplitude and phase of waves to be exposed as a UPROPERTY to the editor for ease of editing. And an Array of USamplableWaveTable (TStrongObjectPtr)s used to actually calculate the samples.
I decided to separate these as I want to expose all of the waves being added to the super position to the user of this class, while having the liberty to optimize the waves being stored for actually calculating the samples (we do this by merging waves, for details on how see the Pre-Processing Waves section below).

The implementation of UHarmonicOscillator's Sample function is very simple, it runs through all of its USamplableWaveTables and sums the results from sampling them at the provided timestep (this is why mergine the waves is an optimization, we're reducing the number of waves we actually have to sample).

When we add a wave, we add it to our list of FSineWaves and attempt to create a new USamplableWaveTable for the wave (potentially merging it with an existing wave instead).

When we remove a wave, we can simply remove it from the list of FSineWaves, but that wave may have been merged with any number of other waves, so there might not be a matching USamplableWaveTable to remove. Instead we can use a trick of wave functions and maths, removing a number from a sum is the same as adding the negation of that number, the same goes for summing waves. If we negate the amplitude of the wave we're removing and add this new wave to the super-position it will destructively interfere with the wave we want to remove, cancelling out its effect on the result.

Pre-Processing Waves

When super-imposing sine waves we can use equivalences to reduce the work we will need to do when sampling the superposition.

Lets take the equation for a wave: $A \sin (f x + p 2 π) + O$
Where $A$ is the amplitude of the wave; $f$ is the frequency; $p$ is the phase (scaled from 0 - 1); and $O$ is the focus of the wave. Adding two waves together looks like: $A_{1} \sin (f_{1} x + p_{1} 2 π) + O_{1} + A_{2} \sin (f_{2} x + p_{2} 2 π) + O_{2}$
Immediately we can see a way to simplify this equation a little: Adding two waves together looks like: $A_{1} \sin (f_{1} x + p_{1} 2 π) + A_{2} \sin (f_{2} x + p_{2} 2 π) + O Where O = O_{1} + O_{2}$
However, if we make some assumptions, then we can simplify the equation further.
Note: I will only mention the assumtions that have results that simplify the coefficients of the wave.

If we assume that $f_{1} = f_{2} = F$ and $p_{1} = p_{2} = P$ then $\begin{matrix} A_{1} \sin (F x + P 2 π) + A_{2} \sin (F x + P 2 π) + O \\ ⇒ & (A_{1} + A_{2}) \sin (F x + P 2 π) + O \end{matrix}$
We have re-arranged the sum of the two waves into the equation for a new wave with the same frequency and phase but with amplitude $A = A_{1} + A_{2}$

If we assume that $f_{1} = f_{2} = F$ and $A_{1} = A_{2} = A$ then $\begin{matrix} A \sin (F x + p_{1} 2 π) + A \sin (F x + p_{2} 2 π) + O \\ ⇒ & A (\sin (F x + p_{1} 2 π) + \sin (F x + p_{2} 2 π)) + O \\ Given & \sin (A) + \sin (B) = 2 \sin (\frac{A + B}{2}) \cos (\frac{A - B}{2}) \\ ⇒ & A 2 \sin (\frac{F x + p_{1} 2 π + F x + p_{2} 2 π}{2}) \cos (\frac{F x + p_{1} 2 π - (F x + p_{2} 2 π)}{2}) + O \\ ⇒ & A 2 \sin (\frac{2 F x + 2 (p_{1} + p_{2}) π}{2}) \cos (\frac{p_{1} 2 π - p_{2} 2 π}{2}) + O \\ ⇒ & 2 A \sin (F x + \frac{(p_{1} + p_{2}) π}{2}) \cos (\frac{(p_{1} - p_{2}) π}{2}) + O \\ ⇒ & 2 A \cos (\frac{(p_{1} - p_{2})}{2} π) \sin (F x + \frac{(p_{1} + p_{2})}{2} π) + O \end{matrix}$
We have again re-arranged the sum into the equation for a new wave with the same frequency; amplitude $A = 2 A \cos (\frac{(p_{1} - p_{2})}{2} π)$ ; and phase $P = \frac{(p_{1} + p_{2})}{2}$

We can use these equivalences when we add a new wave to the super-position by:
1. Finding the first wave with the same frequency.
2. If that wave has the same phase then we use the first equivalence to update the wave's amplitude.
3. Otherwise, if that wave has the same amplitude then we use the second equivalence to update both the wave's amplitude and phase.
4. If we did merge the new wave with an existing wave, then if the new amplitude is 0, we can remove the altered wave as it won't contribute to the result.

This way, instead of adding a new wave to be evaluated, we merge the new wave into an exisiting one, reducing the number of sine waves we need to sample when evaluating the super-position. Additionally, if we start with a super-position of 0 waves and apply this process everytime a wave is added, our end result will use the minimum number of waves required to evaluate the super-position.