Morphing between two FFTs

Question

I'm a musician and I'm making a script that takes a wave file and snaps each of its frequencies from the fourier transforms to the nearest musical harmonic. Thanks to help from another question I posted here, that part works, but what I need to do now is to experiment with some interesting ways to interpolate between the before and after version. I don't want to just do a crossfade because I want it to sound like one sound continuously morphing rather than one fading out and in. One idea I have is to snap one frequency at a time, maybe either in order based on the volume of that frequency by the end, from the lowest frequencies to the highest, highest to lowest, or random. Another idea was to interpolate them like image interpolation where if I have a frequency at 100hz with a volume of 255 and another at 120hz with a volume of 215, halfway between the interpolation I should have a frequency at 110hz with a volume of 235. If anyone knows of a function library that lets me do these then please let me know, otherwise I'll be making them and will need some help on getting started. How would I accomplish either of these?

Another thought I had was to reduce the number of harmonics over time so that in the beginning there are enough harmonics for the sound to be indistinguishable from the original but by the end there are only a couple. To get what I mean by the number of harmonics see the part of the script below where I calculate valid_ratios and valid_frequencies.

This is the sample I'm testing with (rename it to missile.wav): https://my.mixtape.moe/iltlos.wav

Here is the script so far:

import numpy as np
import scipy.io.wavfile


samples, data = scipy.io.wavfile.read("missile.wav")


# fast fourier transform makes an array of the frequencies of sine waves that comprise the sound
data_fft = np.fft.rfft(data)


# the higher MAX_HARMONIC is, the more it sounds like the original, 
# the lower it is, the more it sounds like an instrument
MAX_HARMONIC = 2

# generate list of ratios that can be used for tuning (not octave reduced)
valid_ratios = []
for i in range(1, MAX_HARMONIC + 1):
    for j in range(1, MAX_HARMONIC + 1):
        if i % 2 != 0 and j % 2 != 0:
            valid_ratios.append(i/float(j))
            valid_ratios.append(j/float(i))


# remove dupes
valid_ratios = list(set(valid_ratios))


# find all the frequencies with the valid ratios
valid_frequencies = []
multiple = 2
while(multiple < samples):
    multiple *= 2

    for ratio in valid_ratios:
        frequency = ratio * multiple

        if frequency < samples:
            valid_frequencies.append(frequency)


# remove dupes and sort and turn into a numpy array
valid_frequencies = np.sort(np.array(list(set(valid_frequencies))))


# bin the data_fft into the nearest valid frequency
valid_frequencies = valid_frequencies.astype(np.int64)
boundaries = np.concatenate([[0], np.round(np.sqrt(0.25 + valid_frequencies[:-1] * valid_frequencies[1:])).astype(np.int64)])
select = np.abs(data_fft) > 1
filtered_data_fft = np.zeros_like(data_fft)
filtered_data_fft[valid_frequencies] = np.add.reduceat(np.where(select, data_fft, 0), boundaries)


# do the inverse fourier transform to get a sound wave back
recovered_signal = np.fft.irfft(filtered_data_fft)


# amplify to ceiling
recovered_signal /= np.amax(recovered_signal)

scipy.io.wavfile.write("missile_output.wav", samples, recovered_signal)

Answer 1

maybe you can do it from bottom to top or top to bottom using something like the following.

---

a notation about notation:

$x[n]$ is the input sample for the $n$-th sample. this originally comes from an A/D converter, but may have come from a sound file. $x_1[n]$ and $x_2[n]$ are two separate inputs.

$y[n]$ is the output sample at the $n$-th time. it will eventually go to a D/A converter.

