パワーの正規化が必要だと考えられます。
単にFFTして2乗にしただけでは意味をなしません。
ですのでブロックサイズで割って2倍する必要があります。
dt = 0.01; % 時間分解能
L = 512; % ブロックサイズ
t = 0:dt:dt*(L-1); % 時間
Fs = 1/dt; % サンプリング周波数
y = .5*sin(2*pi*10*t) + 2*sin(2*pi*30*t); % 波形
plot(t,y,'Color',[.3 .3 .3])
xlim([0 t(end)])
xlabel '時間[sec]'
ylabel '信号'
これをFFTしてパワー値にします
plot(abs(fft(y)).^2,'Color',[.3 .3 .3])
xlim([0 L])
xline(L/2,':')
パーシバルの定理は入力信号の振幅二乗和とパワースペクトルの平均が一致するというようなものです。
式は で表されます。
つまりブロックサイズで割ってあげる必要があります。
また折り返しの部分も足す必要があります。ですので2倍にします。
横軸は以下の通りです。
f = Fs*(0:ceil(L/2)-1)/L;
縦軸配下の通りです。
P = abs(fft(y)).^2;
P = P(1:ceil(L/2))./ceil(L/2); % 片側スペクトル
plot(f,P,'Color',[.3 .3 .3])
xlabel '周波数[Hz]'
ylabel '片側スペクトル'
詳しくは手前味噌ですがこちらご覧ください。
■懸念点
ここからは余談なのですが、
音なのでdBで表すと思うのですが、パワー値のままでいいのでしょうか?