This post has some details:
http://tolstoy.newcastle.edu.au/~rking/R/help/04/10/5581.html
On 1/31/06, Leif Kirschenbaum <[EMAIL PROTECTED]> wrote:
> 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: [email protected]
> > 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
> > >
> > > ______________________________________________
> > > [email protected] mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide!
> > http://www.R-project.org/posting-guide.html
> >
>
> ______________________________________________
> [email protected] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
