derive an average waveform from an FFT spectrum, where all cycles have same shape and period
I found the answer here:

https://www.dsprelated.com/thread/3253/simple-and-accurate-phase-estimation-of-harmonics-in-a-signal

Just change atan() to atan2() to account for the correct phase in all quadrants.