FFTした後、周波数​と時間ごとのパワー強​度のデータが欲しい。

FFTをして　Ｘ軸：周波数 , Y軸 : パワーのプロットはできたのですが、

1秒後毎の、周波数毎のパワー強度の求めからを教えてください。 よろしくお願いします。

サンプリングレート：　250

サンプリング時間：　　0.004sec

切り出すデータ時間：　1sec

オーバーラップ：　0

トータル時間：　6min

周波数：　0-125 Hz

Answer 1

[Shunichi Kusano]

短時間フーリエ変換のことかと思います。関数がありますので、任意のパラメータ設定してお使いください。

https://jp.mathworks.com/help/signal/ref/stft.html

Answer 2

[Hernia Baby]

二つアプローチがありますが、やってることはだいたい同じです。

■信号作成

サンプルに10Hzの信号を作ります

clc,clear;
sec = 6;
Fs = 1/4e-3;
t = (0:1/Fs:sec)';
f1 = 10;
x = sin(2*pi*f1.*t);

figure
hold on
plot(t,x,'Color',[.5 .5 .5])
plot(t,ones(size(t)).*(max(x)/sqrt(2)),'--r')
xlim([0 .2])
xlabel 'time'

■stft

0.5秒分ごとに区切り、リーク誤差の軽減のためにハニングウィンドウを使いました

L = 125;
w = hann(L,"periodic");
[s,f,tt]=stft(x,Fs,'Window',w,...
    'Overlaplength',0,...
    'FrequencyRange','onesided');

窓の正規化をして図示します

figure
c = sum(w)/length(w);
waterfall(f,tt,(((abs(s)/(L/2))/c/sqrt(2)).^2)')
xlim([0 125])
ylim([0 sec])
xlabel('frequency')
ylabel('time')

■pspectrum

やっていることはstftと同じです

違うのはブロックサイズの正規化はしなくていい点ですかね

窓の正規化は必要です

[p,fp,tp] = pspectrum(x,Fs,'spectrogram','TimeResolution',.5,'OverlapPercent',0,'Leakage',0.75);
figure
waterfall(fp,tp,p'/c^2/4)
xlim([0 125])
ylim([0 sec])
xlabel('frequency')
ylabel('time')

■結論

1secのデータを中央値である0.5secだと考えれば、

ブロックサイズをサンプリング周波数として取り出すことが可能です

今回だと見栄え上 0.5sec ずつで刻んだので、その刻みごとのデータがとれるというわけですね

[友一 安藤]

ありがとうございます。

もらった情報を元に、試してみましたが、上手く表示できません。

x軸、Y軸の設定がおかしいのでしょうか？

[Hernia Baby]

ウィンドウサイズの設定かもしれません 私の書いたコードでいうLはどのようにしましたか？

[友一 安藤]

L は　時間ではなく、ポイント数にしていたため、　上手く表示出来なかったようです。

アドバイス、ありがとうございました。