# Finding the mode of propagation of a progressive wave and the velocity of that mode

Injecting a sinusoid say, $$Z = A \sin(k_x X + k_y Y - \omega t)$$, $$K = 2 \pi f/c$$, for example, assume $$f = 10$$ Hz or $$50$$ Hz, and $$c = 50$$ m/s, now if I draw the 2D K-space plot with the following script

n=40
f = 50  # Temporal frequency in Hz
c = 50  # wave speed 50  m/s
k = 2.0 * np.pi * f / c  # wavenumber
kx = np.linspace(-k, k, n * 10)  # space vector
ky = np.linspace(-k, k, n * 10)  # space vector
X, Y = np.meshgrid(kx,ky)
t= np.arange(0, n, 1)

#2D sine wave
d2_sig = np.sin(-(kx*X+ky*Y)+2.0 * np.pi *f*2)

plt.figure()
c = plt.imshow(d2_sig, cmap='seismic', vmin=d2_sig.min(), vmax=d2_sig.max(),
extent=[kx.min(), kx.max(), ky.min(), ky.max()],
interpolation='nearest', origin='lower')
plt.colorbar(c)
plt.xlabel("Wavenumber, $$K_x$$ [rad/m]", fontsize=18)
plt.ylabel("Wavenumber,$$K_y$$ [rad/m]", fontsize=18)
plt.title(f'Progressive wave', weight="bold")


The plot I get:

I'm worried to find the modes of propagation of the wave vector and also the speed of the wave.

How do I interpret the image in figure 03 representing spatial spectrum of a Wave, where the spatial spectrum is expressed by the equation

Here, $$C(r_i,r_j,\omega)$$ represents the cross-spectral density between every pair of sensors. $$k=2\pi f/c$$ is the wavenumber $$r_i, r_j$$ be the position of sensors.

N.B. I've pointed out what was going wrong on the above script to generate a 2D spatial plane wave and the answer has provided in the following.

However, I've experienced another problem to inject a planar Sine wave to spatially positioned sensors, the sensor position has known to me and I provide the sensor position as follows:

x = [2.1, 2.1, -0.7, -2.1, -2.1, -0.7, -0.7, 0.6, -5.7, -8.5, -11.4, -7.7, -6.3, -3.5, -2.1, -3.4, 5.4, -5.2, -8.9, -10,
-10, 5.4, 5.4, -0.8, -3.6, -6.2, -6.8, -12.2, -17.1, -19, -18.6, -13.5, 14.8, 14.8]
y = [6.65, 4.15, 3.65, 5.05, 7.25, 8.95, 11.85, 8.95, -2, -0.6, -0.9, 1.25, 2.9, 0.9, -0.1, -1.4, 9.2, 5.2, 4.8, 6.1,
8.9, 13.3, 17.1, 17.9, 13.8, -9.3, -5.2, -3.6, -3.6, -0.9, 3.7, 3.7, -1.8, 5.7]


The sensors position is at irregular intervals and spatial resolution has not been maintained here. What would be the possible approach to inject a plane wave to these sensors structure?

• In 2D K-space, I think $k_x$ and $k_y$ should satisfy the relationship of $k_x^2+k_y^2=k^2$ where $k=2\pi f/c$ is the wave number. In general the solution of wave equation is given by complex exponential function $e^{j(k_xx+k_yy-\omega t)}$, so when $k_x$ or $k_y$ is greater than $k$, the other one is an imaginary number which means along that dimension the wave has an exponential decay. Commented Feb 9, 2022 at 23:21
• Thanks. What is the primary and secondary mode of propagation? Is it seen from the plot? How to distinguish in which direction the wave propagates? Also how to find the speed of the wave? Commented Feb 10, 2022 at 4:10
• What do you mean by mode? Does it mean discretization of $k_x$ and $k_y$? I think these modes exist in waveguides, i.e. boundary conditions are required. There is no discrete modes for a wave propagating in free field. In free field the speed of the wave is equal to $c$ which is 50 m/s in your case, while in waveguides higher order waves have different propagation speed. Commented Feb 10, 2022 at 4:18
• From the figure, how to fix the direction of propagation of the wave along which direction? I mean where is the source and destination of the traveling wave in the K-space ? Commented Feb 10, 2022 at 4:30

