python 信号処理　エラー

前提

信号に対して周波数解析を行っています。
kaggleのコードと同じことを実行してもエラーが出てしまいます。
kaggleのコードをコピペすると、実行できましたが、信号の中身を私の信号にするとエラーが出てしまいます。

実現したいこと

エラーを解決し、実行可能状態にしたいです。

kaggleのコード

python
1def frequency_domain(r, fs=4):
2    # Estimate the spectral density using Welch's method
3    fxx, pxx = signal.welch(x=rr_interpolated, fs=fs)
4    
5    '''
6    Segement found frequencies in the bands 
7     - Very Low Frequency (VLF): 0-0.04Hz 
8     - Low Frequency (LF): 0.04-0.15Hz 
9     - High Frequency (HF): 0.15-0.4Hz
10    '''
11    cond_vlf = (fxx >= 0) & (fxx < 0.04)
12    cond_lf = (fxx >= 0.04) & (fxx < 0.15)
13    cond_hf = (fxx >= 0.15) & (fxx < 0.4)
14    
15    # calculate power in each band by integrating the spectral density 
16    vlf = trapz(pxx[cond_vlf], fxx[cond_vlf])
17    lf = trapz(pxx[cond_lf], fxx[cond_lf])
18    hf = trapz(pxx[cond_hf], fxx[cond_hf])
19    
20    # sum these up to get total power
21    total_power = vlf + lf + hf
22     # find which frequency has the most power in each band
23    peak_vlf = fxx[cond_vlf][np.argmax(pxx[cond_vlf])]
24    peak_lf = fxx[cond_lf][np.argmax(pxx[cond_lf])]
25    peak_hf = fxx[cond_hf][np.argmax(pxx[cond_hf])]
26
27    # fraction of lf and hf
28    lf_nu = 100 * lf / (lf + hf)
29    hf_nu = 100 * hf / (lf + hf)
30    
31    results = {}
32    results['Power VLF (ms2)'] = vlf
33    results['Power LF (ms2)'] = lf
34    results['Power HF (ms2)'] = hf   
35    results['Power Total (ms2)'] = total_power
36
37    results['LF/HF'] = (lf/hf)
38    results['Peak VLF (Hz)'] = peak_vlf
39    results['Peak LF (Hz)'] = peak_lf
40    results['Peak HF (Hz)'] = peak_hf
41
42    results['Fraction LF (nu)'] = lf_nu
43    results['Fraction HF (nu)'] = hf_nu
44    return results, fxx, pxx
45
46
47# sample rate for interpolation
48fs = 4.0
49steps = 1 / fs
50
51# now we can sample from interpolation function
52xx_2 = np.arange(1, np.max(x_2), steps)
53
54rr_interpolated = f_2(xx_2)
55print("Frequency domain metrics:")
56results_2, fxx_2, pxx_2 = frequency_domain(rr_interpolated)
57
58for k, v in results_2.items():
59    print("- %s: %.2f" % (k, v))
60#結果
61Frequency domain metrics:
62- Power VLF (ms2): 3.64
63- Power LF (ms2): 7.07
64- Power HF (ms2): 7.81
65- Power Total (ms2): 18.53
66- LF/HF: 0.91
67- Peak VLF (Hz): 0.02
68- Peak LF (Hz): 0.09
69- Peak HF (Hz): 0.28
70- Fraction LF (nu): 47.51
71- Fraction HF (nu): 52.49
72

kaggleの信号

python
1print(rr_interpolated)
2#[59.12206572 58.903625   58.67572032 ... 75.11356874 74.51171895
3 70.1799672 ]
4print(xx_2)
5#[  1.     1.25   1.5  ... 298.75 299.   299.25]

私のコード

python
1def frequency_domain(rri, fs=40):
2    # Estimate the spectral density using Welch's method
3    fxx, pxx = signal.welch(x=f_sci(xx), fs=fs)
4    
5    '''
6    Segement found frequencies in the bands 
7     - Very Low Frequency (VLF): 0-0.04Hz 
8     - Low Frequency (LF): 0.04-0.15Hz 
9     - High Frequency (HF): 0.15-0.4Hz
10    '''
11    #cond_vlf = (fxx >= 0) & (fxx < 0.04)
12    cond_lf = (fxx >= 0.04) & (fxx < 0.15)
13    cond_hf = (fxx >= 0.15) & (fxx < 0.4)
14    
15    # calculate power in each band by integrating the spectral density 
16    #vlf = trapz(pxx[cond_vlf], fxx[cond_vlf])
17    lf = trapz(pxx[cond_lf], fxx[cond_lf])
18    hf = trapz(pxx[cond_hf], fxx[cond_hf])
19      # sum these up to get total power
20    #total_power = vlf + lf + hf
21
22    # find which frequency has the most power in each band
23    #peak_vlf = fxx[cond_vlf][np.argmax(pxx[cond_vlf])]
24    peak_lf = fxx[cond_lf][np.argmax(pxx[cond_lf])]
25    peak_hf = fxx[cond_hf][np.argmax(pxx[cond_hf])]
26
27    # fraction of lf and hf
28    lf_nu = 100 * lf / (lf + hf)
29    hf_nu = 100 * hf / (lf + hf)
30    
31    results = {}
32    #results['Power VLF (ms2)'] = vlf
33    results['Power LF (ms2)'] = lf
34    results['Power HF (ms2)'] = hf   
35    #results['Power Total (ms2)'] = total_power
36    results['LF/HF'] = (lf/hf)
37    #results['Peak VLF (Hz)'] = peak_vlf
38    results['Peak LF (Hz)'] = peak_lf
39    results['Peak HF (Hz)'] = peak_hf
40
41    results['Fraction LF (nu)'] = lf_nu
42    results['Fraction HF (nu)'] = hf_nu
43    return results, fxx, pxx
44
45print("Frequency domain metrics:")
46results_2, fxx_2, pxx_2 = frequency_domain(f_sci(xx))
47
48for k, v in results_2.items():
49    print("- %s: %.2f" % (k, v))

