[R] How do I "normalise" a power spectral density
Gabor Grothendieck
ggrothendieck at gmail.com
Tue Jan 31 22:44:21 CET 2006
This post has some details:
http://tolstoy.newcastle.edu.au/~rking/R/help/04/10/5581.html
On 1/31/06, Leif Kirschenbaum <leif at reflectivity.com> 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 <spencer.graves at pdf.com>
> > Subject: Re: [R] How do I "normalise" a power spectral density
> > analysis?
> > To: Tom C Cameron <bgytcc at leeds.ac.uk>
> > Cc: r-help at stat.math.ethz.ch
> > Message-ID: <43DEC156.2090900 at pdf.com>
> > 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:
> > t.c.cameron at leeds.ac.uk
> > > Leeds LS2 9JT
> > > LS2 9JT
> > >
> > > ______________________________________________
> > > R-help at 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 at 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
>
More information about the R-help
mailing list