I have done a fair bit of spectral analysis, and hadn't finished collecting my
thoughts for a reply, so hadn't replied yet.
What exactly do you mean by normalize?
I have not used the functons periodogram or spectrum, however from the
description for periodogram it appears that it returns the spectral density,
which is already normalized by frequency, so you don't have to worry about
changing the appearance of your periodogram or power spectrum if you change the
time intervals of your data. With a normal Fourier Transform, not only do you
need to complex square the terms but you also need to divide by a normalizing
factor to give power-per-frequency-bin, which the periodogram appears to do.
However, if you look in various textbooks, the definition of the Fourier
Transform (FT) varies from author to author in the magnitude of the prepended
scaling factor. Since the periodogram is related to the FT (periodogram
ultimately uses the function mvfft), without examination of the code for
periodogram you cannot know the scaling factor, which is almost always one of
the following:
1.0
1/2
1/(2*pi)
SQRT(1/2)
SQRT(1/(2*pi))
In fact, if you obtain an FT (or FFT or DFT) from a piece of electronics (say
an electronic spectrum analyzer), the prepending factor can vary from
manufacturer to manufacturer.
Fortunately, there is a strict relationship between the variance of your
signal and the integrated spectral density. If your time signal is x(t), the
spectral density is S(f), and fc = frequency(x)/2 the Nyquist cutoff frequency,
then this may be expressed as:
variance(x(t)) = constant * {Integral from 0 to fc of S(f)}
In R-code: let x be your time series, and "constant" be the unknown scaling
factor (1/2, 1/2pi, etc.)
p <- periodogram(x)
var(x) == constant * sum(p[[1]])/length(p[[1]])
Or:
constant = sum(p[[1]])/length(p[[1]])/var(x)
and we find that the appropriate scaling constant is 1.0.
As regards plotting versus period, the periodogram returns arrays of spectral
amplitude and frequencies. The frequencies are in inverse units of the
intervals of your time series. i.e. if your time series is 1-point per day,
then the frequencies are in 1/day units. Thus if you wish to plot amplitude
versus period in weeks you have a little math to do.
I believe that plotting is usually versus frequency since most observers are
interested in how things vary versus frequency: multiple evenly spaced peaks on
a linear frequency scale indicates the presence of harmonics; this is not so
simply seen in a plot versus period. ex. peaks of 10 Hz, 20 Hz, 30 Hz, 40
Hz,... in period would be at periods of 100 ms, 50 ms, 25 ms, 12.5 ms, and the
peaks are not evenly spaced.
Additionally, there are all kinds of typical responses versus frequency (1/f,
1/f^2) which are clearly seen in plots versus frequency as straight lines (log
power vs log frequency), but would come out as curves in plots versus period.
I can see how ecological studies may indeed be more interested in the periods.
However, I would be wary of using the periodogram function, for if I
calculate periodograms of the same sinewave but for differing lengths of the
sample, the spectral density does not come out the same. All 4 of the plots
produced by the code below should overlay, and yet as the time series becomes
longer there appears to be an increasing offset of the magnitudes returned.
(black-brown-blue-red)
x0<-ts(data=sin(2*pi*1.1*(1:1000)/10),frequency=10)
p0<-periodogram(x0)
var(x0)
x1<-ts(data=sin(2*pi*1.1*(1:10000)/10),frequency=10)
p1<-periodogram(x1)
var(x1)
x2<-ts(data=sin(2*pi*1.1*(1:100000)/10),frequency=10)
p2<-periodogram(x2)
var(x2)
x3<-ts(data=sin(2*pi*1.1*(1:1000000)/10),frequency=10)
p3<-periodogram(x3)
var(x3)
plot(p3[[2]],p3[[1]],col="red",type="p",pch=19,cex=0.05,log="y",
xlim=c(0.1,0.12),ylim=c(1e-30,1e-15))
points(p2[[2]],p2[[1]],type="l",col="blue")
points(p1[[2]],p1[[1]],type="l",col="brown")
points(p0[[2]],p0[[1]],type="o",col="black")
> Message: 113
> Date: Mon, 30 Jan 2006 17:45:58 -0800
> From: Spencer Graves <[EMAIL PROTECTED]>
> Subject: Re: [R] How do I "normalise" a power spectral density
> analysis?
> To: Tom C Cameron <[EMAIL PROTECTED]>
> Cc: r-help@stat.math.ethz.ch
> Message-ID: <[EMAIL PROTECTED]>
> Content-Type: text/plain; charset=us-ascii; format=flowed
>
> Since I have not seen a reply to this post, I will
> offer a comment,
> even though I have not used spectral analysis myself and
> therefore have
> you intuition about it. First, from the definitions I read in the
> results from, e.g., RSiteSearch("time series power spectral density")
> [e.g.,
> http://finzi.psych.upenn.edu/R/library/GeneTS/html/periodogram
> .html] and
> "spectral analysis" in Venables and Ripley (2002) Modern Applied
> Statistics with S (Springer), I see no reason why you
> couldn't plot the
> spectrum vs. the period rather than the frequency. Someone else may
> help us understand why it is usually plotted vs. the frequency; I'd
> guess that the standard plot looks more like the integrand in the
> standard Fourier inversion formula, but I'm not sure.
>
> If you'd like more help from this listserve, you
> might briefly
> describe the problem you are trying to solve, why you think spectral
> analysis analysis should help, and include a toy example with some
> self-contained R code to illustrate what you tried and what you don't
> understand about it. (And PLEASE do read the posting guide!
> "www.R-project.org/posting-guide.html". Nothing is certain but
> following that posting guide will, I believe, tend to
> increase the speed
> and utility of response.)
>
> hope this helps.
> spencer graves
>
> Tom C Cameron wrote:
>
> > Hi everyone
> >
> > Can anyone tell me how I normalise a power spectral density
> (PSD) plot of a
> > periodical time-series. At present I get the graphical
> output of spectrum VS
> > frequency.
> >
> > What I want to acheive is period VS spectrum? Are these the
> same things but the
> > x-axis scale needs transformed ?
> >
> > Any help would be greatly appreciated
> >
> > Tom
> >
> ..............................................................
> .............
> > Dr Tom C Cameron office: 0113 34
> 32837 (10.23 Miall)
> > Ecology & Evolution Res. Group. lab: 0113 34 32884
> (10.20 Miall)
> > School of Biological Sciences Mobile: 07966160266
> > University of Leeds email:
> [EMAIL PROTECTED]
> > Leeds LS2 9JT
> > LS2 9JT
> >
> > ______________________________________________
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
