How to define the spatial frequency in Matlab?

Hello all,

I am working on the simulation of 1-D nonlinear elastic waves using Iteration contrast source technique. In which, I have some questions on a correct definition for spatial frequency and wavenumber. The spatial differential operator in the nonlinear force term (dx^2 and dx) will become (-k^2 and -j*k) in the frequency domain (j-imaginary number). I understand that wrong definition of k vector has a huge impact on the results because of the multiplication of square terms (-k^2) to u (displacement). I have shown my code which has different definitions of k (where M-spatial sample points). Also, whether the spatial frequency of a longitudinal wave (k_x) and shear wave (k_y and k_z) are same or different? Please share your knowledge. Many thanks in advance.

 %%Spatial frequency
% dk=1/(dx);                      % Wavenumber increment
% kp=([(0:M/2-1) (-M/2:-1)]./M) * (dk); % == k_x      %with/without 2*pi
% ks=([(0:M/2-1) (-M/2:-1)]./M) * (dk); % == k_x
% kp=(0:M-1)*(dk);
% ks=(0:M-1)*(dk);
% Zero padded
% dk=1/(dx);                      % Wavenumber increment
% kp=([(0:(2*M)/2-1) (-(2*M)/2:-1)]./(M)) * (dk); % == k_x%with/without 2*pi
% ks=([(0:(2*M)/2-1) (-(2*M)/2:-1)]./(M)) * (dk); % == k_x%with/without 2*pi

[David Goodmanson]

Hi Elango,

Will you be using Matlab fft to do this?

[Elango Selvam]

Yes. I use Matlab fft and ifft.

Answer 1

[David Goodmanson]

Hi Elango,

I want to get the terminology well defined, since in space there is a bit of ambiguity as compared to time, where f and w are not in question.

For the wave number I assume you mean k = 2*pi/lambda. A wave is exp(j*k*x).

For spacial frequency I assume you mean 1/lambda. (No well established symbol for this? I will call it inv**). A wave is exp(2*pi*j*inv*x).

The rules for fft and ifft are

let dx be the spacing of the x array, similarly for k and dk, inv and dinv.
let M be the number of points in the fft or ifft.
then

dx*dk = 2*pi/M
dx*dinv = 1/M

so it looks like with the k defined above, your dk needs a factor of 2*pi/M.

These arrays cover f(x) going to the transform g(k) or g(inv).

Beware of using fft(f,N) where N is usually obtained from nextpow2(M). This can be done but it changes the M used in the formulas above, so if you had already set up the k array you would have to go back and change it.

** spectroscopists at least have a name for it, inverse cm, but I am not implying specific units here.

[Elango Selvam]

Thank you for your answer. The factor 2*pi plays role when I multiply my angular wavenumber (which is 2*pi*inv). But k_p and k_s in my code denotes spatial frequency in x and y,z directions where the total wavenumber is k=sqrt(kp^2+ks^2). So I am trying to define array for kp and ks alone in this case, which follows:

dk=1/(M*dx);                      % spatial frequency increment
kp=([(0:M/2-1) (-M/2:-1)]) * (dk); % == k_x      %with/without 2*pi
ks=([(0:M/2-1) (-M/2:-1)]) * (dk); % == k_y & k_z

Is that right? Also, I have another question here, by physics will wave number (or spatial frequencies) for a pressure wave (k_p) and a shear wave (k_s) will be equal? Because my array for k_p and k_s are same by above definition. But both waves have different speeds. And the difference between k_p and k_s will be taken care within total wave number k=sqrt(kp^2+ks^2), is that correct?.

[David Goodmanson]

Hi Elango,

You should verify this in some references on seismology, but since the p and s are independent waves believe there are two independent k's, with

ks = sqrt(ks_x^2 + ks_y^2 + ks_z^2)
kp = sqrt(kp_x^2 + kp_y^2 + kp_z^2)

and as you say they have different speeds, so for circular frequency w,

ks = w/vs
kp = w/vp

[Elango Selvam]

Thanks for your explanation. I am considering the 1-D wave propagation (x-direction) in 3-D system (an Aluminium block), so wavenumbers for p and s- wave becomes,

 ks = sqrt(ks_x^2 + 0 + 0)
kp = sqrt(kp_x^2 + 0 + 0)

The pressure and shear waves both travel in the x-direction only, at different speeds, so kp_y and kp_z, similarly ky_y and ky_z become zero. So, I have two waves Pressure and shear wave, traveling in x-direction, correspondingly I get two wave numbers (kp=kp_x and ks=ks_x) for each of them. Until now everything seems ideal to me. Well now I should have two different spatial (frequency) sampling distance (dk) that corresponds to both waves, i.e., dk_p, dk_s related to dx_p and dx_s. So I understand that I need to define two x-axes, for each waves. Please correct me if I am wrong.

[David Goodmanson]

Isn't there just one kx for each wave, ks_x for the shear wave and ks_x for the primary wave? Of course the displacement direction of the rock is different for the two of them, x direction for p and some linear combination of y and z direction for s (for an isotropic material).

[Elango Selvam]

Yes, there is only one kx for each waves, because wave propagates only in x-direction. The component of the pressure wave oscillates in x-direction whereas the components of Shear wave oscillates in y- & z- direction, and both the waves travels in x-direction only.