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.

**Re: Re: HTMLSave and GIF size**

**Re: How can I combinate two strings into one?**

**Re: Re: HTMLSave and GIF size**

**Converting to QuickTime**