私の信号

python
1print(xx)
2#[ 20.5    20.525  20.55  ... 379.925 379.95  379.975]
3print(f_sci(xx))
4#[5.61056341 5.58633306 5.57034141 ... 5.93681271 5.93674881 5.93688724]
5

発生している問題・エラーメッセージ

Output exceeds the size limit. Open the full output data in a text editor
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_13212/220404614.py in <module>
      1 print("Frequency domain metrics:")
----> 2 results_2, fxx_2, pxx_2 = frequency_domain(f_sci(xx))
      3 
      4 for k, v in results_2.items():
      5     print("- %s: %.2f" % (k, v))

~\AppData\Local\Temp/ipykernel_13212/92220601.py in frequency_domain(rri, fs)
     22     # find which frequency has the most power in each band
     23     #peak_vlf = fxx[cond_vlf][np.argmax(pxx[cond_vlf])]
---> 24     peak_lf = fxx[cond_lf][np.argmax(pxx[cond_lf])]
     25     peak_hf = fxx[cond_hf][np.argmax(pxx[cond_hf])]
     26 

<__array_function__ internals> in argmax(*args, **kwargs)

c:\Users\ryous\anacon\lib\site-packages\numpy\core\fromnumeric.py in argmax(a, axis, out)
   1191 
   1192     """
-> 1193     return _wrapfunc(a, 'argmax', axis=axis, out=out)
   1194 
   1195 
...
---> 58         return bound(*args, **kwds)
     59     except TypeError:
     60         # A TypeError occurs if the object does have such a method in its

ValueError: attempt to get argmax of an empty sequence

行動規範の内容に同意します

回答1件

ベストアンサー

原因

LFの帯域値が存在しなかったため，抽出条件cond_lf = (fxx >= 0.04) & (fxx < 0.15)でTrueになっているindexが0個となり，pxx[cond_lf]として抽出条件を適用した抽出結果の配列サイズが0，結果np.argmax([])という処理が行われ，ValueError: attempt to get argmax of an empty sequenceとなったようです．

このエラーを治すには，スペクトル周波数の配列fxxにおいてLFの帯域値が必要ということになります．配列fxxの中身を確認しながらデバッグするか，周波数スペクトラムのグラフを表示しながらデバッグする必要があります．

原因の再現/解析

そちらで解析のために持っているデータの配列サイズは不明だったので，適当な整数Nを用いて弊環境でfxx, pxx = signal.welch(x = np.random.randint(0, 100, N), fs = 40.0)として実行すると，Nの値に依存せず配列fxxの長さlen(fxx)は必ず129になりました．signal.welch()の引数npersegのデフォルト値256の半分+1の値がスペクトル周波数fxxの配列長len(fxx)になるからです．

ところでこの配列fxxの中を見ると，0番目から順に[0, 0.15625, 0.3125, 0.46875, 0.625, 0.78125...となっており，コードで定義した抽出条件cond_lf = (fxx >= 0.04) & (fxx < 0.15)に合致する値は存在しません．

改善案

したがって抽出条件に合う値を存在させるために，サンプリング周波数fsを小さくするよう変更するか，signal.welch()の引数npersegの値を大きくするよう変更してください．前者は物理的な定数ですので，実際変更可能なのは後者かと思われます．

セグメント長npersegは最大で入力データの長さlen(f_sci(xx))までの値を設定することができます．Kaggleの方のコードではfs=4.0,nperseg=256であるので，相対的に考えてfs=40.0と10倍にしていることからnpersegも10倍にして2560にすると同じスペクトル密度になるでしょう．

投稿2022/10/05 15:01

編集2022/10/05 22:34

ps_aux_grep

総合スコア1579

ps_aux_grep

2022/10/05 15:50 編集

Kaggleのデータとあなたのデータのうち，fsに差があり上のようなエラーに出会うことになっています．ナイキスト周波数までの周波数を，セグメント長npersegの半分の長さを持つ配列に押し込めたため，fxxの解像度が相対的に低くなった．ということです．あまりにもnpersegを大きくしすぎても，精度に影響が出ると思われるのでパワースペクトラムを描きながら変更することをお勧めします．