# MUSIC Algorithm for 2D Power Spectrum

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[0], 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()