# What is going wrong with the plot of 2D spatial spectrum at a specific frequency?

I've a set of 09 sensors in the following arrangement and the script for the sensor positions as follows:

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

arr_x = []
arr_y = []
sensor_count = 9 # No of sensors [mock network of sensors]
for i in range(1,sensor_count):
x = np.around(12*np.cos(2*np.pi*i/8), decimals=2)
y = np.around(12*np.sin(2*np.pi*i/8), decimals=2)
arr_x.append(x)
arr_y.append(y)
print('Numbers List X: ', arr_x)
print('Numbers List Y: ', arr_y)
arr_x.insert(0,0) # num_list.insert(index, num)
arr_y.insert(0,0)
print(f'Updated Numbers List X: {arr_x}')
print(f'Updated Numbers List Y: {arr_y}')

s_count = len(arr_x)
arr_x = np.array(arr_x) # Make list to an array
arr_y = np.array(arr_y)

if (len(arr_x) != len(arr_y)):
raise Exception('X and Y lengthd differ')

# Plot sensor position
labels = ["1st_S", "2nd_S", "3rd_S", "4th_S" , "5th_S", "6th_S", "7th_S", "8th_S", "9th_S"]
plt.figure()
for i, label in enumerate(labels):
plt.scatter(arr_x[i], arr_y[i],s=60, c='r', marker='o', cmap=None, norm=None)
text_object = plt.annotate(label, xy=(arr_x[i], arr_y[i]), xytext=(1,25), textcoords='offset points',ha='center')
for x,y in zip(arr_x, arr_y):
label = f"({x},{y})"
plt.annotate(label, xy=(x,y), xytext=(8,5), textcoords='offset points',ha='center')
plt.xlim(np.min(arr_x)-2,np.max(arr_x)+2)
plt.ylim(np.min(arr_y)-2,np.max(arr_y)+2)
plt.xlabel(" $$X$$-axis", fontsize=18)
plt.ylabel("$$Y$$-axis", fontsize=18)
plt.title(f'Sensor position in $$X-Y$$ coordinates',  loc='left', weight="bold")


The above code will give the sensor position as

Now, I want to inject a 2D sinusoidal wave into the sensor network say, $$Z = A \sin(k_x X + k_y Y - \omega t)$$, $$K = 2 \pi /\lambda$$, for example, assume $$f = 10$$ Hz, and $$\lambda = 12$$ Meter,

so the code for signal generator is as follows:

    # Signal Generation Parameters
#theta in degrees:direction of propagation, lamda in Meters, f=Frequency[Hz] & percentage is noise Percentage
[theta, lamda, f, percentage]=[90*1, 12*1, 10, 0.1]
amp=1   # Signal Amplitude
fs = 500;       # sample frequency
Ts = 1/fs;      # sample period
t = np.arange(0,10, Ts)  # time index
speed = f*lamda       # speed of wave
omega = 2*np.pi *f            # angular frequency
ki = omega/speed             # wave number
# Component of ki
kxx = ki*np.around(np.cos(math.radians(theta)), decimals=3) # use this method to find exact Zero
kyy = ki*np.around(np.sin(math.radians(theta)), decimals=3)

M = len(arr_x)
N = len(arr_y)
[xx, yy] = np.meshgrid(arr_x, arr_y)

arr = []
for count in range(len(t)):
p = amp*np.sin(kxx* xx + kyy*yy - omega * t[count])  # plane wave
arr.append(p)
sig = np.array(arr)
print('The shape of the array', sig.shape)

# Add Noise to the Signal
np.random.seed(123456)
noise = np.random.normal(0, sig.std(), sig.shape) * percentage
raw_data =  sig+noise # Only Signal present
print('The shape of the Noisy Signal', raw_data.shape)


Now, I want to find the spatial spectrum at 10 Hz where the equation of a 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.

So, the script for the cross spectral density (CSD) and the 2D DFT are as follows:

    # Function of CSD
