I am working on a small radar project that can measure the Doppler shift created by the heart and chest. Since I know the number of sources in advance, I decided to choose the MUSIC Algorithm for spectral analysis. I am acquiring data and sending it to Python for analysis. However, my Python code is saying that the power for ALL frequencies of a signal with two mixed sinusoids of frequency 1 Hz and 2 Hz is equal. My code is linked here with a sample output:
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
"""
Here is the formula for the MUSIC Algorithm:
https://drive.google.com/file/d/0B5EG2FEWlIZwYmkteUludHNXS0k/view?usp=sharing
Mathematically, the problem is that i and f are both integers. Thus, 2*π*i*f is an integral multiple of 2π. Allowing for a tiny bit of round-off error, this gives you a cosine very close to 1.0 and a sin very close to 0.0. These values yield virtually no variation in frequencyVector from one iteration to the next.
I also see a problem in that you set up your signal_eigen matrix, but never use it. Isn't the signal itself required by this algorithm? As a result, all you're doing is sampling the noise at intervals of 2πi.
Let's try chopping up one cycle into sRate evenly-spaced sampling points. This results in spikes at 0.24 and 0.76 (out of the range 0.0 - 0.99). Does this match your intuition about how this should work?
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
Output
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
I'm also unsure of the correct implementation; having more of the paper for formula context would help. I'm not certain about the range and sampling of the f values. When I worked on FFT software, f was swept over the wave form in small increments, typically 2π/sRate.
I'm not getting those distinctive spikes now -- not sure what I did before. I made a small parametrized change, adding a num_slice variable:
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
You can compute it however you like, of course, but the ensuing loop runs through just the one cycle. Here's my output:
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