[Scilab-users] FFT - upward continuation of a 2D grid

classic Classic list List threaded Threaded
7 messages Options
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] FFT - upward continuation of a 2D grid

Hello all,

A quick query regarding fft and wavenumber calculations. I am looking to do
an upward contuation of a gravity data grid using fft and the basic
equation:

fft_Grav_up = fft2(gravity_input) * e^(-kz) // think it also needs fftshif
to centre Zero?

... Upward_continued_gravity=ifft(fft_out_Grav_up)

Where z = upward continuation distance and k = wavenumber

k =sqrt(kx^2 + ky^2)

Is there an easy way to get the wavenumber from the input?

Before doing the fft, the grid would need padding on all edges, so I think
this works unless there is a better method:

dimx=size(xt)
dimy=size(yt)

// For FFT pad X,Y to 2*ncols, 2*nrows (large enough to avoid boundary
effects)

addedRows = dimy(1); addedCols = dimx(1);

tmp = [Boug_corr; Boug_corr($, :).*.ones(addedRows, 1)];
tmp = [tmp  tmp(:, $).*.ones(1,addedCols)];
tmp = [repmat(tmp(1, :), addedRows, 1) ; tmp];
Boug_corr_padded = [repmat(tmp(:,1), 1, addedCols) tmp]

After running fft, the grid would then be clipped back to the correct size.

If there is a simple code for doing this great!

Thanks
Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: FFT - upward continuation of a 2D grid

Hi Lester,

If I understood properly, for an even number of spatial input grid points along Nx and Ny, with regular spacing dx and dy, I think k should be computed from:

kx = dkx*(-Nx/2+1:Nx/2);
ky = dky*(-Ny/2+1:Ny/2);

with:
   dkx = 1/Nx/dx;
   dky = 1/Ny/dy;

Regards,
Rafael

-----Original Message-----
From: users <[hidden email]> On Behalf Of arctica1963
Sent: Thursday, July 26, 2018 3:54 PM
To: [hidden email]
Subject: [Scilab-users] FFT - upward continuation of a 2D grid

Hello all,

A quick query regarding fft and wavenumber calculations. I am looking to do an upward contuation of a gravity data grid using fft and the basic
equation:

fft_Grav_up = fft2(gravity_input) * e^(-kz) // think it also needs fftshif to centre Zero?

... Upward_continued_gravity=ifft(fft_out_Grav_up)

Where z = upward continuation distance and k = wavenumber

k =sqrt(kx^2 + ky^2)

Is there an easy way to get the wavenumber from the input?

Before doing the fft, the grid would need padding on all edges, so I think this works unless there is a better method:

dimx=size(xt)
dimy=size(yt)

// For FFT pad X,Y to 2*ncols, 2*nrows (large enough to avoid boundary
effects)

addedRows = dimy(1); addedCols = dimx(1);

tmp = [Boug_corr; Boug_corr($, :).*.ones(addedRows, 1)]; tmp = [tmp  tmp(:, $).*.ones(1,addedCols)]; tmp = [repmat(tmp(1, :), addedRows, 1) ; tmp]; Boug_corr_padded = [repmat(tmp(:,1), 1, addedCols) tmp]

After running fft, the grid would then be clipped back to the correct size.

If there is a simple code for doing this great!

Thanks
Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: FFT - upward continuation of a 2D grid

You are welcome

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: FFT - upward continuation of a 2D grid

In reply to this post by Rafael Guerra
Hi Rafael,

Thanks for the pointers. Still having some issues getting it to work on real
data. Upward continuation, should act like a low-pass and give a smoother
view of the data.

dimx=size(xt);
dimy=size(yt);

nx=dimx(1);
ny=dimy(1);

dxy = 0.016666666;
dkx=1/(nx*dxy);
dky=1/(ny*dxy);

kx=dkx*(-(nx/2)+1:nx/2);
ky=dky*(-(ny/2)+1:ny/2);

[KX,KY]=meshgrid(kx,ky);
k=sqrt(KX.^2 + KY.^2);

uc=4000; // continue up 4000m

//fft_Boug_corr=fft2(Boug_corr);
//shift_f=fftshift(fft_Boug_corr);
//fft_clean=clean(shift_f);

