Re: power spectrum/frequencies
- To: mathgroup at smc.vnet.net
- Subject: [mg15410] Re: power spectrum/frequencies
- From: "Atul Sharma" <mdsa at musica.mcgill.ca>
- Date: Wed, 13 Jan 1999 20:57:24 -0500
- Organization: McGill University Computing Centre
- Sender: owner-wri-mathgroup at wolfram.com
There are several points to be made about the Fourier transform you wish
to do:
You are right that the Fourier power is given by the product of the
complex Fourier series and its conjugate; by the definition of complex
conjugate, this product will be real, so you needn't use the Abs
function:
In[83]:=
pwrspec = Abs[nft cft];
pwrspec1 = nft cft;
pwrspec==pwrspec1
Out[83]=
True
However, if you are throwing out the phase data in any case, you can use
the Abs function built into Mathematica and speed up calculation of the
complex conjugate,again giving the same result:
In[84]:=
pwrspec1 = nft cft;
pwrspec2=Abs[nft]^2;
pwrspec1==pwrspec2
Out[84]=
True
While this may seem pedantic, for larger data sets, the time saving can
be significant. For example, lets create a random data set and
comparing Timing:
In[87]:=
(pwr =Fourier[rndData]*Conjugate[Fourier[rndData] ]); //Timing
pwr2=Abs[Fourier[rndData]]^2;//Timing Out[87]=
{8.4 Second,Null}
Out[88]=
{5.77 Second,Null}
(if I include the superfluous Abs function, the time is > 10s).
You are quite right about the need to discard the higher frequencies, on
account of aliasing. A quick and dirty way to do so and scale the
frequency axis of the data is given below (you will have to change the
'delta' = the frequency steps used based on your own sampling
frequency).
npts = Length[Fourier[data]]; delta = N[2 Pi/npts];
pwr=Abs[Fourier[data]]^2;
plotdata = Table[{k delta, pwr[[k]]},
{k, Floor[npts/2]}];
ListPlot[plotdata,PlotJoined->True,PlotRange->{{0,Pi},{0,2.5 10^7}}]
A. Sharma
----------------------------------------------------------------------------
----------------
Atul Sharma MD, FRCP(C)
Assistant Professor,
Department of Pediatrics,
Division of Pediatric Nephrology,
McGill University/Montreal Children's Hospital, Montreal, QC, Canada
------------------
From: jim leddon <jleddon at cyberramp.net> To: mathgroup at smc.vnet.net
Subject: [mg15410] FW: Resend: power spectrum/freqs.
Hi,
I have a question regarding calculating the power spectrum and
frequencies of a set of intensities values. The data set is a small
one.
In the program below , I have calculated and plotted not only a power
spectrum but also a power spectrum versus frequencies of a data set of
intensity values.I am trying to determine if the intensities contain
any dominant frequencies or if the power spectrum has a chaotic type
behavior shape. From what I understand, the Fourier Transform (FT)
takes a time domain set of values and converts these to a frequency
domain set of values. It seemed reasonable that once I calculated the
FT, I could account for the aliasing by removing the second part of
the spectrum and then use the remaining FT values to get the power
spectrum which is by definition the absolute value of the the
FT(remaining values) multiplied by it's conjugate (remaining values)FT
- divided by the square root of the the length of the data list. I
then calculated the frequencies based on the FT, although I've noticed
that alot of people who have used Mathematica will use the
"InverseFourier" function to calculate the power spectrum and
frequencies. My question is; in looking at the attached program, am I
correct about how to calculate the power spectrum and frequencies? Is
it also correct to use the entire FT spectrum times its' conjugate FT
(or just the square of the FTs- since the absolute value of the
conjugate is the same as the absolute value of the FT) to calculate
the power spectrum?
ClearAll;
<<Calculus`FourierTransform`
(* Program to determine power spectrum of a set of solar flux data;
Intensities in SFUs*)
data = {643,783,642,742,683,765,683,789,676,623,456,743,
783,646,666,673,873,675,673,646,676,666,763,785,783,783,
657,643,781,789,786,678,673,673,654,642,675,782,652,783,
782,767,673,672,673,675,672,673,674,675,784,783,783,666};
n = Length[data]
ft = Fourier[data];
m = Length[ft]
ListPlot[Abs[ft], PlotJoined->True]
ListPlot[data, PlotJoined->True, PlotRange->All]
(*In order to account for aliasing of fourier(FT) of data, calculate the
Power Spectrum using only the first half of the FT data*) nft =
Take[ft,27];
cft = Conjugate[nft];
pwrspec = Abs[nft cft]/ Sqrt[27];
p = Length[pwrspec]
freq = Abs[Sqrt[27] nft];
l = Length[freq]
Transpose[{freq, pwrspec}];
ListPlot[%, PlotJoined->True, PlotRange->All]
ListPlot[pwrspec, PlotJoined ->True, PlotRange->All]
I also would like to ask if there is a way to use the "Frequencies"
function under the "Statistics`DataManipulation`" package in order to
calculate the frequencies. I have used this function before, however,
there must be some way of ordering the output to correspond with a
power spectrum if one wanted to plot both the power spectrum vs.
frequencies.
Thank you for your time and consideration in reading this email message.
I'd look forward to hearing from you if you can help.
Regards,
Debbie Leddon