Spectral Analysis in R - the periodogram

Question

I am doing a spectral analysis in R using the spec.pgram() function. Suppose I have observations for $y_1$ to $y_n$ which are a time series with annual observations.
From my understanding if I run: per <- spec.pgram(y,taper=0,log="no")
a) per$spec stores the periodogram values for frequencies $1/n, 2/n, \ldots, 1/2$.
b) per$freq appears to give these frequencies.
The frequencies with large per$spec values are the dominant frequencies.

My question is:
(1) why do we use the frequencies $1/n, 2/n, \ldots, 1/2$ instead of using all the frequencies from $1/n$ to $1$. Is this simply because the periodogram is symmetric?
(2) I noticed that if the data I used is a time series with freq= 2 (has two observations per year) then the periodogram in R shows frequencies from $1/n$ to $1$ instead of $1/n$ to $0.5$. Why is this? I am a little confused what R is doing. A detailed explanation would be appreciated (for example this happens when I run spec.pgram for the sunspotz data in the astsa package).

Answer 1

The periodogram function $P(\cdot)$ has the property that $$ P(j/n) = P(1 - j/n) $$ for $j=0,\ldots,n-1$. The frequency $.5$ is a folding frequency. This all can be verified by the definition of a periodogram: $$ P(j/n) = \left(\frac{2}{n} \sum_{t=1}^n x_t \cos(2 \pi t j/n) \right)^2 + \left(\frac{2}{n} \sum_{t=1}^n x_t \sin(2 \pi t j/n )\right)^2. $$

References: http://www.stat.pitt.edu/stoffer/tsa4/

Answer 2

The Fourier transform of real-valued data is Hermitian symmetric (I usually just say conjugate symmetric), i.e., $X(f) = X^{*}(-f)$ where $*$ is conjugation and $X(f)$ is the Fourier transform of the real-valued data, $x(t)$. Since the spectrum is $S(f) = |X(f)|^{2} = X(f)X^{*}(f)$, the spectrum is symmetric about $f = 0$.

Now, when you take use the FFT in R, it is storing the positive frequencies in the 1 through (n/2)+1 indices (the 1st is $f = 0$ and the last, (n/2)+1, is the Nyquist frequency, $f_{N} = \frac{1}{2 \Delta t}$ where $\Delta t$ is the sampling rate or time between samples).

The "negative" frequencies are thus stored in the (n/2)+2 through n indices. As mentioned, the negative frequencies in the spectrum are mirrored in the positive frequencies (except 0 and Nyquist) and so normally we only plot the first (n/2) + 1 frequencies as no more information is gained. You would plot all the frequencies if your data were complex-valued however.

So, the frequencies you are estimating the spectrum at depend on the sampling rate.

I say "negative" because since the Spectrum / Fourier transform are periodic.

A couple of notes if you're interested. The spec.pgram function in R isn't actually calculating the periodogram, it's calculating the spectrum using a 10% (by default) cosine taper. There are much better tapers around, specifically the zero-th order Slepian sequence is "the best" (optimal in a broadband bias sense). The Slepian sequences are what are used in the multitaper package (more below).

It is reasonably standard practice to zero-pad your series before estimating the spectrum. This helps to take care of circular autocorrelation and takes advantage of FFT efficiency (series length with small primes as factors are optimal). You also obtain a finer frequency mesh at which the spectrum is estimated (however, no information is gained as these points are basically interpolated). A by-product of zero-padding is that the plots can often look nicer as well.

I normally just zero-pad out to a factor of two, but this will depend on application. Say M is my zero-padded length, then: $$M = 2^{\lfloor\log_{2}(n)\rfloor + 2}$$

Also, if you would like to use an optimal spectral estimator (in the bias / variance sense), there is a package called multitaper that can be used (the periodogram is horrendously biased (in a broad-band sense)). I posted a bit about it here: Power density spectrum formula in R

Hopefully that at least addressed your questions, but if not, please let me know and I'd be happy to try to help out.

Answer 3

The short answer is that yes, the periodogram is symmetric in a sense that if you get "all frequencies" instead of half, then the other half will have essentially the same information.

I'm not sure there's a simple way to explain this without some Fourier transform theory. Here's one way to think about this. The orthogonal spectral decomposition such as Fourier is used to construct the periodograms. You represent your signal as a sum of perfect waves. Each wave has an amplitude $A_i$ and a phase $\phi_i$ in addition to the frequency $\omega_i$: $$f(t)=\sum_iA_i\sin(\omega_i t+\phi_i)$$

If you have N observations of real function, you have a vector with N elements, each with one real value. When you apply something like FFT (fast Fourier transformation) you output a vector with 2 values per frequency: amplitude and phase. Hence, your output vector must be half the size of the input vector. In other words you can only go from 1/N to 1/2, not 1 because each element is now twice the size.

Why is this happening? Because Fourier transform was developed for complex functions, where each input vector entry has two values: real and imaginary parts. In fact you can think of a complex number as having an amplitude and phase, like shown on this page from Wolfram's web site. So, you need N element output to explain the complex vector with N elements, unlike N/2 elements output for real vector.