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

27 views (last 30 days)
友一 安藤
友一 安藤 on 22 Mar 2022
Commented: 友一 安藤 on 24 Mar 2022
FFTをして X軸:周波数 , Y軸 : パワーのプロットはできたのですが、
1秒後毎の、周波数毎のパワー強度の求めからを教えてください。 よろしくお願いします。
サンプリングレート: 250
サンプリング時間:  0.004sec
切り出すデータ時間: 1sec
オーバーラップ: 0
トータル時間: 6min
周波数: 0-125 Hz

Answers (2)

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

Hernia Baby
Hernia Baby on 22 Mar 2022
Edited: Hernia Baby on 22 Mar 2022
二つアプローチがありますが、やってることはだいたい同じです。
■信号作成
サンプルに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')
やっていることは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 ずつで刻んだので、その刻みごとのデータがとれるというわけですね
  3 Comments
友一 安藤
友一 安藤 on 24 Mar 2022
L は 時間ではなく、ポイント数にしていたため、 上手く表示出来なかったようです。
アドバイス、ありがとうございました。

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!