0
$\begingroup$

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 Fig. 01

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 fig.02

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 Fig. 03. 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 02 sensor position and spatial spectrum for 02 sensors @0 degree is fig. spatial spectrum for 02 sensors @90 degree as in fig Fig.: 2 sensors & propagation angle=90 degree

@180 degree as in fig Fig: 2 sensors & angle of propagation = 180 degree

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. 03 sensor position 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 Fig. : spectrum of 1st & 2nd sensor 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 Fig.: spatial spectrum of 1st & 3rd sensor 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 Fig.: spatial spectrum of 1st & 3rd sensor 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.

Improve this question
$\endgroup$
10
  • 2
    $\begingroup$ 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 $\endgroup$
    – Hilmar
    May 5 at 11:49
  • $\begingroup$ @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. $\endgroup$
    – Alan22
    May 6 at 5:05
  • $\begingroup$ 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.... $\endgroup$
    – ZaellixA
    May 6 at 13:13
  • $\begingroup$ @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/… $\endgroup$
    – Alan22
    May 6 at 13:31
  • $\begingroup$ @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)? $\endgroup$
    – ZaellixA
    May 6 at 13:54

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.