fs = 500;       # sample frequency
nperseg = fs * 1  #fs = 200
nfft = fs * 1  # No. of FFT point
def csdMat(data):
dat, rows, cols = data.shape # For 3D data
#dat, cols = data.shape   # For 2D data
total_csd = []
for i in range(cols):
col_csd =[]

for j in range( cols):
freq, Pxy = signal.csd(data[:,i,  i], data[:,i, j], fs=fs, window='hann', nperseg=100, noverlap=70, nfft=6200) # nfft=6200
# real_csd = np.real(Pxy)
#print(Pxy.shape)
col_csd.append(Pxy)  # output as list
#col_csd = np.array(col_csd, dtype=np.ndarray)
total_csd.append(col_csd)
pxy = np.array(total_csd)
return freq, pxy

# Finding cross spectral density (CSD)
freq, csd = csdMat(raw_data)
print('The shape of the csd data', csd.shape)

kf=2*omega/speed #2*2*np.pi/lamda
kx = np.linspace(-kf, kf, N * 15)  # space vector
ky = np.linspace(-kf, kf, N * 15)  # space vector

def DFT2D(data):

P=len(kx)
Q=len(ky)
dft2d = np.zeros((P,Q), dtype=complex)
for k in range(P):
for l in range(Q):
sum_matrix = 0.0
for m in range(M):
for n in range(N): #
e = cmath.exp(-1j*(float(kx[k]*(dx[m]-dx[n])+ float(ky[l]*(dy[m]-dy[n])))))#* cmath.exp(-1j*w*t[n]))
sum_matrix += data[m, n] * e
#print('sum matrix would be', sum_matrix)
#print('sum matrix would be', sum_matrix)
dft2d[k,l] = sum_matrix
return dft2d

dx = arr_x[:]; dy = arr_y[:]
# jj = freq[124]    # Index 124 corresponds to 10Hz


now if I draw the 2D K-space plot with the following script

    # Plot of Spatial Spectrum
for fr_count in range(124, 148*1, 24):
jj = freq[fr_count]    # Index 124 corresponds to 10Hz
jj= np.around(jj) #fcsd = np.reshape(csd_dat, (-1, N))
spec_csd = csd[:,:, fr_count]
print('The shape of the freq csd data', spec_csd.shape)
dft = DFT2D(spec_csd) # For s0,s1 dft = DFT2D(csd[0:2,0:2,124]); M=2 & N=2
spec1 = np.real(dft)  # Spectrum or 2D_DFT of data[real part]
np.seterr(invalid='ignore')
spec = spec1 / spec1.max() # Normalization

plt.figure()
c = plt.imshow(spec, cmap='seismic', vmin=0, vmax=spec.max(),
extent=[kx.min(), kx.max(), ky.min(), ky.max()],
interpolation='nearest', origin='lower') #
plt.rcParams.update({'font.size': 18})
cbar = plt.colorbar(c)
cbar.set_label(f'Spatial Spectrum at {jj}Hz', rotation=90, fontsize=16,labelpad=15)#, weight="bold")
plt.xlabel("Wavenumber, $$K_x$$ [rad/m]", fontsize=18)
plt.ylabel("Wavenumber,$$K_y$$ [rad/m]", fontsize=18)
plt.title(f'Sig_L={lamda}M_a={theta}deg@ freq_{jj}Hz_$$S_1$$ to $$S_({s_count})$$', weight="bold")  # 'Channel %d' %i
figure = plt.gcf()
figure.set_size_inches(13, 8)

cc = 2*np.pi*(f/speed) *np.cos(np.linspace(0, 2*np.pi, 180))
cs = 2*np.pi*(f/speed) *np.sin(np.linspace(0, 2*np.pi, 180)) # speed = f*lamda
plt.plot(cc,cs)


