MUSIC 算法谱 Python 实现答案

【问题标题】：MUSIC Algorithm Spectrum Python ImplementationMUSIC 算法谱 Python 实现
【发布时间】：2016-06-15 17:10:51
【问题描述】：

我正在开展一个小型雷达项目，该项目可以测量由心脏和胸部产生的多普勒频移。由于我事先知道了声源的数量，所以我决定选择 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


"""

MUSIC 算法的公式如下：

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

【问题讨论】：

标签： python signal-processing fft linear-algebra frequency-analysis

【解决方案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

【讨论】：

感谢您的回复。你是 100% 正确的，频率向量只是导致 cos 和 sin 分别为 1 和 0 的值。但是，我不确定如何根据图片中的公式来解决这个问题。此外，未使用 signal_eigen 矩阵，因为 MUSIC 算法基于信号与噪声正交的事实。通过在噪声向量和潜在频率之间进行点积，我们可以看到结果是正交还是零，这意味着它是信号子空间的一部分。为此，将信号子空间与噪声子空间分开。
当然，这里有一个使用 MUSIC 算法的书的链接：dsp-book.narod.ru/302.pdf，它从第 288 页开始
谢谢。我跑了一遍；看来 f 应该是 F-sub-i / F-sub-k，这似乎转化为 2πf / sRate ...？这意味着从 0 到 2π 进行一次运行，取自 sRate 采样点（省略任一端点）。
是的，这些结果看起来要好得多，但是基于初始正弦曲线，不应该只在 1 和 0.5 处出现尖峰吗？还有你为什么从 0 变成 sRate * p(2) + 1？
有点意思，我发现如果把sum1 += (abs(np.dot(np.asarray(frequencyVector).transpose(), np.asarray(noise_eigen[u] ) )))**2 通过特征值得到更清晰的图。看到 FFT 的表现比这好得多令人不安，所以我想我们一定忽略了一些错误。