1

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.

giopico
  • 67
  • 6
  • 2
    R has a Fast Fourier Transform built in: see the help for fft. https://stats.stackexchange.com/questions/16117 has a little more about this. You basically read the dominant frequencies off the peaks in a (smoothed) periodogram--that's all the information you need. The rest is routine coding and plotting. – whuber May 18 '22 at 21:50
  • +1 Thank you very much @whuber. I’m clear on how to find the dominant frequencies, and your comment helped a lot. However, I’m still unsure on how to retrieve the amplitude that corresponds to the frequency for the constituent harmonics. Would the amplitude for each harmonic relate to the value on the y-axis of the periodogram? – giopico May 19 '22 at 01:25
  • 1
    Yes: that's part of the definition of the periodogram. – whuber May 19 '22 at 13:43
  • Much appreciated, @whuber. Not to belabour the question, but applying the FFT to time series yt <- cos(2*pi*t/16) + 0.75*cos(2*pi*t/5) should result in a periodogram with dominant frequencies 1/16 and 1/5, and amplitudes 1 and 0.75. However, FFT <-fft(yt)and then FFT <- abs(FFT/FFT[1]) plotted into a periodogram plot(FFT, type = "l”)has none of these values. – giopico May 19 '22 at 15:07
  • http://svn.r-project.org/R/trunk/src/library/stats/R/spectrum.R has some code that looks relevant. – whuber May 19 '22 at 15:33

1 Answers1

0

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((2pi41t)/(20010.1)) + bh$bh[41]sin((2pi41t)/(20010.1))
y2 = ah$ah[14]cos((2pi14t)/(20010.1)) + bh$bh[14]sin((2pi14t)/(20010.1))
y0 = ah$ah[1]cos((2pi2t)/(20010.1)) - bh$bh[1]sin((2pi2t)/(20010.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”)
giopico
  • 67
  • 6