> When I use GRDFFT to compute the power spectrum in the x and y directions:
>
> grdfft temp_topo.grd -M -Exw > fftx.out
> grdfft temp_topo.grd -M -Eyw > ffty.out
>
> 1. I still get the output in frequency not wavelength (not a big problem since
> I can compute the wavelength from freq, just wondering why the flag doesn't
> work). Using -Ewx or -Ewy both return the radial power spectrum. Anyone else
> have this problem?
I guess I do. I checked grdfft.c and found that "give_wavelength" is
BOOLEAN in the main program but int in the do_spectrum subroutine. Changing
it to BOOLEAN didn't make a difference.
give_wavelength = FALSE;
j = 2;
while(argv[i][j]) {
if (argv[i][2] == 'x' || argv[i][2] == 'X')
par[par_count] = 1.0;
else if (argv[i][2] == 'y' || argv[i][2] == 'Y')
par[par_count] = -1.0;
else if (argv[i][j] == 'w' || argv[i][j] == 'W')
give_wavelength = TRUE;
j++;
}
The above fragment of code looks like it never sees the second else if, because
for -Exw or -Eyw, the check for x or y is true. The second else if should just
be an "if". It won't be true for -Ex or -Ey, but would for -Exw or -Eyw.
Changing it worked for me.
> 2. A bigger problem: the wavelength spacing in x and y are different (despite
> the fact that temp_topo.grd is square in x and y = 240x240 units)
> might this have something to do with the fact that temp_topo.grd is centered
> at 40N ???
Probably. You could use grdproject to convert it to a rectangular
coordinate system.
> 3. How is "1 standard deviation in power[f]" computed and how might this
> information be useful?
1-D case: n/2 power spectral estimates at frequencies i/L, where
i=1,n/2 (length L=n*dt). Total ndata, nest = ndata/n independent estimates.
Standard error is then 1/sqrt(nest).
> Also, has anyone modified GRDFFT to compute the power spectrum in a given
> azimuth, not just x and y. Or, alternatively, modified GRDFFT so that it
> returns the transform coefficents
Not that I know of.
Ben Horner-Johnson
ben@earth.nwu.edu