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

# 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:

# 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:

# 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)
  • $\begingroup$ Yap. Already considered :) good idea $\endgroup$ May 20, 2018 at 21:56
  • 1
    $\begingroup$ well, you got my interest. $\endgroup$ May 21, 2018 at 2:21
  • $\begingroup$ Very interesting technical question, but if you aim to do anything creative with it based on any data that is not 100% your own.. it's most likely not really worth it. Copyright dragons will hunt down and eat their imagined share and you'd be sued broke before you know it. That's why I haven't put any of my ideas in this area into any computer. It's kind of sad because I am sure there are many people who could do lots of cool stuff with this kind of tech. $\endgroup$ May 21, 2018 at 18:22
  • $\begingroup$ @mathreadler Using the audio I linked is just for testing purposes, I intend to make music out of some random sound effects from nature I record in the forest. $\endgroup$
    – halbe
    May 21, 2018 at 20:54
  • $\begingroup$ @halbe Ok then you may be fine. Just be really careful you don't catch anyones voice on there. $\endgroup$ May 21, 2018 at 21:11

1 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.

  • $\begingroup$ +1 for interesting approach $\endgroup$ May 21, 2018 at 18:25
  • $\begingroup$ it's not done. but i wanted to spell out what you would have to do to set up the problem from the POV of sinusoidal modelling. $\endgroup$ May 21, 2018 at 19:32
  • $\begingroup$ yes i know it's not finished, but the approach is worth the upvote anyway. $\endgroup$ May 21, 2018 at 19:35
  • $\begingroup$ I don't know much math or physics outside of what you'd learn from the basic college classes so I'm having a lot of trouble following the variables here, do you know of a source I can use to understand what those greek letters mean in English in this context? Some of them seem to have quite a few scientific uses. I know how sound works and enough trigonometry to understand what I've made so far, but mostly through intuition rather than a formal understanding. $\endgroup$
    – halbe
    May 21, 2018 at 21:04
  • $\begingroup$ Admittedly I don't even really understand the fourier transform, only the inverse one (or at least that I could get the same result by adding together a bunch of sine waves which is what I did originally). I didn't think understanding it was essential to make this program but after looking at what you wrote I might be wrong about that. $\endgroup$
    – halbe
    May 21, 2018 at 21:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.