Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding idft #129

Merged
merged 35 commits into from
Jan 28, 2021
Merged

Adding idft #129

merged 35 commits into from
Jan 28, 2021

Conversation

lanougue
Copy link
Contributor

@lanougue lanougue commented Dec 4, 2020

Adding new function idft to do inverse Fourier Transform

lanougue and others added 4 commits November 10, 2020 16:10
update for coming fft and ifft modifications
Moving spacing as coords attribute + some update (xgcm#126)
@lanougue
Copy link
Contributor Author

lanougue commented Dec 4, 2020

Here is an addition of the idft function. True_phase option is also available to do correct phase inversion. A new argument "lag" is also added to allow the user to have arbitrary output coordinates.

I can add tests. What kind of tests do you think are relevant ?
I think that the two tests below should be done:

  • idft(dft(S)) == S
  • idft(dft(S, true_phase=True), true_phase=True) == S

@roxyboy
Copy link
Member

roxyboy commented Dec 7, 2020

I think that the two tests below should be done:

  • idft(dft(S)) == S
  • idft(dft(S, true_phase=True), true_phase=True) == S

I agree with this.

xrft/xrft.py Outdated Show resolved Hide resolved
@rabernat rabernat mentioned this pull request Dec 9, 2020
@lanougue
Copy link
Contributor Author

lanougue commented Dec 10, 2020

I am almost done,
I still have a question on how to resolve some tests. I want to check what happen if a user try to do: ifft(fft(s))

If the user try idft(dft(s)) with true_phase and true_amplitude set to True for both transformations, he will recover the original signal "s". No problem.

However, ifft(fft(s)) will almost never be equal to s. This is because initial coordinates of "s" (and thus the phase) have been lost during the operation. Thus, coordinates of "s" and "ifft(fft(s))" will never match (it will match only if "s" coordinates are centered (meaning phase=0). However, s.data and ifft(fft(s)).data match as with numpy.
(shifting operation can yet alter this equality).

I think this is a totally normal behaviour, but can be disturbing for a user who always tough that ifft(fft(s))==s. It is true for amplitude but this is not the case anymore for the coordinates unless they are correctly taken into account during both transformations (as with dft() and true_phase flag).

In an xarray framework, the unique option to have ifft(fft(s))==s for both amplitude and coordinates is to convince people to abandon fft to go to dft.

@roxyboy
Copy link
Member

roxyboy commented Dec 10, 2020

However, s.data and ifft(fft(s)).data match as with numpy.

Thanks for your work @lanougue ! I'm in favor of s.data = ifft(fft(s)).data as that is what one would also expect from numpy.

@roxyboy
Copy link
Member

roxyboy commented Dec 10, 2020

@lanougue I think you can have xrft.ifft raise a warning about the coordinate information getting lost in the operation and recommend the usage of xrft.idft.

@lanougue
Copy link
Contributor Author

Done.
In all functions I change dft calls into fft calls to keep consistency with current xrft behaviour.
It could be interesting to let the user choose between fft or dft call in other functions such as power_spectrum, cross_spectrum, ...

Copy link
Member

@roxyboy roxyboy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lanougue I think these changes are great! My suggestions are only minor. Also could you add the testing of Parseval's identity with the amplitude and phase being set to True? There's already a test function (test_parseval) that checks for it using xrft.power_spectrum so you could probably add it there.

xrft/xrft.py Outdated Show resolved Hide resolved
xrft/xrft.py Outdated Show resolved Hide resolved
xrft/xrft.py Outdated Show resolved Hide resolved
xrft/xrft.py Outdated Show resolved Hide resolved
@roxyboy
Copy link
Member

roxyboy commented Dec 11, 2020

It could be interesting to let the user choose between fft or dft call in other functions such as power_spectrum, cross_spectrum, ...

I agree with this. Could we just add the true_phase and true_amplitude flags to power_spectrum and cross_spectrum by also keeping the default as False?

xrft/tests/test_xrft.py Show resolved Hide resolved
@lanougue
Copy link
Contributor Author

Hi @roxyboy ,
Thanks for your review. I think I did all the changes you have suggested.
I recreated two very small tests (1D and 2D) for the parseval identity checking with dft and true flags.
I did not use existing test_parseval function because this later tests several other things at the same time (windowing, ...)

Adding true_phase and true_amplitude flags in power_spectrum and cross_spectrum is feasible (fft call should be transformed into dft call). However, output phase and amplitude will of course be modified and we should evaluate the impact of these modifications in these functions (which is probably also desirable too).

xrft/xrft.py Show resolved Hide resolved
@lanougue
Copy link
Contributor Author

lanougue commented Dec 14, 2020

I am simplifying and adding true_flags possibilities to (co, cross and phase of) spectrum.

I have remaing questions:

  • Why cross_spectrum is returning only the real part ? I think that the imaginary part is as important as the real one. Shouln't we return the complex spectrum ?

  • In computing cross_phase, I used xr.ufuncs.angle() to compute the phase. I think it automattically handle all dask and chunk stuff.

  • What the return of power_spectrum with density flags set to False means ? Shouldn't we use the same syntax of "scipy.signal.periodogram" with a scaling flag set to "density" or "spectrum" ?

  • if cross_spectrum would return the complex value, the cross-phase would also only be xr.ufuncs.angle(cross_spectrum())

(@roxyboy : sorry for the failing checks, my local pytest is not working for now)

@lanougue
Copy link
Contributor Author

Even if there is no need of compatibility of idft with current xrft version, I just restore idft default behaviour with true_flags to False to be consistent with dft. (I set them to true in the commit before). I changed the message of the warning to be more explicit

xrft/xrft.py Outdated Show resolved Hide resolved
@roxyboy
Copy link
Member

roxyboy commented Dec 15, 2020

I think this is ready to merge. In hindsight, we probably should've split up this massive PR into multiple short PRs but "c'est la vie". @rabernat Can you give this PR a final approval?

@roxyboy
Copy link
Member

roxyboy commented Dec 15, 2020

Thanks for your patience and all the work @lanougue !

@lanougue
Copy link
Contributor Author

Thank you too for reviewing all these modifications !

A last thing. In cross_phase, I set true_phase=False as default to keep the default behaviour across version
In beleive that in future we should probably change it to true because it takes into account the coordinates to evaluate the correct cross_phase.
I think we should had a warning here too.

@roxyboy
Copy link
Member

roxyboy commented Dec 15, 2020

A last thing. In cross_phase, I set true_phase=False as default to keep the default behaviour across version
In beleive that in future we should probably change it to true because it takes into account the coordinates to evaluate the correct cross_phase.
I think we should had a warning here too.

I agree, a warning would be nice.

@lanougue
Copy link
Contributor Author

I added warnings for both cross_spectrum and cross_phase. All true_flags default behaviour is now only supported by dft and idft functions ensuring compatibility across all functions which depend on them.
Code looks more clear to me.

xrft/xrft.py Outdated Show resolved Hide resolved
xrft/xrft.py Outdated Show resolved Hide resolved
@lanougue
Copy link
Contributor Author

Loooong PR ! :) Happy that it reaches the end

xrft/xrft.py Outdated
)
warnings.warn(msg, FutureWarning)
if scaling != "density":
if "scaling" in kwargs:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem with hiding scaling in kwargs in my opinion is that the user will have a hard time knowing that the flag even exists as it's not in xrft.dft...

