Hi,all
grdfft is a good tool to do the processing and computing of grid files. I used it very frequent with no technical problem. But I also have some questions similar to the topic of references for grdfft.
I test the spectrum computation using the grdfft and spectrum1D.
The result from spectrum1d is the power spectrum density and it follows specific rules as Parseval’s theorem. Generally speaking, the variance of the input data equals the integral of the spectrum density.
But, what is spectrum output from the grdfft? Is it the power spectrum density or not?
I simulated a simple grid file and tested the two different tools. The spectrum1D follows the above theorem correctly. The grdfft dose not.
Also, I am not sure what is relation between the spectrums from x,y and the radial directions. Are they follow the formula:r^2=x^2+y^2 or not? From my test, it dose not.
Here is the simulated figure
Here is the code:
#!/usr/bin/env bash
ps=testgrd.ps
# gmt grdmath -R-50/50/-50/50 -I0.2 X Y ADD SIN = x.nc
gmt grdmath -R-50/50/-50/50 -I0.2 X SIN = x.nc
gmt makecpt -Crainbow -T-1/1/0.1 > t.cpt
gmt grdimage x.nc -R-30/30/-30/30 -JX3i -Ct.cpt -Bxa -Bya -BWSne -P -K -Y7i > $ps
gmt grd2xyz x.nc | awk '{if ($2 == 0) print $0}' > x.txt
gmt psxy -R -J -O -K x.txt -Sc0.03i -Gblack >> $ps
# Do grdfft
gmt grdfft x.nc -Er+w+n -Gfft.txt
gmt grdfft x.nc -Ex+w+n -Gfftx.txt
gmt grdfft x.nc -Ey+w+n -Gffty.txt
# Plot specrum result from grdfft
gmt psbasemap -R1/20/1e-14/1e-2 -JX-2.8il/3il -Bxa1pf2g3+l"Wavelength" -Bya1pg+l"power spectrum" -BwNsE -K -X3.2i -O >> $ps
gmt psxy -R fft.txt -J -W0.5p,red, -O -K >> $ps
gmt psxy -R fftx.txt -J -W0.5p,blue, -O -K >> $ps
gmt psxy -R ffty.txt -J -W0.5p,green, -O -K >> $ps
gmt psbasemap -R0.05/1e0/1e-14/1e-2 -JX2.8il/3il -BSE -Bxa1p+l"wavenumber/cycle" -O -K>> $ps
gmt pslegend -D+w1.2i+jBL+o0.1i/0.1i -R -J -O -K -F+p1p+glightgray --FONT_ANNOT_PRIMARY=10p,1,black << EOF >> $ps
S 0.2i - 0.5i - 1p,red, 0.5i r
S 0.2i - 0.5i - 1p,blue 0.5i x
S 0.2i - 0.5i - 1p,green, 0.5i y
EOF
# Do spectrum1D
gmt psxy -R-30/30/-1/1 -JX3i/3i x.txt -i0,2 -Bxaf -Byaf -BWSne -O -K -X-3.2i -Y-4.5i >> $ps
gmt psxy -R -J x.txt -i0,2 -Sc0.03i -Gred -O -K >> $ps
awk '{print $1,$3}' x.txt | gmt spectrum1d -S256 -W --GMT_FFT=brenner -N -i1 -D0.2 > fft_spec.txt
gmt psbasemap -R1/20/1e-12/1e2 -JX-2.8il/3il -Bxa1pf2g3+l"Wavelength" -Bya1pg+l"power spectrum" -BwNsE -K -X3.2i -O >> $ps
gmt psxy -R fft_spec.txt -J -W0.5p,red, -O -K >> $ps
gmt psbasemap -R0.05/1e0/1e-12/1e2 -JX2.8il/3il -BSE -Bxa1p+l"wavenumber/cycle" -O >> $ps
awk '{print $3}' x.txt| gmt gmtmath STDIN -Sl STD SQR =
gmt psconvert $ps -P -Tg -A