パワースペクトルグラフ化 x軸とy軸の要素数が合わない

Question

前提・実現したいこと

Pythonで動画から脈拍測定を目的としたシステムを作っています。
その中で、パワースペクトルグラフ化をしたいのですが、実装中に以下のエラーメッセージが発生しました。

初心者で分からないことが多いのですが、よろしくお願い致します。

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

Traceback (most recent call last):
  File "HR_Video.py", line 116, in <module>
    plt.plot(result[1][1:-1], result[0][1:])
  File "C:\Users\〇〇\Anaconda3\lib\site-packages\matplotlib\pyplot.py", line 2795, in plot
    is not None else {}), **kwargs)
  File "C:\Users\〇〇\Anaconda3\lib\site-packages\matplotlib\axes\_axes.py", line 1666, in plot
    lines = [*self._get_lines(*args, data=data, **kwargs)]
  File "C:\Users\〇〇\Anaconda3\lib\site-packages\matplotlib\axes\_base.py", line 225, in __call__
    yield from self._plot_args(this, kwargs)
  File "C:\Users\〇〇\Anaconda3\lib\site-packages\matplotlib\axes\_base.py", line 391, in _plot_args
    x, y = self._xy_from_xy(x, y)
  File "C:\Users\〇〇\Anaconda3\lib\site-packages\matplotlib\axes\_base.py", line 270, in _xy_from_xy
    "have shapes {} and {}".format(x.shape, y.shape))
ValueError: x and y must have same first dimension, but have shapes (448,) and (449,)

```Python

該当のソースコード

#パワースペクトルを求める関数を定義
def signal_fft (data,fs):
n = len(data)-1
y = fftpack.fft(data)/n
y = y[0:round(n/2)]
power = np.abs(y)
# power = 10*np.log10(power)
f = np.arange(0,fs/2,fs/n)
return power , f

######################################################追加部分です
#15sの窓でパワースペクトル算出　　　　　　　　　　　　　　　　　　　　　　　　　
time_interval = 15
loop_len = list(range(round(len(img_mean)/(fstime_interval))))
fft_result = []
for n,i in enumerate(loop_len):
if n ==0:
img = img_mean[:time_intervalfs]
fft_result.append(signal_fft(img,fs))
else :
img = img_mean[time_intervalfsn : time_intervalfs(n+1)]
fft_result.append(signal_fft(img,fs))

パワースペクトルグラフ化

#print(fft_result)
time = np.arange(0, len(fft_result[0][0]), 1)
for i, result in enumerate(fft_result):
plt.plot(result[1][1:-1], result[0][1:])
plt.xlabel("Frequency(Hz)",size=12)
plt.ylabel("Power(dB)",size=12)
fig_path = "./Cropped_Output_figs/{0}.png".format(i) #ディレクトリ変更可能
if os.path.isfile(fig_path):
os.remove(fig_path)
plt.savefig(fig_path)
plt.clf()
plt.show()

avg_result = []
N = 5
b = np.ones(N) / N
for i, result in enumerate(fft_result):
avg_result.append([np.convolve(result[0], b, mode="same").tolist(), result[1]])
plt.plot(result[1][1:-1], avg_result[i][0][1:])
plt.xlabel("Frequency(Hz)",size=12)
plt.ylabel("Power(dB)",size=12)
fig_path = "./Cropped_Output_figs/avg_{0}.png".format(i) #ディレクトリ変更可能
if os.path.isfile(fig_path):
os.remove(fig_path)
plt.savefig(fig_path)
plt.clf()

### 試したこと

plt.plot(result[1][1:-1], result[0][1:])
の部分が問題だと思ったので、
result[1][1:-1], result[0][1:-1]
result[1][1:-1], result[0][1:0]
のようにしてみると、何もプロットされない状態でグラフが表示されました。

### 補足情報（FW/ツールのバージョンなど）

Comments

[meg_]

result[1][1:-1], result[0][1:-1]にはプロットしたい値が入っていますか？

[jbpb0]

signal_fft()が返すpowerとfの数が違うのがエラーの原因なら、 round()が悪さしてるのかもしれません

下記を実行したら分かりますが、Pythonのround()は普通の四捨五入ではありませんので、気を付けてください
round(2.5)
round(3.5)

詳細は、下記を参照してください
https://note.nkmk.me/python-round-decimal-quantize/

[jbpb0]

fft_resultには何が入っているのですか?
signal_fft()の実行結果ですか?

[jbpb0]

長い1次元データ img_mean を fs*time_interval 個ずつに分割して、スペクトルを計算しているようですね
挙げられているコードだと、fs*time_interval が奇数の場合と偶数の場合で、状況が変わります

>>> print(len(img_mean))
90
>>> print(fs)
3
>>> print(time_interval)
15
の場合は、fs*time_interval=45 と奇数で、下記のように周波数とスペクトルの数が一致します
>>> print(len(fft_result[0][0]))
22
>>> print(len(fft_result[0][1]))
22

>>> print(len(img_mean))
90
>>> print(fs)
2
>>> print(time_interval)
15
の場合は、fs*time_interval=30 と偶数で、下記のように周波数とスペクトルの数が一致せず、周波数の方が一つ多くなります
>>> print(len(fft_result[0][0]))
14
>>> print(len(fft_result[0][1]))
15

その結果、グラフを書くforループ内で、周波数(result[1])とスペクトル(result[0])の数が合ったり合わなかったりします
そうなる原因は、pythonのround()関数の仕様にあります(私の最初の書き込みを見てください)

[jbpb0]

対策を回答に書きましたが、img_meanが1次元データの場合でしか手元で検証してません
もしそうでない場合は、ちゃんと動かないかもしれません

Answer 1

nが整数だと分かっているので、n/2は整数か整数+0.5しかないので、下記のように0.1を足せば必ず四捨五入されます

# パワースペクトルを求める関数を定義
def signal_fft (data,fs):
    n = len(data)-1
    y = fftpack.fft(data)/n
    y = y[0:round(n/2+0.1)]
    ## Pythonのround()の仕様(偶数への丸め)を回避するため0.1を足す (nは必ず整数)
    power = np.abs(y)
    # power = 10*np.log10(power)
    f = np.arange(0,fs/2,fs/n)
    return power , f

そうすれば周波数とスペクトルの数が合うので、要素数の調整は不要になります

# パワースペクトルグラフ化
#import pprint
#pprint.pprint(fft_result)
time = np.arange(0, len(fft_result[0][0]), 1)
for i, result in enumerate(fft_result):
   plt.plot(result[1], result[0])
   plt.xlabel("Frequency(Hz)",size=12)
   plt.ylabel("Power(dB)",size=12)
   #fig_path = "./Cropped_Output_figs/{0}.png".format(i)          #ディレクトリ変更可能
   #if os.path.isfile(fig_path):
   #    os.remove(fig_path)
   #plt.savefig(fig_path)
   #plt.clf()
   plt.show()

もちろん、どんな場合でも0.1足せば良い、というわけではありません
あくまでも、今回の問題の場合は大丈夫、というだけです

Comments

[M_77]

ありがとうございます。
fs*time_intervalが偶数の時と奇数の時で周波数とスペクトルの数が一致しなくなってしまうことに驚きました。そんなところまでは考えが及ばなかったので、感謝です。round関数についても学べました。

初めて質問してみて不安だったのですが、非常に丁寧に回答して頂き、嬉しかったです。ありがとうございました。