fft_clean=clean(fftshift(fft2(Boug_corr))); // merge commands

Fup=fft_clean.*exp(-k*uc);

F_real=real(ifft(Fup));

It does generate output, but certainly not as expected for upward
continuation of data. Not sure what I'm missing, assuming I have done the
wavenumber thing correct. Normally I would do this via GMT, but want to get
it going under Scilab for convenience.

Clearly, before doing a proper fft, we need to pad the data to expand the
dimensions over which the fft is done to avoid edge effects. Just running it
with the grid dimensions should still produce a recognisable result.

Matbe testing with a synthetic grid would help.

Cheers
Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: FFT - upward continuation of a 2D grid

Hi Lester,

The spatial resolution of 0.016 m seems unrealistically too-finely sampled... Could you confirm the units?
Also, it is possible that your upward continuation requires an angular wavenumber:   exp(-2*%pi*k*uc) ?

Regards,
Rafael

-----Original Message-----
From: users <[hidden email]> On Behalf Of arctica1963
Sent: Saturday, July 28, 2018 2:23 PM
To: [hidden email]
Subject: Re: [Scilab-users] FFT - upward continuation of a 2D grid

Hi Rafael,

Thanks for the pointers. Still having some issues getting it to work on real data. Upward continuation, should act like a low-pass and give a smoother view of the data.

dimx=size(xt);
dimy=size(yt);

nx=dimx(1);
ny=dimy(1);

dxy = 0.016666666;
dkx=1/(nx*dxy);
dky=1/(ny*dxy);

kx=dkx*(-(nx/2)+1:nx/2);
ky=dky*(-(ny/2)+1:ny/2);

[KX,KY]=meshgrid(kx,ky);
k=sqrt(KX.^2 + KY.^2);

uc=4000; // continue up 4000m

//fft_Boug_corr=fft2(Boug_corr);
//shift_f=fftshift(fft_Boug_corr);
//fft_clean=clean(shift_f);

fft_clean=clean(fftshift(fft2(Boug_corr))); // merge commands

Fup=fft_clean.*exp(-k*uc);

F_real=real(ifft(Fup));

It does generate output, but certainly not as expected for upward continuation of data. Not sure what I'm missing, assuming I have done the wavenumber thing correct. Normally I would do this via GMT, but want to get it going under Scilab for convenience.

Clearly, before doing a proper fft, we need to pad the data to expand the dimensions over which the fft is done to avoid edge effects. Just running it with the grid dimensions should still produce a recognisable result.

Matbe testing with a synthetic grid would help.

Cheers
Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: FFT - upward continuation of a 2D grid

One additional remark.
Fftshift might have been misused. Try instead:

   fft_clean=clean(fft2(Boug_corr)); // merge commands
   k= fftshift(k);
   Fup=fft_clean.*exp(-k*uc);

Regards,
Rafael

-----Original Message-----
From: users <[hidden email]> On Behalf Of Rafael Guerra
Sent: Saturday, July 28, 2018 2:44 PM
To: Users mailing list for Scilab <[hidden email]>
Subject: Re: [Scilab-users] FFT - upward continuation of a 2D grid

Hi Lester,

The spatial resolution of 0.016 m seems unrealistically too-finely sampled... Could you confirm the units?
Also, it is possible that your upward continuation requires an angular wavenumber:   exp(-2*%pi*k*uc) ?

Regards,
Rafael

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: FFT - upward continuation of a 2D grid

Hi Rafael,

Problem solved with your guidance! The cellsize was in decimal degrees, so
converted to linear units 0.016666 degrees = 1850m approx.

//dxy = 0.016666666; decimal degrees - 1 arc minute
dxy = 1850;
dkx=1/(nx*dxy);
dky=1/(ny*dxy);

...

fft_clean=clean(fft2(Boug_corr));
k=fftshift(k);

Fup=fft_clean.*exp(-2*%pi*k*uc); // needed to be angular k

F_real=real(ifft(Fup));

fftshift was applied wrong, so thanks for spotting that!

All I need to do to finish this is pad the grid, run the fft, and clip it
back to the original size. A rough check of the output compared to that from
GMT is minor

Thanks again!



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users