I get the following figure . From the plot I've following questions: Q1. Do the plot really represent the spatial spectrum at 10 Hz? Q2. How do I find the K-space vector from the figure? Q3. If the figure is not correct, then how do I get a correct figure which represent the $$K_x$$ and $$K_y$$ values? Please guide if I do anything wrong.

The position of the red blob at the center means $$K_x=0.0$$ , $$K_y=0.0$$ while the injected wave has $$K_x=0.52* cos(90 deg.)=0.0$$ and $$K_y=0.52* sin(90 deg.)=0.52$$. Can anybody make me understand? Thanks.

The following link could be supportive if anyone need it, link given https://iopscience.iop.org/article/10.1088/1361-6382/ac348a/pdf

As @Hilmer suggested, I work for only 02 sensors [1st and 2nd sensor] the figure is and spatial spectrum for 02 sensors @0 degree is @90 degree as in fig

@180 degree as in fig

Again, from the figure why the value of $$K_x$$ and $$K_y$$ don't match as of the injected wave number $$K_x$$ and $$K_y$$? Can anybody please explain where did I do the mistake?

@Hilmar suggestion, so I go for a combination of pair of sensors [sensor 1 & sensor 2, sensor 1 & sensor 3, and sensor 2 & sensor 3] from a 3 sensor network. position of 3 sensors arr_x = $$[0, 8.49, 8.49]$$ ; arr_y = $$[0, 8.49, -8.49]$$

Case A) First I consider 1st & 2nd sensor, position at arr_x = $$[0, 8.49]$$ ; arr_y = $$[0, 8.49]$$, in that case the spatial spectrum looks like where Kx value is matched with the injected value [Kx=0.52 & Ky = 0.0]

Case B) Then I consider 1st & 3rd sensor position at arr_x = $$[0, 8.49]$$ ; arr_y = $$[0, -8.49]$$, in that case the spatial spectrum looks like where Kx value does not match with the injected value [Kx=0.52 & Ky = 0.0]

Case C) Then I consider 2nd & 3rd sensor position at arr_x = $$[8.49, 8.49]$$ ; arr_y = $$[8.49, -8.49]$$, in that case the spatial spectrum looks like where Kx value does not match with the injected value [Kx=0.52 & Ky = 0.0]

My question is why the wavenumber in Case (B) & Case (C) do not match with the injected value? Thanks in advance.

• That's a lot to digest in one go. Suggestion: start with two sensors and a bunch easy test signals where you already know the answer (0 deg incident, low frequency). Than work your way up one step at a time May 5 at 11:49
• @Hilmar, thanks for your comment. I've attached the plot for 02 sensors (0 degree, 90 degree and 180 degree), the plot looks similar to a plane wave. Still, how do I cross-check the value of $K_x$ and $K_y$? Could you please tell me how to match the value of wavenumber from the injected wave to the 2D plots? For me, the approach looks right but don't get expected results, look for a solution. Thanks again. May 6 at 5:05
• I am not sure what exactly you do with the meshgrid command. I thought that you should somehow get the coordinates of the sensors if you want to calculate the signal in those positions but you already have those. Furthermore, the output of the meshgrid command does not return the position of the sensors. You can see that yourself if you type plt.scatter(xx, yy) (of course accompanied by the appropriate commands to show the plot). Would you care to clarify that point a bit? I will try to look further into your code and see what I can find.... May 6 at 13:13
• @ZaellixA, "meshgrid function is used to create a rectangular grid out of two given one-dimensional arrays representing the Cartesian indexing or Matrix indexing". As I want to plot a 2D Sine wave, so I use meshgrid. Also, you can find a nice example of meshgrid here stackoverflow.com/questions/23768618/… May 6 at 13:31
• @Alan22 I do know what meshgrid does. My point is that in the way you use it, it provides neither the position of the sensors nor a uniform grid. So, I am wondering exactly why you use it the way you do. I am trying to follow your way of thinking to understand your code. In addition to that, why exactly do you calculate the $4 \pi/\lambda$ value (just before you define the 2D DFT)? May 6 at 13:54