How to calculate harmonic equations from a complex wave I’m working through an online example that explains spectral analysis in R. Let’s say you have a complex wave:
t <- seq(0,200,by=0.1)
x <- cos(2*pi*t/16) + 0.75*sin(2*pi*t/5)

then the periodogram:
del<-0.1 # sampling interval
x.spec <- spectrum(x,log="no",span=10,plot=FALSE)
spx <- x.spec$freq/del
spy <- 2*x.spec$spec
plot(spy~spx,xlab="frequency",ylab="spectral density",type="l")

tells you that the frequency from the $cos(2*pi*t/16)$ harmonic dominates the wave. So, from the periodogram, you can know the frequency, and you can have a rough idea of the magnitude of the amplitude that corresponds to the frequency.
Let’s say you were only given the output of x (a time series) and you did not know the underlying harmonic formulae that compose the complex wave. Is there a way to extract the formula for, say, the 5 most dominant waveforms from a given complex time series?
I’d like to recreate something like this figure 1, where the complex waveform from the time series is visualized as the sum of individual harmonic waves. However, I’m stuck on how to retrieve the harmonic equations to visualize each harmonic individually.
 A: I think I’ve found an answer to my question - at least, an answer I’m satisfied with. Thank you to @whuber for the direction.
Given a complex wave that follows the time series (from the original post) ... :
t <- seq(0,200,by=0.1)
y <- cos(2*pi*t/16) + 0.75*sin(2*pi*t/5)

... find the formula for n harmonics that compose the signal. Let’s keep it to the first 2 harmonics for the sake of brevity.
First, we can keep the following equation in mind:
$$ y(t) = \sum_{n = 0}^{n<N/2} \left [ a_h \cdot cos(\tfrac{2 \pi n t}{N \Delta t}) + b_h \cdot sin(\tfrac{2 \pi n t}{N \Delta t})  \right ] $$
where h is the harmonic number, and t is time. For the purposes of this example, we are not interested in finding every harmonic that composes the complex wave, so we can limit the sum of the time-domain signal to the first 2 harmonics.
Now we can apply an fft() to the data (thanks again, @whuber), and proceed to extract the Fourier series sine and cosine coefficients from the FFT object using the formulae $A_h = 2*(ReFFT/N)$ and $B_h = -2*(ImFFT/N)$, where $ReFFT$ is the real component from the FFT algorithm and $ImFFT$ is the imaginary component.
The reconstructed wave kind of resembles the original complex wave, and the individual harmonics can be plotted as functions y1 and y2 in the code below.
My full workflow is shown here as a reprex. Open to suggestions or comments on any improvements that can be made!!
FFT <- fft(y)
FFT2 <- abs(FFT)/(length(FFT)/2)
plot(FFT2, xlim = c(0, length(FFT2)/2)) # Frequency domain with amplitude on y-axis, frequency on x-axis; peaks at 41Hz and 14Hz
ah = data.frame(ah = 2*Re(FFT)/length(FFT))
bh = data.frame(bh = -2*Im(FFT)/length(FFT))

y1 = ah$ah[41]*cos((2*pi*41*t)/(2001*0.1)) + bh$bh[41]*sin((2*pi*41*t)/(2001*0.1))
y2 = ah$ah[14]*cos((2*pi*14*t)/(2001*0.1)) + bh$bh[14]*sin((2*pi*14*t)/(2001*0.1))
y0 = ah$ah[1]*cos((2*pi*2*t)/(2001*0.1)) - bh$bh[1]*sin((2*pi*2*t)/(2001*0.1))

library(ggplot2)
ggplot() +
  geom_line(mapping = aes(x = t, y = y), color = "black") +
  geom_line(mapping = aes(x = t, y = y0+y1+y2), color = “blue”)

# The above ggplot looks similar to plot(t, y, “l”)

