MUSIC Algorithm Spectrum Python Implementation

Question

Я работаю над небольшим проектом радара, который может измерять допплеровский сдвиг, создаваемый сердцем и сундуком. Поскольку я знаю количество источников заранее, я решил выбрать МУЗЫКАЛЬНЫЙ алгоритм спектрального анализа. Я получаю данные и отправляю их на Python для анализа. Однако мой код Python говорит, что мощность для ВСЕХ частот сигнала с двумя смешанными синусоидами частотой 1 Гц и 2 Гц равна. Мой код здесь связан с выходом образца:

from scipy import signal 
import numpy as np 
from numpy import linalg as LA 
import matplotlib.pyplot as plt 
import cmath 
import scipy 


N = 5 

z = np.linspace(0,2*np.pi, num=N) 
x = np.sin(2*np.pi * z) + np.sin(1 * np.pi * z) + np.random.random(N) * 0.3 # sample signal 

conj = np.conj(x); 

l = len(conj) 

sRate = 25 # sampling rate 

p = 2 
flipped = [0 for h in range(0, l)] 

flipped = conj[::-1] 


acf = signal.convolve(x,flipped,'full') 


a1 = scipy.linalg.toeplitz(c=np.asarray(acf),r=np.asarray(acf))#autocorrelation matrix that will be decomposed into eigenvectors 

eigenValues,eigenVectors = LA.eig(a1) 

idx = eigenValues.argsort()[::-1] 
eigenValues = eigenValues[idx] 
eigenVectors = eigenVectors[:,idx] 

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

eigenValues = eigenValues[idx]# soriting the eigenvectors and eigenvalues from greatest to least eigenvalue 
eigenVectors = eigenVectors[:,idx] 

signal_eigen = eigenVectors[0:p]#these vectors make up the signal subspace, by using the number of principal compoenets, 2 to split the eigenvectors 
noise_eigen = eigenVectors[p:len(eigenVectors)]# noise subspace 

for f in range(0, sRate): 
    sum1 = 0 

    frequencyVector = np.zeros(len(noise_eigen[0]), dtype=np.complex_) 

    for i in range(0,len(noise_eigen[0])): 
     frequencyVector[i] = np.conjugate(complex(np.cos(2 * np.pi * i * f), np.sin(2 * np.pi * i * f)))#creating a frequency vector with e to the 2pi *k *f and taking the conjugate of the each component 

    for u in range(0,len(noise_eigen)): 
     sum1 += (abs(np.dot(np.asarray(frequencyVector).transpose(), np.asarray( noise_eigen[u]))))**2 # summing the dot product of each noise eigenvector and frequency vector taking the absolute value and squaring 

    print(1/sum1) 
    print("\n") 


""" 
(OUTPUT OF THE ABOVE CODE) 
0.120681885992 
0 
0.120681885992 
1 
0.120681885992 
2 
0.120681885992 
3 
0.120681885992 
4 
0.120681885992 
5 
0.120681885992 
6 
0.120681885992 
7 
0.120681885992 
8 
0.120681885992 
9 
0.120681885992 
10 
0.120681885992 
11 
0.120681885992 
12 
0.120681885992 
13 
0.120681885992 
14 
0.120681885992 
15 
0.120681885992 
16 
0.120681885992 
17 
0.120681885992 
18 
0.120681885992 
19 
0.120681885992 
20 
0.120681885992 
21 
0.120681885992 
22 
0.120681885992 
23 
0.120681885992 
24 

Process finished with exit code 0 


""" 

Вот формула для MUSIC алгоритма:

https://drive.google.com/file/d/0B5EG2FEWlIZwYmkteUludHNXS0k/view?usp=sharing

Answer 1

Математически, проблема заключается в том, что я и е оба целые числа. Таким образом, 2 * π * i * f является целым кратным 2π. Учитывая крошечную ошибку округления, это дает вам косинус, очень близкий к 1.0, и грех, очень близкий к 0.0. Эти значения практически не изменяются в диапазоне frequencyVector от одной итерации к следующей.

У меня также возникла проблема в том, что вы создали свою signal_eigen матрицу, но никогда не используете ее. Не требуется ли этот сигнал для этого алгоритма? В результате все, что вы делаете, это выборка шума с интервалом в 2πi.

---

