Frequency Analysis w/ Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg25886] Frequency Analysis w/ Mathematica?
- From: mathprof at bigfoot.com
- Date: Tue, 7 Nov 2000 02:55:57 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
I'm trying to do some frequency analysis w/ Mathematica, and am
running into some problems. I'm not sure if this is because I
don't understand Mathematica well enough, mathematics well
enough, or even if the problem I'm trying to solve is
fundamentally difficult. Consider:
x1=Table[N[197*Sin[(x-113)*2*Pi/522]+622],{x,1,10000}]
This samples a sine wave centered at y=622, amplitude 197, period
522, and phase shifted 113 units. Notice the number of samples
(10000) is not a multiple of 522, so there are a non-integral
number of periods in the sample. Similarly, the sine wave is
phase shifted, so the first sample isn't any special value (ie,
it not either extrema nor the midpoint of the sine wave). The
values 197, 113, 522, and 622 were chosen "at random" to create
an arbitrary sine wave and have no special relation/meaning
either. I also chose integers for convenience, but any/all of
these numbers could be real as well.
My problem: I'd like to reconstruct the sine wave above using
just the sampled data in the list x1. Is this possible? I would
think the problem would yield to Fourier analysis, but I haven't
been able to quite get this to work, details below. Can anyone
help?
What I've tried so far: my first approach was to use a discrete
Fourier transform:
x2=Fourier[x1]
Plotting the absolute value of the elements of x2 shows a peak at
0 (which is just the constant term 622 I think) and another peak
at 20 (which is the "frequency 19" term because of Mathematica's
convention):
ListPlot[Abs[Take[x2,{15,25}]],PlotRange->All];
This means there are 19 periods (to the nearest integer) in my
sample list; since the list has 10000 elements, there's a period
every 10000/19=526.316 elements, which is roughly correct. Of
course, since the discrete Fourier transform only works with
integers, I can't expect the result to be too accurate. I
noticed:
x2[[1]] == 62082.4 + 0. I
which I assume is from the constant term 622 in my original
function. However, it appears to be 100 times too big (do I need
to divide by the square root of 10000? I thought the discrete
Fourier transform already did this?). It also appears to be
slightly inaccurate even after dividing by 100, the result is
620.842, not 622?
I also noticed:
x2[[20]] == -7162.03 + 6181.79 I
(in other words, the peak frequency appears with this amplitude)
but wasn't sure exactly how these terms related to the 197
amplitude and 113 phase shift terms in my original function?
Since the discrete Fourier transform wasn't sufficient for my
needs, I decided to try a continuous one-- oddly, I couldn't find
a continuous Fourier transform for discrete data in Mathematica's
built-in function list (?!), so I created my own:
cft[l_,s_] := Sum[l[[x]]*Exp[2*Pi*I*(x-1)*(s-1)/Length[l]],
{x,1,Length[l]}]/Sqrt[Length[l]]
This identical to the definition of the discrete Fourier
transform given in the Mathematica book, except that s can now be
real. Plotting the absolute value of this function is slow (10000
additions per iteration!), but shows a sharp peak near 20:
Plot[Abs[cft[x1,s]],{s,15,25},PlotPoints->2,PlotRange->All]
To find the exact location of this peak, we use FindMinimum (with
2 parameters, since cft() isn't really differentiable):
FindMinimum[-Abs[cft[x1,s]],{s,19.5,21}]
This yields {s -> 20.0851} which translates to 19.0851 periods
per 10000 samples or 523.969 samples per period. Since the actual
value is 522, this isn't too bad, though not as accurate as I'd
like. I suspect the fact I don't have an integral number of
periods in my sample is throwing this off?
At s=20.0851, I show:
cft[x1,20.0851] == -8566.16 + 4261.68 I
but again, I'm not sure how to interpret this information (I'm
sure that this is somehow telling me the amplitude is 197 and the
phase shift is 113, but I'm not sure how to get from the complex
number -8566.16 + 4261.68 to the real numbers 197 and 113).
I tried several other things as well (eg, TrigFit) to no avail. I
suspect my problem is either with the phase shift or with not
having an integral number of periods in my sample. In my real
application (the above is just a demonstration of course), my
data will be coming from existing observations, so I can't
guarentee zero phase shift and/or an integral number of periods
(is there any way to take a list of data that's quasi-periodic
and get rid of phase shift and/or find an integral number of
periods? -- if the data is 100% purely periodic, this is easy,
but I suspect it's harder with real data?).
In case anyone's interested, my eventual goal is to curve-fit
observable planetary and lunar data (declination, phase,
magnitude, (modified) right ascension, etc) as a function of time
to a series of sine waves (but by sine waves, I mean phase
shifted sine waves as well, so cosine waves are implicitally
included), to provide a computationally accurate (though bizarre)
method of calculating planetary positions. Since I want these
representations to be compact, I'd like to use phase-shifted sine
waves with non-integral periods. In other words, I'd like to
write a short series of terms of the form:
197*Sin[(x-113)*2*Pi/522]
instead of a longer series of terms with just Sin[n*x], Cos[n*x].
Sincerely,
Math Prof (mathprof at bigfoot.com)
Sent via Deja.com http://www.deja.com/
Before you buy.