1
$\begingroup$

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

enter image description here

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Can you explicitly explain the part you found strange and your corresponding comments so that we can help you better :) $\endgroup$
    – Okan Erturk
    Commented Jan 11, 2023 at 7:49
  • $\begingroup$ @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. $\endgroup$
    – Dmitry Sokolov
    Commented Jan 11, 2023 at 10:31
  • $\begingroup$ Can you link to the EV method? $\endgroup$
    – Royi
    Commented Jan 11, 2023 at 12:21
  • $\begingroup$ @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 $\endgroup$
    – Dmitry Sokolov
    Commented Jan 12, 2023 at 15:38
  • $\begingroup$ So are you after reproducing the results of that package in your code? $\endgroup$
    – Royi
    Commented Jan 12, 2023 at 16:44

0

Reset to default

Your Answer

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