given subscripts $_1$ or $_2$ for inputs $x_1[n]$ and $x_2[n]$, $a_k[n]$ is the amplitude of the $k$-th harmonic or the $k$-th sinusoidal partial (if the sound is not very periodic and has non-harmonic overtones or partials. $\theta_k[n]$ is the instantaneous angle of the sinusoid that is the $k$-th partial, it's dimensionless but we can think of it as in units of "radians". it is incremented at a rate of $\omega_k[n]$ (radians per sample), which is the instantaneous frequency of the $k$-th partial.

likewise for the output $y[n]$, $b_k[n]$ is the amplitude of the $k$-th harmonic or the $k$-th sinusoidal partial, $\phi_k[n]$ is the instantaneous angle of the sinusoid that is the $k$-th partial, and $\nu_k[n]$ is the instantaneous frequency of the $k$-th partial of the output $y[n]$.

---

i am being lazy and not specifying the limits of $k$. how many sinuosoidal harmonics or partials that you keep tabs on is a performance issue that is up to you. too few and your sounds will be incomplete. too many and you might have computational issues. you don't want to treat every little spike or blip in the FFT spectrum as a bona fide "partial". deal with the biggest peaks first.

for sound 1:

$$ x_1[n] = \sum_k a_{1,k}[n] \cos(\theta_{1,k}[n]) $$ where $$ \theta_{1,k}[n] = \theta_{1,k}[n-1] + \omega_{1,k}[n] $$

and sound 2:

$$ x_2[n] = \sum_k a_{2,k}[n] \cos(\theta_{2,k}[n]) $$ where $$ \theta_{2,k}[n] = \theta_{2,k}[n-1] + \omega_{2,k}[n] $$

and for the output:

$$ y[n] = \sum_k b_k[n] \cos(\phi_k[n]) $$ where $$ \phi_k[n] = \phi_k[n-1] + \nu_k[n] $$

and let's say you have the harmonics ordered with $k$ the same for both inputs and output. so you want $a_{1,k}[n]$ to morph to $a_{2,k}[n]$ and $\omega_{1,k}[n]$ to morph to $\omega_{2,k}[n]$ for all of the $k$.

so the amplitude envelope and frequency contours would have crossfade functions might look like:

$$b_k[n] \ = \ (1-g_k[n])\,a_{1,k}[n] \ + \ g_k[n] \,a_{2,k}[n]$$ or perhap you wanna do it to the logarithms of the amplitude, i dunno: $$\log(b_k[n]) \ = \ (1-g_k[n])\,\log(a_{1,k}[n]) \ + \ g_k[n] \,\log(a_{2,k}[n])$$

similarly for frequency:

$$\nu_k[n] \ = \ (1-h_k[n])\,\omega_{1,k}[n] \ + \ h_k[n] \,\omega_{2,k}[n]$$

or log of frequency:

$$\log(\nu_k[n]) \ = \ (1-h_k[n])\,\log(\omega_{1,k}[n]) \ + \ h_k[n] \,\log(\omega_{2,k}[n])$$

so your question is, assumeably, how to define $g_k[n]$ and $h_k[n]$ so that the morph sounds like what you want it to. is that correct?

---

now we'll introduce you to another simple (but conceptually salient) function that we audio guys call "clipping" or "limiting" but the computer graphics people call clamping:

$$ \operatorname{clamp}(x) \triangleq \begin{cases} 0 \qquad & x < 0 \\ x \qquad & 0 \le x < 1 \\ 1 \qquad & 1 \le x \\ \end{cases} $$

---

okay, so a possible way to morph this that might sound cool is to define the crossfading functions like this:

$$ g_k[n] = \operatorname{clamp}\Big( \alpha(n-m) - \beta(k-1) \Big)$$

where $\alpha>0$ is your change rate and $\beta\ge0$ is the delay for the next harmonic to kick in.. $m$ is the sample index when the morph begins. component number 1 will begin to change at sample $n=m$, but component number 2 will not change until sample $n=m+\beta$ and component 3 will not begin to change until time $n=m+2\beta$. component $k$ will not begin to change until time $n=m+(k-1)\beta$. do the same thing for $h_k[n]$. this is the bottom frequencies changing before the higher frequencies. you can make the top frequencies change first before the bottom frequencies by changing the sign on $\beta$.

you can also apply the Smoothstep function to the result of the $\operatorname{clamp}(\cdot)$ function.