In your case the wave number $$k = 2\pi f/c=2\pi$$, $$k_x$$ and $$k_y$$ should satisfy $$k^2 = k_x^2+k_y^2$$ which means they are not always real. But the direction of propagation of wave is easy to obtain given $$k_x$$ and $$k_y$$. The wave vector contains the information of both direction and the wave number: $$\vec{k}=k_x\vec{e}_x + k_y\vec{e}_y$$ where $$\vec{e}_x$$ and $$\vec{e}_y$$ are unit vector along $$x$$ and $$y$$ axis, respectively.

Edit:

@ZaellixA already gives you a visualization of plane wave. Let me give you a animation plot to have a further look.

f = 10;         % frequency
fs = 100;       % sample frequency
Ts = 1/fs;      % sample period
t = 0:Ts:0.2;   % time index
c = 50;         % speed of wave
w = 2*pi*f;     % angular frequency
k = w/c;        % wave number
x = linspace(-10, 10);
y = linspace(-10, 10);
[xx, yy] = meshgrid(x, y);
theta = pi / 6; % direction of propagation
kx = k * cos(theta);
ky = k * sin(theta);

h = figure;
for n = 1:length(t)
p = exp(-1j*(kx * xx + ky * yy - w * t(n))); % plane wave
imagesc(x, y, real(p));
xlabel('$$x /$$m', 'Interpreter', 'latex')
ylabel('$$y /$$m', 'Interpreter', 'latex')
c = colorbar;
c.Label.String = 'Ampitude';
caxis([-1, 1])
drawnow

frame = getframe(h);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if n == 1
imwrite(imind,cm,'planewave.gif','gif', 'Loopcount',inf, 'DelayTime', Ts);
else
imwrite(imind,cm,'planewave.gif','gif','WriteMode','append', 'DelayTime', Ts);
end
end


The above MATLAB code generates an animation

Note that the axis of this figure is spatial axis $$x$$ and $$y$$, and the color axis denotes amplitude. In your picture the axis is $$k_x$$ and $$k_y$$, where every point represents a pair of $$(k_x, k_y)$$ and a corresponding plane wave with a direction of $$\arctan2(k_y, k_x)$$.

• Thanks for your patience. The signal travels in both directions in k-space. There are two options [-3, 0] & [0, 3]. Can I say the signal starts travel from [-3, 0] or [0, 3]? Mathematical point of view, I understand your point, but what does physically meaning? Commented Feb 10, 2022 at 4:51
• @Alan22 No you can't say it. The wave travels in the spatial space not k-space. $\vec{e}_x$ and $\vec{e}_y$ are unit vector in spatial space. Say a plane wave $e^{j(\omega t-\vec{k}\cdot\vec{r})}$ where $\vec{k}$ is wave vector I already defined in the answer, and $\vec{r}=x\vec{e}_x+y\vec{e}_y$ is spatial variable. This wave travels along the direction of $\vec{k}$. Commented Feb 10, 2022 at 5:01
• So in k-space, every point represents a pair of $k_x$ and $k_y$, which yields the wave vector $\vec{k}$ and the direction of wave as well. Commented Feb 10, 2022 at 5:04
• Probably you are right, but still I'm not fully influenced from the physical point of view (sorry for my understanding). Let a 1D signal ' d1_sig = np.sin(-(kx*X)+2.0 * np.pi *f*2)' then it gives vertical lines which represents the wave, traveling to X-direction. Commented Feb 10, 2022 at 5:11
• @Alan22 Yes. For 2D signal, $\sin(x+y)$ gives you a plane wave in the direction of $45^\circ$, where $k_x=k_y=1$, and $k=2\pi f/c=\sqrt{k_x^2+k_y^2}$. Does it make sense? Commented Feb 10, 2022 at 5:36

