Convert MUSIC algorithm to EV

I write some code, but result looks very strange:

import numpy as np
import scipy.linalg as LA

def music(X):
# create a signal
N = 64
m = np.arange(N)
P = 4
x = np.zeros((N, P+1), dtype='complex')

x[:,0] = 1.0*np.exp(1j*3.0*2.0*np.pi*m/N)
x[:,1] = 0.8*np.exp(1j*5.0*2.0*np.pi*m/N)
x[:,2] = 0.3*np.exp(1j*8.0*2.0*np.pi*m/N)
x[:,3] = 0.5*np.exp(1j*12.0*2.0*np.pi*m/N)
x[:,4] = 1.5 * np.random.normal(0, 1, size=N)

X = np.sum(x, axis=1)

# MUSIC

XT = np.mat(LA.toeplitz(X))
Rxx = XT*XT.H

en, ev = LA.eig(Rxx)
_, S, _ = LA.svd(Rxx)

idx = en.argsort()[::-1]
en = en[idx]

EV = np.mat(ev[:,P:]).T
EVSignal = np.mat(ev[:,:P])

ref = np.zeros(N, dtype=complex)

f = np.arange(0, 32, 0.5)

Mmusic = np.zeros(f.size)
Ev_psd_alg = np.zeros(f.size)

for i in range(f.size):
ref = np.exp(1j*2.0*np.pi*m*f[i]/N)
REF = np.mat(ref)
Mult=REF*EV.H*EV*REF.H
Mmusic[i]=np.absolute(1/Mult)

Ev_psd_alg[i] = np.absolute(1/( Mult / 1/(S[i]) )) # EV I'm do not sure about it

Ev_psd_alg[i] = np.abs(1/Mult / S[i])

Mmusic = 10*np.log10(Mmusic/np.max(Mmusic))
Ev_psd_alg = 10*np.log10(Ev_psd_alg/np.max(Ev_psd_alg))

return Mmusic, Ev_psd_alg


• Can you explicitly explain the part you found strange and your corresponding comments so that we can help you better :) Jan 11 at 7:49
• @OkanErturk EV spectrum looks like MUSIC rotated spectrum. I have tried music and ev methods from spectrum python library and with them no was such problem. Jan 11 at 10:31
• Can you link to the EV method?
– Royi
Jan 11 at 12:21
• @Royi pyspectrum.readthedocs.io/en/latest/_modules/spectrum/… if method == 'music': PSD = PSD + abs(Z)**2. elif method == 'ev' : PSD = PSD + abs(Z)**2. / S[I] I'm doubt that's code correct Jan 12 at 15:38
• So are you after reproducing the results of that package in your code?
– Royi
Jan 12 at 16:44