I use ifft2, sorry for the mistake

# derivative using FFT properties

53 views (last 30 days)

Show older comments

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.

### Accepted Answer

Dr. Seis
on 21 Jun 2012

Based on the website below:

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:

[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?

##### 8 Comments

Dr. Seis
on 28 Jun 2012

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):

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.

### More Answers (1)

Francesco
on 18 Jan 2013

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

##### 0 Comments

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!