As already commented so far, for the linear wave equation to be satisfied in any number of dimensions the condition (here shown for three dimensions)

$$k^{2} = k_{x}^{2} + k_{y}^{2} + k_{z}^{2}$$

must hold. Now, consider the complex exponential solution of the linearised wave equation of the form

$$p(x, y, z, t) = A e^{-j \left( k_{x} x + k_{y} y + k_{z} z\right)} e^{-j \omega t}$$

The last term of above equation is the temporal term and will be skipped for the moment. The above solution represents a plane wave travelling in the direction specified by the first exponential term. To see that better consider the exponent (omitting the imaginary unit) to be of the equivalent form

$$\mathbf{k_{x}} + \mathbf{k_{y}} + \mathbf{k_{z}} = k_{x} \mathbf{e}_{x} + k_{y} \mathbf{e}_{y} + k_{z} \mathbf{e}_{z}$$

where with bold are indicated vector values and $$\mathbf{e}_{x}$$, $$\mathbf{e}_{y}$$ and $$\mathbf{e}_{z}$$ represent unit vectors pointing along the specified directions. Thus, their weighted vector addition (weights being the $$k$$ values) forms a vector pointing towards a specific direction, albeit not a unit vector anymore. Nevertheless, the resultant vector points towards the direction of propagation of the plane wave in 3D. Of course, if any of the $$k_{i}$$ values is identically 0 then the propagation happens in only two dimensions, like in your case that happens in the $$x-y$$ plane.

Now, moving to your case specifically. In order to "create" the direction of your plane wave you have to create the components of the resultant unit vector pointing towards the direction of propagation. This is quite easy for the 2D case. In order to create a unit vector pointing towards a specific direction two specific conditions must be satisfied. First, the magnitude of the vector to be 1 and second, its angle with the $$x$$ axis must equal the angle of propagation. An easy choice that satisfies both conditions is

$$\mathbf{e}_{2} = \left[ \cos \left( \theta \right), \sin \left( \theta \right) \right]$$

where the subscript denotes that this is vector in $$\mathbb{R}^{2}$$. So, the resulting equation (for two dimensions this time) becomes

$$p(x, y, t) = A e^{-j \left( k_{x} \cos \left(\theta\right) x + k_{y} \sin \left(\theta\right) y \right)} e^{-j \omega t}$$

Now the good part. Observing the above equation you may notice that on each dimension (Cartesian direction) you end up with a (co)sinusoidal component. Of course their "amplitude" is $$k_{i}$$. Do not be tempted to consider this to be the actual amplitude of each component, it's just a scaling factor and the final quantity ($$k_{i} \cos \left( \theta \right)$$, or similarly $$k_{i} \sin \left( \theta \right)$$), which is a constant for specific temporal frequency (keep in mind that $$k = \frac{\omega}{c} = \frac{2 \pi f}{c}$$ with $$c$$ being the speed of sound and $$f$$ the temporal frequency) and direction (this is dictated by the result of $$\cos \left( \theta \right)$$ or $$\sin \left( \theta \right)$$ which depends only on the angle of propagation $$\theta$$), is just a scalar that declares the spatial frequency of the plane wave on each direction of propagation. A simple example of this would be a plane wave traveling parallel to the $$x$$ axis. The expression is

$$p(x, y, t) = A e^{-j \left( k_{x} \cos \left(\theta\right) x + k_{y} \sin \left(\theta\right) y \right)} e^{-j \omega t} \implies p(x, y, t) = A e^{-j \left( k_{x} \cos \left(0 \right) x + k_{y} \sin \left(0 \right) y \right)} e^{-j \omega t} \implies p(x, y, t) = A e^{-j k_{x} x} e^{-j \omega t}$$

Please note that we are still in two dimensions. If you get the real part you end up with the well known (co)sinusoidal solution

