derivative using FFT properties

Hi, I am trying to make a calculation by taking the first deriviative without success, I've tried the next couple of methods:

    %%%the data is a seismic shot gather
    %%%size(data_TX,1)=1001 (number of time samples)
    %%%siz (data_TX,2)=256     (number of increments)
        data_FK=fft2(data_TX);
    % setting the the time and spacial freq vectors
        fq=linspace(0,1/dx,size(data_TX,1)-1/2*dt);
        kx=linspace(0,1/dy,size(data_TX,2)-1/2*dx);
    %1 method
        [FQ,KX]=meshgrid(ifftshift(fq),ifftshift(kx));
        Dt_data_FK=2i*pi*FQ.*data_FK;
        Dx_data_FK=2i*pi*KX.*data_FK;
    %in this case I need to transpose the matrices
        FQ=FQ';
        KX=KX';
        Dt_data_FK=2i*pi*FQ.*data_FK;
        Dx_data_FK=2i*pi*KX.*data_FK;
        Dt_data_TX=ifft2(Dt_data_FK); 
        Dx_data_TX=ifft2(Dx_data_FK); %I get in this case a complex data set ?

    %2nd method
        data_FK=fftshift(fft2(data_TX));
        Dt_data_FK=2*pi*1i*diag(fq)*data_FK;
        Dx_data_FK=2*pi*1i*data_FK*diag(kx);
        Dt_data_TX=ifft2(ifftshift(Dt_data_FK),'symmetric');
        Dx_data_TX=ifft2(ifftshift(Dt_data_TX),'symmetric');

note: I need to make the calculation for a cube CUBE(k,j,i), k=1001 time samples, j=256 recievers,i=256 shots.

I tried at first to use fftn but all ways produced me with wrong results. the theory is right but i'm not sure exactly how to take the derivative.

[Jonathan Ron]

I use ifft2, sorry for the mistake

Answer 1

[Dr. Seis]

Based on the website below:

http://www.cl.cam.ac.uk/teaching/2000/ContMaths/JGD-notes/node11.html

I believe you would do something similar to the example below:

Nx = 64; % Number of samples collected along first dimension
Ny = 32; % Number of samples collected along second dimension
dx = 1;  % x increment (i.e., Spacing between each column)
dy = .1; % y increment (i.e., Spacing between each row)
x = 0 : dx : (Nx-1)*dx; % x-distance
y = 0 : dy : (Ny-1)*dy; % y-distance
data_spacedomain = randn(Ny,Nx); % random 2D matrix
Nyq_kx = 1/(2*dx); % Nyquist of data in first dimension
Nyq_ky = 1/(2*dy); % Nyquist of data in second dimension
dkx = 1/(Nx*dx);   % x-Wavenumber increment
dky = 1/(Ny*dy);   % y-Wavenumber increment
kx = -Nyq_kx : dkx : Nyq_kx-dkx; % x-wavenumber
ky = -Nyq_ky : dky : Nyq_ky-dky; % y-wavenumber
data_wavenumberdomain = zeros(size(data_spacedomain)); % initialize data
% Compute 2D Discrete Fourier Transform
for i1 = 1 : Nx
for j1 = 1 : Ny
for i2 = 1 : Nx
for j2 = 1 : Ny
data_wavenumberdomain(j1,i1) = data_wavenumberdomain(j1,i1) + ...
(2i*pi*ky(j1))*(2i*pi*kx(i1))* ...
(data_spacedomain(j2,i2)*exp(-1i*(2*pi)*(kx(i1)*x(i2)+ky(j1)*y(j2))));
end
end
end
end

differentiated_data = ifft2(ifftshift(data_wavenumberdomain));

I have a program for performing this on a 1D dataset that you might leverage (in your own way) to produce the same effect on 2D datasets - see my answer in the post below:

http://www.mathworks.com/matlabcentral/answers/21700-finding-the-velocity-from-displacement

[EDIT 6/22]

Practically speaking, if you were going to compute the derivative of a 2D image, where your dimensions are either space-space or time-time, then you would need to use fft2, ifftshift, and meshgrid:

Nx = 64; % Number of samples collected along first dimension
Ny = 32; % Number of samples collected along second dimension
dx = 1;  % x increment (i.e., Spacing between each column)
dy = .1; % y increment (i.e., Spacing between each row)
x = 0 : dx : (Nx-1)*dx; % x-distance
y = 0 : dy : (Ny-1)*dy; % y-distance
data_spacedomain = randn(Ny,Nx); % random 2D matrix
Nyq_kx = 1/(2*dx); % Nyquist of data in first dimension
Nyq_ky = 1/(2*dy); % Nyquist of data in second dimension
dkx = 1/(Nx*dx);   % x-Wavenumber increment
dky = 1/(Ny*dy);   % y-Wavenumber increment
kx = -Nyq_kx : dkx : Nyq_kx-dkx; % x-wavenumber
ky = -Nyq_ky : dky : Nyq_ky-dky; % y-wavenumber

% Compute 2D FFT
data_wavenumberdomain = fft2(data_spacedomain);

% Compute grid of wavenumbers
[KX, KY] = meshgrid(ifftshift(kx),ifftshift(ky));

% Compute 2D derivative
data_wavenumberdomain_differentiated = (2i*pi)^2*KX.*KY.*data_wavenumberdomain; 

% Convert back to space domain
data_spacedomain_differentiated = ifft2(data_wavenumberdomain_differentiated );

If you only want to take the partial derivative of the the "image" with respect to one of your dimensions, then all you would do is:

