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[[&#39;Pressure_hPa&#39;]])<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 &lt;<a href=
=3D"mailto:[email protected]">[email protected]</a>&gt; 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 &quot;ffttest.py&quo=
t;, line 21, in &lt;module&gt;<br>=C2=A0 =C2=A0 signal.lombscargle(minutes,=
df1.to_numpy().ravel(),freq)<br>=C2=A0 File &quot;/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 &lt;<a href=3D"mailto:[email protected]" target=3D"_=
blank">[email protected]</a>&gt; 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 &lt;<a href=3D"mailto:[email protected]" targ=
et=3D"_blank">[email protected]</a>&gt; 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&#39;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==--