2014-01-18 35 views
0

我想在matplotlib在这个库中创建像一个频谱,它们返回的PSD(功率谱密度),而不是绝对量等。计算PSD

我要实现的PSD到我自己的代码(用C++),我一直在寻找的源代码matplotlib它看起来是这样的:

result, windowVals = apply_window(result, window, axis=0, 
            return_window=True) 
    result = np.fft.fft(result, n=pad_to, axis=0)[:numFreqs, :] 

    result = np.conjugate(result) * result 

    # calculations for PSD 

    result /= (np.abs(windowVals)**2).sum() 

    result[1:-1] *= scaling_factor 

链接,可以发现:Here

现在我想办同样,我达到了平等t我已经应用了窗口,并执行了DFT。从那里,我认为我错了。

这是我计算PSD的代码。

的std ::矢量CalculatePSD(标准::矢量&丘壑) {

std::vector<double> hanning = getHanningWindow(vals.size()); 

std::vector<double> vals2(vals.size()); 

for(unsigned i=0; (i < vals.size()); i++) 
{ 
    double v = vals[i].re * (-1 * vals[i].im); 
    vals2[i] = v * v; 
} 

double vi = 0.0; 
for(unsigned i=0; (i < vals2.size()); i++) 
{ 
    vi += abs(hanning[i]); 
} 

for(unsigned i=0; (i < vals2.size()); i++) 
{ 
    vals2[i] /= vi; 

} 

for(unsigned i=0; (i < vals2.size()); i++) 
{ 
    vals2[i] *= 2; 
} 
return vals2; 
} 

我的结果,我得到:

[ 2.66926000e-10 4.71270000e-10 3.25024000e-10 1.86008000e-10 
    8.68276000e-11 2.53098000e-11 7.92737000e-13 5.29274000e-12 
    1.90057000e-11 3.27736000e-11 4.42093000e-11 4.60725000e-11 
    4.98549000e-11 6.08991000e-11 5.59033000e-11 2.26625000e-11 
    7.17159000e-13 1.65150000e-12 1.31530000e-12 7.86863000e-14 
    1.59211000e-13 2.80960000e-13 1.25471000e-13 1.35225000e-11 
    1.09812000e-10 3.49411000e-10 5.96210000e-10 6.43539000e-10 
    4.85667000e-10 2.46446000e-10 7.61480000e-11 7.70735000e-12 
    3.79924000e-13 4.16633000e-13 9.22683000e-12 6.23747000e-11 
    1.80818000e-11 5.13819000e-11 7.20294000e-10 2.46505000e-09 
    4.42844000e-09 4.77968000e-09 3.06048000e-09 1.12040000e-09 
    1.15313000e-10 4.78101000e-11 2.96048000e-10 3.93216000e-10 
    2.66339000e-10 1.01348000e-10 3.13012000e-11 3.59365000e-11 
    3.00758000e-11 2.44970000e-11 4.19013000e-12 1.92963000e-12 
    7.62176000e-13 9.83650000e-14 4.47749000e-15 3.29484000e-19 
    7.46089000e-14 3.99184000e-12 4.75458000e-11 3.36518000e-10 
    9.34579000e-10 9.14683000e-10 1.34570000e-10 6.47235000e-11 
    1.21893000e-10 1.44874000e-11 5.37149000e-15 1.93477000e-11 
    1.15593000e-10 1.27684000e-10 2.10146000e-11 9.34294000e-13 
    8.85961000e-12 6.35360000e-12 3.88437000e-13 2.75095000e-13 
    1.07261000e-11 9.35213000e-12 4.16280000e-11 1.52995000e-11 
    9.86730000e-12 2.04171000e-11 9.69115000e-12 1.81067000e-12 
    1.64522000e-13 4.01619000e-12 2.81457000e-13 5.86617000e-13 
    1.33937000e-11 3.16438000e-11 8.83417000e-11 1.24155000e-10 
    9.36973000e-11 1.89000000e-12 4.76718000e-14 3.46776000e-13 
    3.88868000e-12 1.01426000e-12 7.87006000e-13 1.02698000e-10 
    6.80669000e-12 3.22014000e-12 9.92309000e-12 1.17853000e-10 
    1.20818000e-11 2.20125000e-12 8.78943000e-12 3.34323000e-10 
    4.28139000e-10 3.57678000e-10 1.57808000e-10 2.05267000e-11 
    8.98399000e-11 1.12894000e-11 3.21320000e-14 7.77191000e-12 
    5.91681000e-13 9.92243000e-16 2.10894000e-11 2.19397000e-12 
    1.14148000e-12 3.41732000e-12 1.81439000e-12 1.46689000e-11] 

而处理结果在Python:

