I am trying to use Frequency-wavenumber analysis for estimating speed of sound using two-dimensional data obtained via channel of arrays. The traditional technique is to use FFT to transform the data from distance-time domain to wavenumber-frequency domain and detect the slope of the bright line in the FK plot which in turn represents the speed of sound (for acoustic data). I have recently read MUSIC algorithm https://msol.people.uic.edu/ECE531/papers/Multiple%20Emitter%20Location%20and%20Signal%20Parameter%20Estimation.pdf as a technique that can be used for estimating power spectrum. Below is my code of MUSIC algorithm which is tested on a simple example. The code works well for detecting the frequencies of the signal. I would like to know how to implement that technique for estimating 2D power spectrum if I have multichannel data.
from scipy import signal import numpy as np from numpy import linalg as LA from matplotlib import pyplot as plt import scipy N = 1000 nfft = 2**16 f1 = 20 f2 = -8 fs = 100 t=np.arange(0,N, 1) c1 = np.exp(1j*2*np.pi*t*f1/fs) c2 = np.exp(1j*2*np.pi*t*f2/fs) snr = 10 stdev = 1/((10**(snr/10))**0.5) data = c1 + c2 + 1/(np.sqrt(2)) * (stdev * np.random.normal(0, 1, N) + 1j * stdev * np.random.normal(0, 1, N)) p = 2 m = len(data) acf = np.convolve(data,np.conj(data)[::-1]) center = int(np.ceil(len(acf)/2) - 1) Rxx=acf[center:] Rx = scipy.linalg.toeplitz(Rxx,np.hstack((Rxx, np.conj(Rxx[1:])))) D, V = LA.eig(Rx) i = np.argsort(D) Px = 0 for index in range(0, m-p): Px = Px + np.abs(np.fft.fftshift(np.fft.fft(V[:, i[index]], nfft))) Px = np.reciprocal(Px) Px = 10*np.log10(Px) w = np.arange(-fs/2, fs/2, fs/nfft) plt.plot(w,Px) plt.title('MUSIC Algorithm for Peak Detection') plt.xlabel('Frequency [Hz.]') plt.show()