フーリエ変換のsin,cos成分について

30 views (last 30 days)
ryu
ryu on 16 Mar 2022
Commented: ryu on 24 Mar 2022
初学者です。初歩的な質問失礼いたします。
フーリエ変換のさいにsin波とcos波に分解されると思いますが,fft関数を使用した後に、それぞれの成分ごとのスペクトルパワーを取得したいです。方法を教えてくださいますでしょうか。
また加速度計などの解析で使用されるharmonicRatioは偶関数成分(cos)と奇関数成分(sin)の比率と考えておりますが,単純にパワーの割り算でよろしいのでしょうか。
お詳しい方がいらっしゃいましたらお願いいたします。

Answers (2)

Hernia Baby
Hernia Baby on 17 Mar 2022
時系列信号でものを言います。
フーリエ変換は時間領域を周波数領域に写像変換します。
例でみてみましょう。10Hzと30Hzの正弦波を合成します。
fft を使う
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'
リーク誤差が発生していますが、上記のように求めることができました。
----
次はもっと簡単な方法を紹介します
pspectrum を使う
ぶっちゃけ計算は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)について
  だと思います。
 こちらの 式 (2) をご参照ください。
 とはいえ、harmonicRatio は相関係数のような形をしており、少し違う気がしますが…。
  2 Comments
Hernia Baby
Hernia Baby on 17 Mar 2022
前半の長々としたものは、僕のQiita記事からとってます。
Hernia Baby
Hernia Baby on 18 Mar 2022
■上記の例では具体的にsin,cos成分はどれになるのでしょうか
 cos成分とsin成分は実部虚部です。
 パワースペクトルは絶対値をどうこうするので r の部分に相当します。
■この時左右対象なのでHRはほぼ1になるという解釈でよろしいでしょうか
 基本周波数に対して奇数次成分と偶数次成分の比較から歪度をみるんじゃなかったでしたっけ?
 あまり詳しくはないのではっきりとした回答はできませんが…
 基本周波数 の考えでいくと 0 だと思っています。

Sign in to comment.


ryu
ryu on 17 Mar 2022
大変丁寧な回答ありがとうございます。
pspectrumは知りませんでした。
上記の例では具体的にsin,cos成分はどれになるのでしょうか
複素共役の図で言うとx軸より左がsin,右がsinという解釈でよろしいでしょうか
またこの時左右対象なのでHRはほぼ1になるという解釈でよろしいでしょうか
  2 Comments
Hernia Baby
Hernia Baby on 18 Mar 2022
自分のコメントに記載しました
ryu
ryu on 24 Mar 2022
大変丁寧に説明して頂きありがとうございました!

Sign in to comment.

Categories

Find more on フーリエ解析とフィルター処理 in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!