Feedback FM with bells and whistles

Single oscillator feedback FM is a most economic technique for producing rich harmonic tones. However, the technique suffers from parasitic oscillations at the Nyquist frequency when the modulation index is turned up sufficiently high. The most obvious thing to try is to modify the original formula

x[n] = sin(ωn + Ix[n-1])

by lowpass filtering the feedback signal with some filter that has a zero at half the sample rate. A two point average increases the range the index can take before the spurious oscillations set in, but it cannot stop them at sufficiently high modulation indices. A complementary trick is to put a filter outside the feedback loop. Again, it helps to a certain extent, but should not be expected to solve the problem in all cases. Finally, there is the Overkill Solution of oversampling the system. Or maybe it's not so overkill after all. In any case, a high sample rate is recommended.

Feedback FM waveform with spurious oscillations and no attempt to squelch them.

Depending on the sample rate, the spurious oscillations will typically ring at such a high frequency that only domestic animals will notice them and possibly object to their presence. Nonetheless, the waveform will be contaminated by a conspicuous wiggle that may annoy any purist, whether or not they hear it. Of course, it is a spurious or parasitic oscillation since it follows the sampling rate rather than the synthesis parameters. Sometimes the parasitic oscillation happens at a third or fourth of the sample rate, or other subharmonics.

Single oscillator feedback FM is limited to harmonic spectra. Much flexibility is gained  by introducing a second oscillator since there are several ways to connect the two oscillators. Rather than listing all cases separately, we introduce coupling parameters c (cross terms) and b (self-modulation) in the coupled system

      x[n+1] = sin(θ[n] + b₁x[n] + c₁y[n])

(*)   y[n+1] = sin(φ[n] + b₂y[n] + c₂x[n])

where the phases θ and φ are incremented by the modulating and carrier frequencies. Which one is which depends on what signal you send to the output. Of course both signals can be used for stereo, but then it makes less sense to call one of the oscillators the carrier and the other one the modulator.

In FM synthesis, the phase variables usually depend only on their respective frequencies. By introducing an interaction term, phase coupling can be used to synchronize the oscillators. Hardsync may have been used with FM before, but the gentler kind of sync used in the Kuramoto model is useful here, as it is also suitable for synchronizing more than two oscillators. Now, the phases are incremented by the oscillators' frequencies as usual, but to that we add a coupling term with strength K:

      θ[n+1] = θ[n] + ω_c - K sin(θ[n] - φ[n])
(#)   φ[n+1] = φ[n] + ω_m - K sin(φ[n] - θ[n])

Turning up K too much will collapse the two oscillators into a single strong team working in perfect sync. The system (*, #) is just a four-dimensional map with seven parameters and may thus be studied with the appropriate methods of dynamic systems.

As often happens with iterated maps, the feedback FM system exhibits typical behaviour such as period two oscillations, and the period doubling route to chaos. The frequency terms may both be set to zero, which means the system (*) becomes autonomous. Then period doublings can be seen more easily, as shown below. The system has to be seeded with nonzero initial values for any oscillation to occur.

Colour legend: gray, period 1; orange, period 2; blue, period 4, bright yellow, period 3, red, chaos. On the horizontal axis, -1 < c₁ = c₂ < 4, vertical axis: - 0.5 < b₁ < 5 and b₂ = 0.5 (constant). The modulator and carrier frequencies are both 0; hence the coupling term does not influence the dynamics.

We are not done with feedback FM!

Other things to try include modifying the feedback signals by waveshapers and filters. Even the phase coupled signal may be filtered. Each filter adds its state variables to the system and increases its dimension. This is complicated territory; suffice it to say that filters plus strong nonlinearities or high feedback gain equals noise!

Oscilloscopes, phase plots and beyond

Similarly to how oscilloscopes trace out two signals x(t) and y(t), one might plot any related variables against each other.

The above illustration shows x(t) against y(t), both of which seem to be periodic, albeit quite complicated waveforms.

Phase plots of, say, position versus momentum are another common way to display orbits of differential equations. Any recorded signal may be plotted with its numerically estimated derivative on the y-axis. Hint: use a higher order differentiating filter, not just a forward or backward difference.

Graphs of the amplitude spectrum of one signal against that of another signal are harder to interpret. The axes then represent amplitudes, and each point is a unique frequency whose coordinates is the amplitude in signal x versus its amplitude in signal y. If the points gather near a thin line along the main diagonal, the two spectra are highly correlated.

Another novel kind of plot is the signed amplitude spectrum. It combines phase information with the amplitude spectrum and distinguishes positive and negative amplitues depending on the sign of the phase. Below, the red line is the signed amplitude spectrum of one variable and the blue line is that of a related variable in the same system.

All the illustrations come from a 3D ODE that is currently under investigation, suspected guilty of exhibiting interesting behaviour.

Filtering with differential equations

For those of us who are more familiar with digital filters than their analog counterparts, a one pole lowpass filter is easy:

y_n = (1 - β)x_n + βy_n-1,   0 < β < 1.

But how do you filter a signal with an ordinary differential equation?



Working backwards, we should have an ODE that says

dy/dt + Ay(t) = Bx(t)   (*)

for some suitable constants A and B yet to be found. The differential may be approximated with a forward difference over a short time interval T, so dy(nT)/dt ≈ (y(nT+T) - y(nT))/T. Setting T equal to one sampling period, the discrete time version of (*) is

(y_n+1 - y_n)/T + ay_n = bx_n. 

Perform some algebraic shuffling to and fro of the variables to obtain

y_n+1 = bTx_n + (1 - aT)y_n

and recall that the filter coefficients are 1 - β = bT and β = 1 - aT, hence bT = aT. Now we introduce a new coefficient τ > 0 for the variables in (*) and set 1/τ equal to both A and B. The system then is

dy/dt = (x - y) / τ

where now τ plays the role of a relaxation time constant. The greater τ is, the slower the response of the filter. Also, when the input equals the output the derivative becomes zero, which is to say that the system has unit DC response as required.