Copy link
Member

@roxyboy roxyboy Dec 15, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we're going to relegate scaling to kwargs, we need to add more explanation to the docstring of the function to explain what it takes as inputs in addition to kwargs accepted in xrft.dft

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I agree. It is in the description of the argument list, but I do not know another way.
Do we have another way to check if the user set an argument explicitely or if it is the default value ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Has this issue been resolved?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think so. "scaling" now appears in argument list as expected. If the user keep using "density", it will work but will raise a deprecation warning

xrft/xrft.py Show resolved Hide resolved
xrft/xrft.py Outdated Show resolved Hide resolved
xrft/tests/test_xrft.py Outdated Show resolved Hide resolved
@lanougue
Copy link
Contributor Author

lanougue commented Jan 5, 2021

Happy new year 2021. May xrft have a long life !
@roxyboy , @rabernat . What about this PR ? followed by a new release ?

@rabernat
Copy link
Collaborator

Sorry @lanougue we have neglected your PR for so long. From my point of view, everything looks good to merge.

@roxyboy do you agree?

@roxyboy
Copy link
Member

roxyboy commented Jan 28, 2021

@rabernat Yes, I agree with merging this PR :)

@rabernat rabernat merged commit 8447e68 into xgcm:master Jan 28, 2021
This was referenced Jan 28, 2021
Closed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants