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)
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.title('MUSIC Algorithm for Peak Detection')
plt.xlabel('Frequency [Hz.]')

enter image description here



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.