我尝试从上面的等式(附上)实现空间谱
其中kX、kY为k空间中的网格点， C(w,r) - 第 i 个和第 j 个传感器之间的交叉光谱密度(这里它是一个大小为 ns * ns > 传感器数量的矩阵)。 x, y 是传感器之间的距离。 (nk - kx, ky 的网格密度)
我正在寻找上述等式的合适 python 实现。我有 34 个传感器，它们生成大小为 [row*column]=[n*34] 的数据。首先，我找到了每个传感器数据中的交叉光谱密度 (CSD)。然后对CSD值进行二维DFT得到空间谱。
*) 我不确定程序是否正确。 **) python 实现过程是否正确？ ***) 另外，如果有人提供一些相关的教程/链接，对我也会有帮助。
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)
abs_csd = np.abs(Pxy)
total_csd.append(abs_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
raw_data=np.reshape(np.random.rand(10000*34),(10000,34))
# 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, :]
#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
n = 34
f = 10 # frequency in Hz
c = 50 # wave speed 50, 80, 100, 200 m/s
k = 2.0*np.pi*f/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]
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")
#c = Wave Speed; 50, 80,100,200
cc = 2*np.pi*f /c *np.cos(np.linspace(0, 2*np.pi, 34))
cs = 2*np.pi*f /c *np.sin(np.linspace(0, 2*np.pi, 34))
plt.plot(cc,cs)
但是，通过使用改进的代码，我得到了与图 01 不同的图 02，其分辨率更高。
我添加了另外两个数字来与图 01 进行比较。当考虑范围 [-k, k] 时，绘图看起来像图 03
这是类似的 [w.r.t. XY-axis] to Fig. 01，我认为这个图是可以的，除了一些 K 空间遗漏。我希望这里存在一个需要解决的问题。
在图 04 中，我们考虑 k 空间范围 [-20k, 20k]，它看起来不错但没有与图 01 相似的轴。
我把更新图放如下:
任何人都可以帮我生成图形 01 或类似的类型吗？我对图 02 感到困惑。任何人都可以帮助我理解吗？提前致谢。
最佳答案
在我看来，您正在放大中央叶。这也可以解释为什么比例不是从 0 到 1。
如果我改变这些行:
kx = np.linspace(-20*k,20*k,nx) # space vector
ky= np.linspace(-20*k,20*k,ny) # space vector
然后我得到
它看起来更接近您要查找的内容。
为了提高分辨率，我做了一些重写以获得这张新图片。请参阅下面的更新代码。
注意:我仍然不确定这样做是否正确。
我使用的代码
# Code from https://stackoverflow.com/questions/70768384/right-method-for-finding-2-d-spatial-spectrum-from-cross-spectral-densities
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import cmath
# 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)
np.random.seed(12345)
raw_data=np.reshape(np.random.rand(10000*n),(10000,n))
f = 10 # frequency in Hz
c = 50 # wave speed 50, 80, 100, 200 m/s
k = 2.0*np.pi*f/c # wavenumber
kx = np.linspace(-20*k,20*k,n*10) # space vector
ky= np.linspace(-20*k,20*k,n*10) # space vector
# 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(Pxy) # output as list
return np.array(total_csd)
## Spatial Spectra:- DFT of the csd along two dimension
def DFT2D(data):
#data = np.asarray(data)
dft2d = np.zeros((len(kx),len(ky)), 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, :]
#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, n))
dft = DFT2D(fcsd) # Data or cross-correlation matrix
spec = np.abs(dft) #dft.real # 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 @10Hz', weight="bold")
关于python - 从 CSD 中寻找二维空间谱的正确方法，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/70768384/