Давайте попробуем нарубить один цикл в sRate равномерно разнесенных точек отбора проб. Это приводит к увеличению пиков на 0,24 и 0,76 (вне диапазона 0,0-0,99). Это соответствует вашей интуиции о том, как это должно работать?

signal_eigen = eigenVectors[0:p] 
noise_eigen = eigenVectors[p:len(eigenVectors)]  # noise subspace 
print "Signal\n", signal_eigen 
print "Noise\n", noise_eigen 

for f_int in range(0, sRate * p + 1): 
    sum1 = 0 
    frequencyVector = np.zeros(len(noise_eigen[0]), dtype=np.complex_) 
    f = float(f_int)/sRate 

    for i in range(0,len(noise_eigen[0])): 
     # create a frequency vector with e to the 2pi *k *f and taking the conjugate of the each component 
     frequencyVector[i] = np.conjugate(complex(np.cos(2 * np.pi * i * f), np.sin(2 * np.pi * i * f))) 
     # print f, i, np.pi, np.cos(2 * np.pi * i * f) 

    # print frequencyVector 

    for u in range(0,len(noise_eigen)): 
     # sum the squared dot product of each noise eigenvector and frequency vector. 
     sum1 += (abs(np.dot(np.asarray(frequencyVector).transpose(), np.asarray(noise_eigen[u]))))**2 

    print f, 1/sum1 

Выход

Signal 
[[ -3.25974386e-01 3.26744322e-01 -5.24205744e-16 -1.84108176e-01 
    -7.07106781e-01 -6.86652798e-17 2.71561652e-01 3.78607948e-16 
    4.23482344e-01] 
[ 3.40976541e-01 5.42419088e-02 -5.00000000e-01 -3.62655793e-01 
    -1.06880232e-16 3.53553391e-01 -3.89304223e-01 -3.53553391e-01 
    3.12595284e-01]] 
Noise 
[[ -3.06261935e-01 -5.16768248e-01 7.82012443e-16 -3.72989138e-01 
    -3.12515753e-16 -5.00000000e-01 5.19589478e-03 -5.00000000e-01 
    -2.51205535e-03] 
[ 3.21775774e-01 8.19916352e-02 5.00000000e-01 -3.70053622e-01 
    1.44550753e-16 3.53553391e-01 4.33613344e-01 -3.53553391e-01 
    -2.54514258e-01] 
[ -4.00349040e-01 4.82750272e-01 -8.71533036e-16 -3.42123880e-01 
    -2.68725150e-16 2.42479504e-16 -4.16290671e-01 -4.89739378e-16 
    -5.62428795e-01] 
[ 3.21775774e-01 8.19916352e-02 -5.00000000e-01 -3.70053622e-01 
    -2.80456498e-16 -3.53553391e-01 4.33613344e-01 3.53553391e-01 
    -2.54514258e-01] 
[ -3.06261935e-01 -5.16768248e-01 1.08027782e-15 -3.72989138e-01 
    -1.25036869e-16 5.00000000e-01 5.19589478e-03 5.00000000e-01 
    -2.51205535e-03] 
[ 3.40976541e-01 5.42419088e-02 5.00000000e-01 -3.62655793e-01 
    -2.64414807e-16 -3.53553391e-01 -3.89304223e-01 3.53553391e-01 
    3.12595284e-01] 
[ -3.25974386e-01 3.26744322e-01 -4.97151703e-16 -1.84108176e-01 
    7.07106781e-01 -1.62796158e-16 2.71561652e-01 2.06561854e-16 
    4.23482344e-01]] 
0.0 0.115397176866 
0.04 0.12355071192 
0.08 0.135377011677 
0.12 0.136669716901 
0.16 0.148772917566 
0.2 0.195742574649 
0.24 0.237792763699 
0.28 0.181921271171 
0.32 0.12959840172 
0.36 0.121070836044 
0.4 0.139075881122 
0.44 0.139216853056 
0.48 0.117815494324 
0.52 0.117815494324 
0.56 0.139216853056 
0.6 0.139075881122 
0.64 0.121070836044 
0.68 0.12959840172 
0.72 0.181921271171 
0.76 0.237792763699 
0.8 0.195742574649 
0.84 0.148772917566 
0.88 0.136669716901 
0.92 0.135377011677 
0.96 0.12355071192 

---

Я также не уверен в правильности реализации; использование большего количества бумаги для контекста формулы помогло бы. Я не уверен в диапазоне и выборке значений f. Когда я работал над программным обеспечением FFT, f был пронесен по волновой форме небольшими приращениями, обычно 2π/sRate.

