I am writing a Direction of Arrival (DOA) estimator using the MUSIC algorithm in python and am having trouble getting the correct spectrogram output.
Specifically, my resulting graph comes out aliased as such:
On the left is the SYNTHESIZED signal, using the same source azimuth and elevation as the real dataset. On the right is the real-world signal.
These graphs looks like this when I am looking for 2 signals. If I am looking for only 1, they look like this:
My guess is that the reason this is happening is that I am trying to use a Hermitian, which should be reserved exclusively for complex numbers, when my signals are represented using only real values, without an imaginary one. The the real world dataset is presented as an NxM .wav where N is the number of channels (8) and M is the length of the signals (174800 or so) sampled at 44.1KHz. The data arrives in int16 format, which I then convert into floating point.
The eigenvectors I obtain from the covariance matrix look realistic, with 2 peaks, 1 twice as high as the other. The synthesized signal has only 1 peak. As seen in the graphs, if I limit to using all but the highest peak eigenvector, the spectrogram doesn't really show anything, while taking M-2 does provide peaks, but theyre aliased. I am certain this is aliasing and not a reverberation of some sort, because setting an extremely low wavelength variable for the steering matrix makes the real-world signal spectrogram perfectly symmetrical. If I shuffle microphone positions, I also never get more or less than 2 peaks.
The following are some code snippets to get an idea of how I implemented the algorithm:
The synthesized signal is created using these two functions as the core:
def emit(shape,C):
x = np.random.randn(shape)
x += np.random.randn(shape)
x *= np.sqrt(0.5)
x *= sqrt(C)
return x
generate(sources):
A = array.steering_matrix(sources, wavelength)
S = source_signal.emit((channel_count,n_snapshots),1.0)
Y = A @ S
N = noise_signal.emit((channel_count,n_snapshots),1.0)
Y += N
return Y
Once Y is returned, I do not tamper with it in any other way.
The full music function is as follows:
def music(A,x,k):
R = (x @ x.conj().T) / x.shape[1]
Ev, E = np.linalg.eigh(R)
En = E[:,:-k]
v = En.T.conj() @ A
return np.reciprocal(np.sum(v * v.conj(), axis=0).real),En
Steering matrix:
def steering_matrix(sensors,sources,wavelength):
sources = np.deg2rad(sources)
s = 2 * np.pi / wavelength
cos_el = np.cos(sources[:, 1])
cc = cos_el * np.cos(sources[:, 0])
cs = cos_el * np.sin(sources[:, 0])
D = s * (np.outer(sensors[:, 0], cc) +
np.outer(sensors[:, 1], cs) +
np.outer(sensors[:, 2], np.sin(sources[:, 1])))
A = np.exp(1j * D)
return A
I am nearly certain that this is related to complex numbers. If I synthesize a signal using the steering matrix, I get a signal that has mixed values on both the real and the imag portions, but my real world data has no imag portion. So if the implementation I am using for MUSIC requires some of the information to be retained on the imag portion, everything would mess up. I am specifically questioning my use of conj().
