pySTEPS / pysteps

Python framework for short-term ensemble prediction systems.
https://pysteps.github.io/
BSD 3-Clause "New" or "Revised" License
466 stars 168 forks source link

rapsd returns a negative frequency for even-sized images #208

Closed mjwillson closed 3 years ago

mjwillson commented 3 years ago
In [109]: _, freqs = rapsd(np.zeros((256, 256)), return_freq=True)
In [110]: freqs[-1]
Out[110]: -0.5

Was this intentional? Note this doesn't happen for odd sizes:

In [113]: _, freqs = rapsd(np.zeros((255, 255)), return_freq=True)
In [114]: all(f >= 0 for f in freqs)
Out[114]: True

It looks like this code is responsible -- was there a reason for the +1 here with even sizes?

    if L % 2 == 0:
        r_range = np.arange(0, int(L / 2) + 1)
    else:
        r_range = np.arange(0, int(L / 2))

Note that plot_spectrum1d is I think ignoring this negative frequency as it'll be mapped to NaN after taking logs.

loforest commented 3 years ago

Tanks a lot for pointing this out! Fourier transform does return both positive and negative frequencies due to symmetry, but the subsequent application of rapsd should only return positive frequencies. It seems that the range above is too wide and "overflew" into the negative frequencies.
I will check the code to see what happens. My first feeling is that it can be solved by stopping the range for even sizes at int(L/2)-1, but I should check where is the DC term to ensure there is the correct correspondence.

pulkkins commented 3 years ago

https://github.com/pySTEPS/pysteps/blob/76324d8f315f63c6723887f4c99d155749a31e83/pysteps/utils/spectral.py#L149 I believe there is a typo. It should be if L % 2 == 1:

loforest commented 3 years ago

Ok. In that case the docstring return section should be changed to: The length of the array is int(L/2) (if L is even) or int(L/2)+1 (if L is odd), where L=max(M,N). Correct?