Re: Wanting to take a FFT transform of a irregularly spaced sample
"ashwin .D" <[email protected]> Wed, 13 Apr 2022 20:02:45 +0530
| Newsgroups | gmane.comp.python.scientific.user |
|---|---|
| Message-ID | <CAH0LXy7bN=Cqp=mVD2Zk3yo0DPep1OVue7rt6JSbGG9+xSzP7w@mail.gmail.com> |
--===============2680244645832295088== Content-Type: multipart/alternative; boundary="000000000000b0396505dc8a0d02" --000000000000b0396505dc8a0d02 Content-Type: text/plain; charset="UTF-8" Based on this github issue I had to make one modification - https://github.com/scipy/scipy/issues/13812 as shown below. I had to remove the zero frequency signal.lombscargle(minutes,df1.to_numpy().ravel(),freq[1:]) and the code works just fine now. The interpretation task awaits me now . For the record here is the code - times = df.view(np.int64)/6E10 df1 = (data[['Pressure_hPa']]) # I also added a 2 to the angular frequency value freq = np.linspace(0, 2*np.pi/300.0,8928//2) signal.lombscargle(times,df1.to_numpy().ravel(),freq[1:]) On Wed, Apr 13, 2022 at 7:46 PM ashwin .D <[email protected]> wrote: > Robert, > I ran that code as you suggested. I get this error now - > File "ffttest.py", line 21, in <module> > signal.lombscargle(minutes,df1.to_numpy().ravel(),freq) > File > "/usr/local/lib/python3.8/dist-packages/scipy-1.7.1-py3.8-linux-x86_64.egg/scipy/signal/spectral.py", > line 150, in lombscargle > pgram = _lombscargle(x, y, freqs) > ZeroDivisionError: () > > Any suggestions would be appreciated. Is this an issue I should open on > github by providing the test data ? > > Best regards, > Ashwin. > > > > On Tue, Apr 12, 2022 at 6:48 PM Robert Kern <[email protected]> wrote: > >> On Tue, Apr 12, 2022 at 3:32 AM ashwin .D <[email protected]> wrote: >> >>> Hi Robert, >>> Thanks for your prompt response. I am going to try >>> both. Regarding this answer that you recommended - >>> https://stackoverflow.com/questions/34428886/discrete-fourier-transformation-from-a-list-of-x-y-points/34432195#34432195 >>> >>> what would be my angular frequencies from this API - >>> https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lombscargle.html >>> ? >>> >>> The x and y arguments are straightforward and are available to me from >>> the CSV file. What about the third one ? >>> >> >> That's the angular frequencies at which you want to evaluate the >> periodogram at. In your case (otherwise-regular time series but with >> missing values), I would recommend using the angular frequencies that you >> would have had if you had computed a normal periodogram using the FFT on >> the whole time series, e.g. `np.linspace(0, np.pi/300.0, 8928//2)` >> (assuming your `x` is in seconds). The running time is >> O(len(x)*len(freqs)), though, so that may take a long time. You may want to >> reduce the number of points you sample at first for visualization, then you >> can zoom in at the full frequency resolution to an area of interest if >> there is lots of dead space. >> >> -- >> Robert Kern >> _______________________________________________ >> SciPy-User mailing list -- [email protected] >> To unsubscribe send an email to [email protected] >> https://mail.python.org/mailman3/lists/scipy-user.python.org/ >> Member address: [email protected] >> > --000000000000b0396505dc8a0d02 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr"><div>Based on this github issue I had to make one modifica= tion -=C2=A0<a href=3D"https://github.com/scipy/scipy/issues/13812">https:/= /github.com/scipy/scipy/issues/13812</a></div><div>as shown below. I had to= remove the zero frequency=C2=A0</div>signal.lombscargle(minutes,df1.to_num= py().ravel(),freq[1:])<div><br></div><div>and the code works just fine now.= The interpretation task awaits me now .=C2=A0</div><div><br></div><div>For= the record here is the code -=C2=A0</div><div><br></div><div>times =3D df.= view(np.int64)/6E10<br><br><br>df1 =3D (data[['Pressure_hPa']])<br>= <br># I also added a 2 to the angular frequency value=C2=A0</div><div><br>f= req =3D np.linspace(0, 2*np.pi/300.0,8928//2)</div><div><br>signal.lombscar= gle(times,df1.to_numpy().ravel(),freq[1:])<br></div><div><br><div><br></div= ><div><br></div></div></div><br><div class=3D"gmail_quote"><div dir=3D"ltr"= class=3D"gmail_attr">On Wed, Apr 13, 2022 at 7:46 PM ashwin .D <<a href= =3D"mailto:[email protected]">[email protected]</a>> wrote:<br></div><= blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-l= eft:1px solid rgb(204,204,204);padding-left:1ex"><div dir=3D"ltr">Robert,<d= iv>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 I ran that code as you = suggested. I get this error now -=C2=A0</div><div>File "ffttest.py&quo= t;, line 21, in <module><br>=C2=A0 =C2=A0 signal.lombscargle(minutes,= df1.to_numpy().ravel(),freq)<br>=C2=A0 File "/usr/local/lib/python3.8/= dist-packages/scipy-1.7.1-py3.8-linux-x86_64.egg/scipy/signal/spectral.py&q= uot;, line 150, in lombscargle<br>=C2=A0 =C2=A0 pgram =3D _lombscargle(x, y= , freqs)<br>ZeroDivisionError: ()<br></div><div><br></div><div>Any suggesti= ons would be appreciated. Is this an issue I should open on github by provi= ding the test data ?=C2=A0</div><div><br></div><div>Best regards,</div><div= >Ashwin.=C2=A0</div><div><br></div><div><br></div></div><br><div class=3D"g= mail_quote"><div dir=3D"ltr" class=3D"gmail_attr">On Tue, Apr 12, 2022 at 6= :48 PM Robert Kern <<a href=3D"mailto:[email protected]" target=3D"_= blank">[email protected]</a>> wrote:<br></div><blockquote class=3D"g= mail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204= ,204,204);padding-left:1ex"><div dir=3D"ltr"><div dir=3D"ltr">On Tue, Apr 1= 2, 2022 at 3:32 AM ashwin .D <<a href=3D"mailto:[email protected]" targ= et=3D"_blank">[email protected]</a>> wrote:<br></div><div class=3D"gmai= l_quote"><blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8e= x;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir=3D"ltr"= >Hi Robert,<div>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2= =A0 Thanks for your prompt response. I am going to try both. Regarding this= answer that you recommended -=C2=A0<a href=3D"https://stackoverflow.com/qu= estions/34428886/discrete-fourier-transformation-from-a-list-of-x-y-points/= 34432195#34432195" target=3D"_blank">https://stackoverflow.com/questions/34= 428886/discrete-fourier-transformation-from-a-list-of-x-y-points/34432195#3= 4432195</a></div><div><br></div><div>what would be my angular frequencies f= rom this API -=C2=A0<a href=3D"https://docs.scipy.org/doc/scipy/reference/g= enerated/scipy.signal.lombscargle.html" target=3D"_blank">https://docs.scip= y.org/doc/scipy/reference/generated/scipy.signal.lombscargle.html</a> ?=C2= =A0</div><div><br></div><div>The x and y arguments are straightforward and = are available to me from the CSV file. What about the third one ?=C2=A0</di= v></div></blockquote><div><br></div><div>That's the angular frequencies= at which you want to evaluate the periodogram at. In your case (otherwise-= regular time series but with missing values), I would recommend using the a= ngular frequencies that you would have had if you had computed a normal per= iodogram using the FFT on the whole time series,=C2=A0e.g. `np.linspace(0, = np.pi/300.0, 8928//2)` (assuming your `x` is in seconds). The running time = is O(len(x)*len(freqs)), though, so that may take a long time. You may want= to reduce the number of points you sample at first for visualization, then= you can zoom in at the full frequency resolution to an area of interest if= there is lots of dead space.</div></div><div><br></div>-- <br><div dir=3D"= ltr">Robert Kern</div></div> _______________________________________________<br> SciPy-User mailing list -- <a href=3D"mailto:[email protected]" target= =3D"_blank">[email protected]</a><br> To unsubscribe send an email to <a href=3D"mailto:[email protected]= rg" target=3D"_blank">[email protected]</a><br> <a href=3D"https://mail.python.org/mailman3/lists/scipy-user.python.org/" r= el=3D"noreferrer" target=3D"_blank">https://mail.python.org/mailman3/lists/= scipy-user.python.org/</a><br> Member address: <a href=3D"mailto:[email protected]" target=3D"_blank">win= [email protected]</a><br> </blockquote></div> </blockquote></div> --000000000000b0396505dc8a0d02-- --===============2680244645832295088== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ SciPy-User mailing list -- [email protected] To unsubscribe send an email to [email protected] https://mail.python.org/mailman3/lists/scipy-user.python.org/ Member address: [email protected] --===============2680244645832295088==--