[ 6.44554713e-04 2.26979569e-02 1.48395306e-02 1.39560086e-02 
    1.70585613e-02 4.24042116e-04 4.10722082e-04 1.77314474e-02 
    5.48046037e-03 6.86724979e-03 1.33342952e-02 5.45918807e-04 
    1.42011959e-06 8.15283041e-03 3.02976247e-02 2.95310636e-02 
    2.69222586e-02 2.70161073e-04 4.27988811e-04 8.22069685e-03 
    1.14550280e-03 5.94684341e-03 5.03412155e-03 2.39065158e-04 
    1.88851349e-03 1.63618611e-02 1.02155767e-02 5.56409334e-03 
    2.03783039e-02 1.30646965e-03 7.83925381e-03 6.58153969e-04 
    8.58222471e-05 4.90329132e-03 9.27321780e-03 2.18878971e-02 
    7.80419597e-03 1.65506496e-05 2.12233732e-03 3.48564618e-02 
    3.04324943e-02 1.14097124e-02 1.83163044e-02 5.53528648e-04 
    1.72024876e-03 1.05496508e-02 1.22350425e-02 6.81764861e-03 
    2.18181750e-02 1.25305967e-04 3.45533908e-04 2.52806605e-02 
    2.79032703e-02 3.30741745e-02 8.92045889e-03 1.43861624e-04 
    1.37729407e-03 4.40048633e-02 4.43466583e-02 3.21348174e-02 
    1.97845126e-02 4.76052263e-05 1.90059116e-03 1.36124930e-02 
    3.08483724e-02 3.18817777e-02 4.22224299e-02 6.48991341e-04 
    3.37298579e-04 1.92796091e-02 3.26384995e-02 7.44582047e-03 
    2.75911372e-02 8.24143406e-05 6.19298800e-04 2.18909904e-02 
    7.85534253e-03 1.35622242e-02 1.64534364e-02 4.24610034e-07 
    1.25296965e-04 3.85049720e-03 1.56208315e-02 1.51067447e-02 
    1.15560295e-02 1.91845524e-02 1.51484986e-02 3.68090803e-05 
    4.56878093e-04 1.32583767e-02 2.67413477e-02 2.12116190e-02 
    1.40731136e-02 7.46782595e-06 7.56130481e-04 1.11894743e-02 
    3.82556474e-02 2.20488800e-02 1.18449472e-02 6.41610843e-05 
    8.93214478e-04 1.44705708e-02 8.95544599e-03 8.24627650e-03 
    1.54125088e-02 3.82922435e-07 1.21567170e-03 3.66207393e-02 
    2.52421164e-02 2.79258696e-02 2.42711875e-02 4.41070028e-04 
    9.18506931e-04 2.29391748e-02 2.93676503e-06 3.51546485e-02 
    1.53622376e-03 3.93588210e-05 3.35935113e-04 1.74232319e-02 
    1.89744096e-02 1.02178421e-02 1.54763125e-02 7.24992746e-05 
    3.46909205e-04 2.41130633e-02 1.96140922e-02 1.94479820e-02 
    1.02416629e-02 1.08582494e-04 5.91329398e-04 8.53889890e-03 
    8.05471223e-03 1.78734494e-02 2.12358089e-02 2.74402258e-02 
    1.86277585e-02 2.97114647e-06 1.58242458e-03 1.77131478e-02 
    1.03301989e-03 1.17236867e-02 2.70723000e-02 2.45157582e-04 
    1.70524416e-04 1.12361676e-03] 

我可能会出现错误的地方转换th e代码,因为我不太了解Python。如果有人发现我可能会出错的地方,它会对我有很大的帮助。

编辑:

complex<double> foo[vals.size()]; 
for(unsigned i=0; (i < vals.size()); i++) 
{ 
    foo[i] = complex<double>(vals[i].re, vals[i].im); 
} 

std::vector<double> vals2(vals.size()); 

for(unsigned i=0; (i < vals2.size()); i++) 
{ 
    double d =(std::conj(foo[i])*std::real(foo[i])); 
} 

回答

0

你的函数返回vals(参照修改它),而你似乎是存储结果vals2,这是获取在函数结束时删除局部变量。

编辑

是的,这似乎是错误可能结合 - 为什么不使用std ::复杂的功能复制Python代码,因为它是而不是重新发明轮子:

for(unsigned i=0; (i < vals.size()); i++) 
{ 
    vals2[i]= (vals[i].conj()*vals[i]).re(); 
} 
+0

道歉,我已经更新了 – Phorce

+0

我认为问题出现在'result = np.conjugate(result)* result'的问题中,我不认为我是以正确的方式对复数向量进行平方。 – Phorce

+0

@ user1326876请参阅编辑... –