$$p(x, y, t) = A \cos \left( \omega t - k_{x} x \right)$$

Below you can see an illustration of such a wave.

The picture shows a monochromatic plane wave with temporal frequency equal to $$100 Hz$$. Similarly, if you use a different $$\theta$$ and plug it in the general equation you end up with a plane wave traveling in the direction dictated by the $$\theta$$ parameter. An example of a monochromatic wave with the same temporal frequency and direction of travel $$60^{o}$$ is shown in the next image.

If you take a closer look, especially easy to notice on the boundaries of the screen, you may see that on each direction ($$x$$ and $$y$$ axes) there are (co)sinusoidal components forming. Those are depicted below for the plane wave of the last (second) picture (please note that the origin of the system lies in the top left corner. This is due to MATLAB's plotting function used. I apologise for the inconvenience)

Now, I believe it is quite clear that those two components constitute the spatial components of the plane wave. These are also utilised in the so called "Plane wave decomposition" methods if further simplification of the plane waves is needed. They become useful in 3D audio and array processing techniques. The frequency of the spatial components is the spatial frequency of each component governed by (actually being identical to) the respective wavenumbers ($$k_{x}$$ and $$k_{y}$$ in this case).

At this point, it is instructive to say (and easy to notice) that if you offset the constant directions by a value, the spatial component will experience a phase shift. For example, if the spatial components were not plotted for $$y = 0$$ and $$x = 0$$ but for some arbitrary $$x$$ and $$y$$ values, the corresponding spatial components would experience a phase shift. You can have a look at the respective spatial components by inspecting the 2D plot of the wave in the middle of the figure.

Although we have covered a lot of the road we haven't touched upon the temporal dimension yet. What you've seen show far is a snapshot of the monochromatic plane waves we have been talking about. Of course, if you were to take many such snapshots at different time instances you could create a "movie" of the propagating wave. The problem is that you cannot tell which way the wave is progressing towards. This is due to the fact that deal with the steady state and we cannot deduce any such information from what we have at hand. It's like floating in the middle of the ocean and moving up and down due to waves but with absolutely no mass transfer from the waves (this is essentially what linearised acoustic waves are). How would you be able to say which direction the waves are moving towards since you only know when you are at the peak and when you reach the lowest possible height? Well, the answer is you can't.

As for the velocity, it cannot be deduced from a temporal snapshot because time is "frozen". You look at the wave at a specific moment in time but in order to calculate speed you need at least two time instances to calculate the difference in position and divide by the duration ($$c = \frac{x}{t}$$). If you are allowed to let time change then you could possibly look for the duration between two maxima or minima (or whichever part of the cycle) and estimate the speed from that.

## Code

The code to create the plots I presented can be seen below

%% Set parameters
c = 343; % Speed of sound
ang = 60; % Angle of propagation
f = [100]; % Temporal frequency
Lx = 10; % length of x dimension of the domain
Ly = 10; % length of y dimension of the domain
res = 1e-2; % Resolution of the lengths

% Calculate Variables
Lx = 0:res:Lx - res; % Calculate the x-axis
Ly = Ly - res:-res:0; % Calculate the y-axis
w = 2 * pi * f; % Radial frequency
k = w/c; % Wavenumber
kx = k .* sin(deg2rad(ang)); % x component of wavenumber
ky = k .* cos(deg2rad(ang)); % y component of wavenumber

%% Calculate signal
sig = zeros(2, length(Lx), length(Ly)); % Pre-allocate for speed

% Go through the positions (this is quite slow but is a straight forward implementation)
for xx = 1:length(Lx) % x dimension
for yy = 1:length(Ly) % y dimension
sig(1, xx, yy) = sin(w(1) - kx(1) * Lx(xx) - ky(1) * Ly(yy)); % Calculate the value of the plane wave on each position
end
end

%% Visualise
figure("Name", "Plane Wave")
imagesc(squeeze(sig(1, :, :)));
xticklabels([])
yticklabels([])

figure("Name", "Components")
subplot(2, 1, 1)
temp = squeeze(sig(1, :, :));
plot(temp(1, :), 'LineWidth', 5);
grid on
xticklabels([])
yticklabels([])
title("f = " + num2str(f) + "Hz, angle = " + num2str(ang) + "^o (x component, y = 0)", 'FontSize', 24)

subplot(2, 1, 2)
plot(temp(:, 1), 'LineWidth', 5);
grid on
xticklabels([])
yticklabels([])
title("f = " + num2str(1) + "Hz, angle = " + num2str(ang) + "^o (y component, x = 0)", 'FontSize', 24)

• P.S.: I believe you question is better suited to the Physics SE (physics.stackexchange.com). Commented Feb 11, 2022 at 2:29

So, for 1D signal, the wave directed at 90 degree while for 2D, the direction of propagation is 45 degree. [ N.B. unfortunately, my comment option is blocked!]

• So, for 1D signal, the wave directed at 90 degree while for 2D, the direction of propagation is 45 degree. [ N.B. unfortunately, my comment option is blocked!] Commented Feb 10, 2022 at 6:22
• $k=2\pi f/c$, the inner circle has smaller $k$ so higher $c$. The direction is given by atan2(ky, kx) as I already said. Draw a line from the origin to any point in the $k_x-k_y$ plane, that's the direction. Commented Feb 10, 2022 at 7:21
• That is not a wave. It's kx-ky-cross spectral density plot, not x-y-amplitude. Commented Feb 10, 2022 at 9:16
• And you are reading a scientific article. The best way to understand the figure is to read the article carefully. The author should have explained the figures. Commented Feb 15, 2022 at 9:52
• @ZaellixA, can you suggest how to find the time continuity while appending the wave at different time instants in the link provided stackoverflow.com/questions/71654062/… Commented Mar 29, 2022 at 5:52

I can change the above python script to plot the 2 dimensional plane wave. The script is as follows:

import numpy as np
import random
import matplotlib.pyplot as plt
import cmath
from scipy import signal
import math

f = 19;         # Temporal frequency
fs = 100;       # sampling frequency
Ts = 1/fs;      # sample period
t = np.arange(0, 25, Ts);   # time index
c = 40;         # speed of wave
w = 2*np.pi *f;     # angular frequency
k = 2.0 * np.pi * f / c  # wavenumber

#*** Here Nyquist sampling frequency w.r.t. space #####
resolution = 0.01   # Spatial resolution
x = np.arange(-1.5, 1.5, resolution)
y = np.arange(-1.5, 1.5, resolution)
X, Y = np.meshgrid(x,y)

psi = np.pi / 6 #% direction of propagation
kx = k * np.cos(psi)
ky = k * np.sin(psi)
t= np.arange(0, 10, 1)

#2D sine wave
d2_sig = np.sin(-(kx*X+ky*Y)+w*t[5])

plt.figure()
c = plt.imshow(d2_sig, cmap='seismic', vmin=d2_sig.min(), vmax=d2_sig.max(), extent=[x.min(), x.max(), y.min(), y.max()],
interpolation='nearest', origin='lower')
plt.colorbar(c)
plt.xlabel("Distance, $$X$$ [Meter]", fontsize=18)
plt.ylabel("Distance,$$Y$$ [Meter]", fontsize=18)
plt.title(f'Progressive wave', weight="bold")


N.B. The above python script has been presented with the help of the above matlab scripts provided by ZaellixA & ZR Han.

• Your code is correct. It is very similar to what both answers (mine and ZR Han's) have presented. If you would like to add more information to the question, please edit it and don't post it as an answer. Otherwise, if this is the exact answer you have been looking for, mark it appropriately. Commented Apr 1, 2022 at 16:09
• @ZaellixA, can you please have a look at the link dsp.stackexchange.com/questions/82903/…. Sorry, the script is in Python. Commented May 5, 2022 at 11:06