search

brmather / pycurious

Python package for computing the Curie depth from the magnetic anomaly
https://brmather.github.io/pycurious/
GNU Lesser General Public License v3.0
38 stars 17 forks source link

radial power spectrum #3

Closed zhongpenggeo closed 5 years ago

zhongpenggeo commented 5 years ago

hello, I have a question about your code to compute radially averaged power spectrum. In grid/curie.py, you compute power spectrum in function _radialspectrum

FT = np.abs(np.fft.fft2(data*vtaper))

But as I known, it shoud be FT = np.abs(np.fft.fft2(data*vtaper))**2, the square of absolute of fft, isn't it?

brmather commented 5 years ago

Thanks for your query. Could you provide a reference for that? To my knowledge the square is used simply to make all values in the power spectrum positive, which is accomplished with np.abs.

Could you provide more context, please?

zhongpenggeo commented 5 years ago

Thanks for your query. Could you provide a reference for that? To my knowledge the square is used simply to make all values in the power spectrum positive, which is accomplished with np.abs.

Could you provide more context, please?

https://github.com/fatiando/fatiando/blob/ac2afbcb2d99b18f145cc1ed40075beb5f92dd5a/fatiando/gravmag/transform.py#L542 I find Power Density Spectra in Fatinado is square. And you userr = 2.0*np.log(FT[mask]) to compute log of Power Density Spectra in function 'radial_spectrum' in pycurious-src/pycurious/grid/curie.py , which is the same as log((FT)**2)

brmather commented 5 years ago

Sure, so it follows then there is no difference in the radial power spectrum. These statements return the same results:

  1. 2*np.log(np.abs(np.fft.fft2(data*vtaper)))
  2. np.log(np.abs(np.fft.fft2(data*vtaper))**2)

Statement 1 is slightly more computationally friendly.