d_MAT_x = 2i*pi*KX.*data_wavenumberdomain;
d_MAT_y = 2i*pi*KY.*data_wavenumberdomain;

Then just transform back into time domain using ifft2

[EDIT 6/27]

1. "Nx" and "Ny" are the number of columns and number of rows, respectively. Thus, you should have:

[KX,FQ]=meshgrid(ifftshift(kx),ifftshift(fq));

2. Your "fq" and "kx" may be defined wrong (I said in an earlier comment that if the number of rows or columns is an odd number then we need to do something different). Also, it seems strange that you would define "fq" in terms of "dx" and "dt"... and "kx" in terms of "dy" and "dx". Lets assume "dt" is the time increment and "dx" is the spatial increment... "Nt" is the number of time samples and "Nx" is the number of spatial samples, then:

Nyq_fq = 1/(2*dt);
Nyq_kx = 1/(2*dx);
dfq = 1/(Nt*dt);
dkx = 1/(Nx*dx);

if mod(Nt,2) == 0 % Nt is even
     fq = -Nyq_fq : dfq : Nyq_fq-dfq;
else            % Nt is odd
     fq = [sort(-1*(dfq:dfq:Nyq_fq)) (0:dfq:Nyq_fq)];
end

if mod(Nx,2) == 0 % Nx is even
     kx = -Nyq_kx : dkx : Nyq_kx-dkx;
else            % Nx is odd
     kx = [sort(-1*(dkx:dkx:Nyq_kx)) (0:dkx:Nyq_kx)];
end

Does this help solve things?

[Jonathan Ron]

Hi Elige,

I'm not looking to build a new fft routine, the matlab one suits me and the 1st part I can obtain by just doing fft2(data_spacedomain).

according to fourier properties multiplication by -2pi will give the deriviative, when I'm doing it I get a worng answer and I don't realize what I'm doing wrong.

let's say I have an Image and I transfrom from tx domain to fk domain,

DATA_FK=fft2(DATA_TX);

the data on DATA_FK is fine, now according to fourier properties, the deriviative is multiplication with -2*pi*1i*k_vec or -2pi*1i*f_vec depands on the domain.

I need for my calculation to get back a matrix witht the time deriviative (d/dt) and another matrix with the spacial deriviative (d/dx).

I don't really understand what the command differentiated_data = fft2(ifftshift(data_wavenumberdomain))*dkx*dky; does

[Dr. Seis]

The "new" fft routine was just to kill two birds with one stone... it demonstrates what the 2D DFT is doing, while also computing the derivative - see the part inside the 4 "for" loops that has "(2i*pi*ky(j1))*(2i*pi*kx(i1))" ? That is where the "differentiation" is taking place. What I am doing is taking the derivative of a 2D "image" whereas it looks like you are trying to generate two images that represent the partial derivative in each dimension. And it should be positive 2*pi*1i*frequencies (not negative).

I made some edits above.

[Dr. Seis]

The above should work provided the number of rows and columns are even numbers. If any of them are odd, then we need to redefine your wavenumber or frequency values.

[Jonathan Ron]

Hi Elige,

I still have problems with the solution

I think it is got to do with this part

[KX, KY] = meshgrid(ifftshift(kx),ifftshift(ky))

as u seen above, I multiply it with the diag(kx) (which does not work)

when I do it your way I need to transpose the mesh

KX=KX'

KY=KY'

and maybe it causes the problem....I still cannot get a correct answer.........

[Dr. Seis]

In my example, "Nx" (i.e., the number of elements in "x" and "kx") corresponds to the number of columns in my 2D "image" while "Ny" (i.e., the number of elements in "y" and "ky") corresponds to the number of rows.

What is "Nx" and "Ny" in your case?

[Dr. Seis]

I have edited my answer above... there appears to be a couple of misunderstandings.

[Jonathan Ron]

Hi Elige,

the fq and kx are fine, I don't have a problem with the odd number the result I'm getting is good now, using the meshgrid, except I have a couple of problems I can't really figure out

1. matlab defines the fft with -2*pi*i and not 2*pi*i. In order to get good results I need in both ways to multiply the derivative in (-1) after the ifft.

2. the partial derivative

     Dx_data_TX=ifft2(Dx_data_FK);

gives me a complex data sets,usually something like (anynumber)+ 0.0000001i so when I take the real part it's ok, but still it is not suppose to be like that.

thanks for the help, if u have any ideas about these two problems let me know

[Dr. Seis]

1. Yes, the forward Fourier transform is "-2*pi*i". However, you are actually performing the derivative on the inverse Fourier transform:

d/dt( f(t) ) = d/dt( ifft(F(w) )

See the link below (about 4/5ths the way down):

http://fourier.eng.hmc.edu/e101/lectures/handout3/node2.html

So we should be using +2*pi*i. You should not need to multiply anything after the ifft... you really should consider using the way I define "fq" and "kx" because it looks like you are only defining the frequencies/wavenumbers associated with half the spectrum! You need to define the full spectrum which includes both negative and positive frequencies/wavenumbers.

2. This may be due to some precision issues to where the real part of the spectrum is not perfectly symmetrical (i.e., F(w) = F(-w) ), or the imaginary part of the spectrum is not perfectly anti-symmetrical (i.e., F(w) = -F(-w) ). But multiplying the frequency/wavenumber amplitudes by the wrong frequencies/wavenumbers will certainly not help.

Answer 2

[Francesco]

Hello.

When I use the script shown after "EDIT 6/22" by Elige Grant, I seem to get complex numbers in the "data_spacedomain_differentiated" matrix. Shouldn't the signal just be real? Should I just extract the real component or am is that script somehow incomplete??

Cheers,

F