Я не получаю эти отличительные всплески сейчас - не знаю, что я делал раньше. Я сделал небольшие параметризованные изменения, добавив переменную num_slice:

num_slice = sRate * N 

for f_int in range(0, num_slice + 1): 
    sum1 = 0 
    frequencyVector = np.zeros(len(noise_eigen[0]), dtype=np.complex_) 
    f = float(f_int)/num_slice 

Вы можете вычислить его, как вам нравится, конечно, но последующий цикл выполняется только через один цикл. Вот мой вывод:

0.0 0.136398199883 
0.008 0.136583829848 
0.016 0.13711117893 
0.024 0.137893463111 
0.032 0.138792904453 
0.04 0.139633157335 
0.048 0.140219450839 
0.056 0.140365986349 
0.064 0.139926689416 
0.072 0.138822121693 
0.08 0.137054535152 
0.088 0.13470609994 
0.096 0.131921188389 
0.104 0.128879079596 
0.112 0.125765649854 
0.12 0.122750994163 
0.128 0.119976226317 
0.136 0.117549199221 
0.144 0.115546862203 
0.152 0.114021482029 
0.16 0.113008398728 
0.168 0.112533730494 
0.176 0.112621097254 
0.184 0.113296863522 
0.192 0.114593615279 
0.2 0.116551634665 
0.208 0.119218062482 
0.216 0.12264326497 
0.224 0.126873674308 
0.232 0.131940131305 
0.24 0.137840727381 
0.248 0.144517728837 
0.256 0.151830000359 
0.264 0.159526062508 
0.272 0.167228413981 
0.28 0.174444818009 
0.288 0.180621604818 
0.296 0.185241411664 
0.304 0.187943197745 
0.312 0.188619481273 
0.32 0.187445977812 
0.328 0.184829467764 
0.336 0.181300320748 
0.344 0.177396490666 
0.352 0.173576190425 
0.36 0.170171993077 
0.368 0.167379359825 
0.376 0.165265454514 
0.384 0.163786582966 
0.392 0.16280869726 
0.4 0.162130870823 
0.408 0.161514399035 
0.416 0.160719375729 
0.424 0.159546457646 
0.432 0.157875982968 
0.44 0.155693319037 
0.448 0.153091632029 
0.456 0.150251065569 
0.464 0.147402137481 
0.472 0.144785618099 
0.48 0.14261932062 
0.488 0.141076562538 
0.496 0.140275496354 
0.504 0.140275496354 
0.512 0.141076562538 
0.52 0.14261932062 
0.528 0.144785618099 
0.536 0.147402137481 
0.544 0.150251065569 
0.552 0.153091632029 
0.56 0.155693319037 
0.568 0.157875982968 
0.576 0.159546457646 
0.584 0.160719375729 
0.592 0.161514399035 
0.6 0.162130870823 
0.608 0.16280869726 
0.616 0.163786582966 
0.624 0.165265454514 
0.632 0.167379359825 
0.64 0.170171993077 
0.648 0.173576190425 
0.656 0.177396490666 
0.664 0.181300320748 
0.672 0.184829467764 
0.68 0.187445977812 
0.688 0.188619481273 
0.696 0.187943197745 
0.704 0.185241411664 
0.712 0.180621604818 
0.72 0.174444818009 
0.728 0.167228413981 
0.736 0.159526062508 
0.744 0.151830000359 
0.752 0.144517728837 
0.76 0.137840727381 
0.768 0.131940131305 
0.776 0.126873674308 
0.784 0.12264326497 
0.792 0.119218062482 
0.8 0.116551634665 
0.808 0.114593615279 
0.816 0.113296863522 
0.824 0.112621097254 
0.832 0.112533730494 
0.84 0.113008398728 
0.848 0.114021482029 
0.856 0.115546862203 
0.864 0.117549199221 
0.872 0.119976226317 
0.88 0.122750994163 
0.888 0.125765649854 
0.896 0.128879079596 
0.904 0.131921188389 
0.912 0.13470609994 
0.92 0.137054535152 
0.928 0.138822121693 
0.936 0.139926689416 
0.944 0.140365986349 
0.952 0.140219450839 
0.96 0.139633157335 
0.968 0.138792904453 
0.976 0.137893463111 
0.984 0.13711117893 
0.992 0.136583829848 
1.0 0.136398199883