MUSIC算法谱Python实现

Question

我正在做一个小型雷达项目，它可以测量心脏和胸部产生的多普勒频移。由于我事先知道了源的数目，所以我决定选择MUSIC算法进行频谱分析。我正在获取数据并将其发送到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

从数学上讲，问题是i和f都是整数。因此，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