フーリエ変換のsin,cos成分について
30 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
初学者です。初歩的な質問失礼いたします。
フーリエ変換のさいにsin波とcos波に分解されると思いますが,fft関数を使用した後に、それぞれの成分ごとのスペクトルパワーを取得したいです。方法を教えてくださいますでしょうか。
また加速度計などの解析で使用されるharmonicRatioは偶関数成分(cos)と奇関数成分(sin)の比率と考えておりますが,単純にパワーの割り算でよろしいのでしょうか。
お詳しい方がいらっしゃいましたらお願いいたします。
0 comentarios
Respuestas (2)
Hernia Baby
el 17 de Mzo. de 2022
時系列信号でものを言います。
フーリエ変換は時間領域を周波数領域に写像変換します。
例でみてみましょう。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)について
だと思います。
2 comentarios
Hernia Baby
el 17 de Mzo. de 2022
前半の長々としたものは、僕のQiita記事からとってます。
Hernia Baby
el 18 de Mzo. de 2022
■上記の例では具体的にsin,cos成分はどれになるのでしょうか
cos成分とsin成分は実部と虚部です。
パワースペクトルは絶対値をどうこうするので r の部分に相当します。
■この時左右対象なのでHRはほぼ1になるという解釈でよろしいでしょうか
基本周波数に対して奇数次成分と偶数次成分の比較から歪度をみるんじゃなかったでしたっけ?
あまり詳しくはないのではっきりとした回答はできませんが…
Ver también
Categorías
Más información sobre フーリエ解析とフィルター処理 en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!