質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.45%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

1回答

613閲覧

python 信号処理 エラー

ryosuke0313

総合スコア65

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2022/10/05 11:37

前提

信号に対して周波数解析を行っています。
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

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

ベストアンサー

原因

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
PondVillege

総合スコア1579

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

PondVillege

2022/10/05 15:50 編集

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

2022/10/06 04:18

回答ありがとうございます。実行できました。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.45%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問