search

Masoud-mk / 2D-DFT

2D dft for image processing in python
4 stars 2 forks source link

Right method for finding 2-D Spatial Spectrum from cross spectral densities #1

Closed ssalam21 closed 2 years ago

ssalam21 commented 2 years ago

I try to implement the spatial spectrum from the above equation (attached)

Where kX, kY are the grid points in k space, C(w,r) - cross spectral densities between the i'th and j'th sensor(here it is a matrix of size ns * ns >no. of sensors). x, y are distances between the sensors. (nk - grid density for kx, ky)

I look for suitable python implementation of the above equation. I've 34 sensors which generates data of size [rowcolumn]=[n34]. At first, I've found the cross spectral densities (CSD) of among the data of each sensor. Then 2D DFT is performed of the CSD values to get the spatial spectrum.

) I'm not sure whether the procedure is correct or not. ) Does the python implementation procedure correct? ) Also, if someone provides some relevant tutorial/link, it will also be helpful for me.

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

Finding cross spectral density (CSD)

fs=500 def csdMat(data): rows, cols = data.shape total_csd = []

for i in range(cols):

    for j in range(cols):
        f, Pxy = signal.csd(data[:,i], data[:,j], fs, nperseg=512)
        real_csd = np.real(Pxy)
        total_csd.append(real_csd)                     # output as list
        csd_mat = np.array(total_csd)
return csd_mat

Spatial Spectra:- DFT of the csd along two dimension

def DFT2D(data):

data = np.asarray(data)

dft2d = np.zeros((M,N), dtype=complex)
for k in range(len(kx)):
    for l in range(len(ky)):
        sum_matrix = 0.0
        for m in range(M):
            for n in range(N):
                e = cmath.exp(- 1j * ((kx[k] * dx[m]) / len(dx) + (ky[l] * dy[n]) / len(dy)))
                sum_matrix +=  data[m,n] * e
        dft2d[k,l] = sum_matrix
return dft2d

Call the seismic array

** Open .NPY files as an array

with open('res_array_1000f_131310.npy', 'rb') as f: arr= np.load(f) raw_data = arr[0:10000, :]

raw_data=np.reshape(np.random.rand(1000*34),(1000,34))

CSD of the seismic data

csd = csdMat(raw_data) print('Shape of CSD data', csd.shape)

CSD data of a specific frequency

csd_dat=csd[:, 11]
fcsd = np.reshape(csd_dat, (-1, 34)) fcsd.shape

data = coh_avg[:, :, 9]

n = 34 f = 10 # frequency in Hz c = 50 # wave speed 50, 80, 100, 200 m/s k = 2.0np.pif/c # wavenumber nx = n # grid density ny = n kx = np.linspace(-k,k,nx) # space vector ky= np.linspace(-k,k,ny) # space vector

Distance[Meter] between sensors

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]

dx = np.array(x); M = len(dx) dy = np.array(y) ; N = len(dy) X,Y = np.meshgrid(kx, ky)

dft = DFT2D(fcsd) # Data or cross-correlation matrix spec = dft.real # Spectrum or 2D_DFT of data[real part]

vmin = spec.min()

spec = spec/spec.max()

plt.figure() c = plt.imshow(spec, cmap ='seismic', vmin = spec.min(), vmax = spec.max(), extent =[kx.min(), kx.max(), ky.min(), ky.max()], interpolation ='nearest', origin ='lower') plt.colorbar(c) plt.rcParams.update({'font.size': 18}) plt.xlabel("Wavenumber, $K_x$ [rad/m]", fontsize=18) plt.ylabel("Wavenumber,$K_y$ [rad/m]", fontsize=18) plt.title(f'Spatial Spectrum @10Hz', weight="bold")`

ssalam21 commented 2 years ago

`import numpy as np from scipy import signal import matplotlib.pyplot as plt import cmath import cv2

Set up data

Distance[Meter] between sensors

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]

if (len(x) != len(y)): raise Exception('X and Y lengthd differ')

n = len(x) dx = np.array(x); M = len(dx) dy = np.array(y); N = len(dy)

Finding cross spectral density (CSD)

fs = 500 nperseg = 512 1 nfft = nperseg 1 # No. of FFT point def csdMat(data): rows, cols = data.shape total_csd = []

for i in range(cols):
    for j in range(cols):
        _, Pxy = signal.csd(data[:, i], data[:, j], fs=fs, window='hann', nperseg=nperseg, noverlap=None, nfft=nfft)
        # real_csd = np.real(Pxy)
        total_csd.append(Pxy)  # output as list
return np.array(total_csd)

Call the seismic array # Open .NPY files as an array

with open('res_array_1000f_131316.npy', 'rb') as f: arr = np.load(f) raw_data = arr[0:18000 * 1, :] # 1-hr data

np.random.seed(12345)

raw_data=np.reshape(np.random.rand(1800*n),(1800,n))

CSD of the seismic data

csd = csdMat(raw_data) print('Shape of CSD data', csd.shape)

k = 3.0 # np.pi f / c # wavenumber kx = np.linspace(-k, k, n 5) # space vector ky = np.linspace(-k, k, n 5) # space vector t = np.linspace(1, 34, 34) P=len(kx) Q=len(ky)

angle1=np.pi/8 angle2=np.pi/2 angle3=np.pi/6

Spatial Spectra:- DFT of the csd along two dimension

S = shift(34//6, 34)

def DFT2D(data): dft2d = np.zeros((P,Q), dtype=complex) A0 = 5 f=10 #Hertz (Hz) w = 2np.pi f w1 = 2np.pi 100 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])+ float(ky[l]dy[n]))))# cmath.exp(-1jwt[n]))

                A = A0 * np.sin(((kx[k]*dx[m])+(ky[l]*dy[n])) - w*t[n])  #+A0 * np.sin( ((ky[l]*dy[n])) - w*t[n])

                sum_matrix += data[m, n] * e
        dft2d[k,l] = sum_matrix
return dft2d

CSD data of a specific frequency

for i in range(5, 10): jj = i +1 csd_dat = csd[:, i] fcsd = np.reshape(csd_dat, (-1, n))

dft = DFT2D(fcsd)  # Data or cross-correlation matrix
spec = np.abs(dft)  # Spectrum or 2D_DFT of data[real part]
spec = spec / spec.max()
plt.figure()
c = plt.imshow(spec, cmap='seismic', vmin=spec.min(), vmax=spec.max(),
               extent=[kx.min(), kx.max(), ky.min(), ky.max()],
               interpolation='nearest', origin='lower')
plt.colorbar(c)
plt.rcParams.update({'font.size': 18})
plt.xlabel("Wavenumber, $K_x$ [rad/m]", fontsize=18)
plt.ylabel("Wavenumber,$K_y$ [rad/m]", fontsize=18)
plt.title(f'Spatial_Spectrum 3hr data f_{jj}hz', weight="bold")  # 'Channel %d' %i
figure = plt.gcf()
figure.set_size_inches(13, 8)`
ssalam21 commented 2 years ago

Now the script works fine as desired.