時系列信号でものを言います。
フーリエ変換は時間領域を周波数領域に写像変換します。
例でみてみましょう。10Hzと30Hzの正弦波を合成します。
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);
figure
plot(t,y)
xlim([0 t(end)])
xlabel '時間[sec]'
ylabel '信号'
FFTをかけるとどうなるか？複素数に変換されます。
cy = fft(y);
n = 100;
cys = cy([n end-n+2]);
tt = linspace(0,2*pi,100);
% 図示
figure
plot(sin(tt),cos(tt),'--')
axis([-2,2,-2,2])
hold on
plot(cys./(abs(cys)),'o:')
xlabel 'Real'
ylabel 'Imag'
xline(0)
yline(0)
hold off
上記のように複素共役の状態で出されます。
パワースペクトルとはフーリエスペクトルを 
とすると 
です。
とすると です。ここで注意すべきは実効値の二乗でパワースペクトルは算出されるということです。
つまり振幅を A とすると 
となります。
となります。 正規化して以下に示します。
f = Fs*(0:(L/2))/L;
f = f(1:end-1);
P = abs(cy(1:ceil(length(cy)/2)))./(length(y)/2);
P = P.^2;
% 図示
figure
stem(f,P)
hold on
idx = knnsearch(f',[10;30]);
scatter(f(idx),P(idx),'o')
name = round(P(idx),1)';
text(f(idx),P(idx)+0.1,cellstr(num2str(name)));
ylim([0 2.5])
xlabel '周波数[Hz]'
ylabel 'P(f)/Hz'
リーク誤差が発生していますが、上記のように求めることができました。
----
次はもっと簡単な方法を紹介します
ぶっちゃけ計算は1行で終わります。めちゃくちゃすごいです。
[p,f] = pspectrum(y,Fs);
% 図示
figure
plot(f,p,'--');
xlabel '周波数[Hz]'
ylabel 'P1(f)/Hz'
ylim([0 2.5])
hold on
idx = knnsearch(f,[10;30]);
scatter(f(idx),p(idx),'o')
name = round(p(idx),1);
text(f(idx),p(idx)+0.1,cellstr(num2str(name)));
hold off
■Harmonic Ratio(HR)について
だと思います。
