1
$\begingroup$

I have written some software which can perform additive synthesis given a set of partials; that is, the frequency, phase and amplitude of a set of sinusoids.

The idea is to perform an FFT on some audio data and use additive synthesis to invert the frequency data back into time-domain samples, but this doesn't seem to happen correctly.

For example, I take a two-second excerpt from a song (88,200 samples), split it into 86 "chunks" of 1024 samples (discarding the remainder samples), perform a real FFT on each chunk (numpy.fft.rfft()) to acquire 513 partials, then synthesize the time-domain signal from them, and stitch the synthesized chunks back together.

Note: I derive the amplitude of a partial by the absolute value of the complex bin value (numpy.abs()), and the phase is calculated from its anti-clockwise angle starting from the position 1+0i (numpy.angle()). Additionally, I divide the partial amplitudes by the total number of bins.

For each output time-domain sample that is being synthesized:

        pi2 = 2.0 * np.pi
        ret = 0.0
        fs1 = 1.0 / self.ctx.fs
        count = len(self.freqs)
        wfreqs = self.wfreqs[:]

        for i in range(0, count):
            ret += self.amps[i] * \
                    np.sin(self.next_phases[i] + self.init_phases[i])

        # For each partial, increment the phase to use next time
        for i in range(0, count):
            self.next_phases[i] += wfreqs[i] * fs1

            # Wrap phases (not strictly necessary I guess)
            if self.next_phases[i] >= pi2:
                self.next_phases[i] -= pi2

This context has the following variables:

  • ret: the sample at this point in time
  • fs1: 1 / sample rate
  • self.freqs: list of frequencies of the partials
  • self.wfreqs: list of angular frequencies (equal to 2.0 * pi * self.freqs)
  • self.amps: list of amplitudes of the partials
  • self.init_phases: list of initial phases (the FFT phases) of each partial
  • self.next_phases: list of phase offsets that represent the current "progress" in each partial's cycle, used by the next iteration of the loop

I was expecting to see the original signal exactly, but this is not the case. The result looks visibly different in the time-domain (top is the original):

Original signal (above) vs resynthesized

My observations:

  • The periodic peaks are the last sample from one chunk and the first of the next being at significantly different amplitudes, resulting in a rhythmic clicking
  • The amplitude is lower, for some reason
  • Ignoring these differences, the audio sounds identical by my reckoning

It's as though my understanding of what the phase data of a Fourier transform actually represents is wrong. Unless I have a bug in my code, this is way above my level of understanding. If you think this should work, and that I have a bug, I'll happily provide more.

Question: why am I not able to reverse frequency-domain data into the original time-domain samples using additive synthesis?

$\endgroup$
6
  • $\begingroup$ This question might be a good fit for dsp.se. $\endgroup$
    – J.G.
    Commented Oct 11, 2019 at 22:14
  • $\begingroup$ @J.G. I considered posting it on DSP; if you think it would get better reception there, please could you move it? (Or I could, if I knew how :) ) $\endgroup$
    – Doddy
    Commented Oct 11, 2019 at 22:17
  • $\begingroup$ I don't know how either. $\endgroup$
    – J.G.
    Commented Oct 11, 2019 at 22:26
  • 1
    $\begingroup$ The phase of the FFT bins is a measure of evenness/oddness of the time domain data for that frequency bin for that chunk, relative to the time domain origin in that chunk. Combining FFT phases across chunks doesn't quite make sense to me. $\endgroup$
    – Andy Walls
    Commented Oct 11, 2019 at 22:47
  • 1
    $\begingroup$ Also, aliasing and your use of a rectangular window (which has an infinitely wide FT) in the time domain is probably a source of some problems with what you're trying to accomplish. $\endgroup$
    – Andy Walls
    Commented Oct 11, 2019 at 22:56

1 Answer 1

1
$\begingroup$

The short answer is that, at least for a chunk, you are performing the "additive synthesis" incorrectly.

The DFT, and hence FFT, are invertible using the IDFT, and will return your original time domain signal chunk. Any additive synthesis algorithm you develop to resynthesize a chunk will be equivalent to the IDFT, and will probably not perform as well as an IFFT. So let's look at the IDFT and see if we can put it into an "additive synthesis" looking form.

Using the IDFT convention here, except using $X_n$ and $x_k$ for the frequency domain and time domain values respectively,

$$\begin{align*} x_k &= \dfrac{1}{N} \sum_{n =0}^{N-1}X_n e^{i 2\pi k n/N} \\ \\ &= \dfrac{1}{N} \sum_{n =0}^{N-1} \left|X_n\right|e^{i \mathrm{Arg}\left(X_n\right)}e^{i 2\pi k n/N}\\ \\ &= \dfrac{1}{N} \sum_{n =0}^{N-1} \left|X_n\right|e^{i \left[2\pi k n/N + \mathrm{Arg}\left(X_n\right)\right]}\\ \\ &= \dfrac{1}{N} \sum_{n =0}^{N-1} \left|X_n\right|\left(\cos\left[2\pi k n/N + \mathrm{Arg}\left(X_n\right)\right]+i\sin\left[2\pi k n/N + \mathrm{Arg}\left(X_n\right)\right]\right)\\ \\ &= \dfrac{1}{N} \sum_{n =0}^{N-1} \left|X_n\right|\cos\left[2\pi k n/N + \mathrm{Arg}\left(X_n\right)\right]+i\dfrac{1}{N} \sum_{n =0}^{N-1} \left|X_n\right|\sin\left[2\pi k n/N + \mathrm{Arg}\left(X_n\right)\right]\\ \\ &\mbox{and for known real signals,$\{x_k\}$, the imaginary terms cancel due to the symmetry of their DFT} \\ \\ x_k &= \dfrac{1}{N} \sum_{n =0}^{N-1} \left|X_n\right|\cos\left[2\pi k n/N + \mathrm{Arg}\left(X_n\right)\right] \\ \end{align*}$$

(Note that the time index $k = 0$ refers to the first sample of the time domain chunk under synthesis.)

The above is the correct formula for additive synthesis to resynthesize the time domain chunk of a real discrete signal from that chunk's DFT/FFT. It is equivalent to performing an IDFT or IFFT, but a naive, straightforward implementation of the formula will have poor performance compared to an IFFT.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .