python - MUSIC算法Spectrum Python实现

Question

我正在开发一个小型雷达项目，可以测量心脏和胸部产生的多普勒频移。由于我事先知道源的数量，因此我决定选择 MUSIC 算法进行频谱分析。我正在获取数据并将其发送到 Python 进行分析。然而，我的 Python 代码表示，具有两个频率为 1 Hz 和 2 Hz 的混合正弦波的信号的所有频率的功率是相等的。我的代码与示例输出链接在此处:

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


"""

以下是音乐算法的公式:

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