フーリエ変換の各成分取得

加速度などのデータをフーリエ変換した後に,sin成分とcos成分を取得する方法を教えてください。

[Atsushi Ueno]

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

＞FFTをかけるとどうなるか？複素数に変換されます。

＞上記の例では具体的にsin,cos成分はどれになるのでしょうか

＞cos成分とsin成分は実部と虚部です。

既に同じ質問/回答済なのに、再度同じ質問をされているのが不可解です。前回質問の回答の何が不明なのか明確にして頂き、また質問内容を他の言葉や図で説明して頂けないでしょうか？

[ryu]

説明不足で申し訳ございません。

下の図はある加速度データをフーリエ変換し,パワースペクトルを求めたものです。

この場合左半分がcos,右半分がsin成分という解釈でよろしいでしょうか。

また少し専門的な話になりますが,加速度データを使用した滑らかさの指標であるHarmonic Ratio(http://www.reha.kobegakuin.ac.jp/~rehgakai/journal/files/no4-2/ronbun03.pdf)を求めたいと思っています。これはcos成分とsin成分の比率(それぞれのパワーを積分したものの比率)と解釈しております。フーリエ変換後にパワースペクトルを求めた際は下の図のように、大体左右対称(sin,cosの成分が同等)になる印象ですが,計算理論上そういうものなのでしょうか。もしそうだとするとHarmonic Ratioはいつも1になってしまうのではないのでしょうか。

長くなり申し訳ございませんが,回答していただけると幸いです。

[Atsushi Ueno]

3. フーリエ変換 (やる夫で学ぶディジタル信号処理) (tohoku.ac.jp) (3.56)の辺りより抜粋

　

>この場合左半分がcos,右半分がsin成分という解釈でよろしいでしょうか。

多分正しい解釈をお持ちだと思いますが、厳密な言葉の定義を要求すればNOです。

例えば上記波形なら「左半分がの成分、右半分がの成分」という解釈です。

>これ(Harmonic Ratio)はcos成分とsin成分の比率(それぞれのパワーを積分したものの比率)と解釈しております。

これ(Harmonic Ratio)は偶数歩と奇数歩の波形のパワー和の比と解釈しました。同じ事でしょうか? ちょっと読みましたが、下記の様に書いてあります。

HRは歩行動作が主として１歩を一周期とする２周期（１歩行周期）の加速度変化で成立していることに注目した指標である。

http://www.reha.kobegakuin.ac.jp/~rehgakai/journal/files/no4-2/ronbun03.pdf

偶/奇数歩の波形が同じならHarmonic Ratioが1で、左右に違いがあれば1から外れるという事ですよね。

1歩行周期(2歩)ではなく1周期(1歩)のHarmonic Ratioが1に近くなるのは自明ではないでしょうか？

[Atsushi Ueno]

つまりここで出てくる「sin成分」「cos成分」が何を意味するのか厳密に合わせないと、回答に辿り着かないという事なのでしょう。

[Hernia Baby]

あー、なるほど！理解しました！

まずやりたいことに関する回答ですが、

ここでいうOdd/Even Harmonicsってたぶん違うと思います。

基本周波数があって、それの偶数倍か奇数倍のことを言ってると思います。

cos/sinは実部虚部に該当しますが、今回求めたいものとは関係ないです。

ちなみに左側右側に関係するのはナイキスト周波数になります。

複素共役の関係上、反対の周波数側にも同じようなスペクトルが生成されます。

なのでサンプリング周波数の半分しか見れませんよって感じのお話です。

それ以上みると折り返し雑音(エイリアシング)が生じます。

なので今回のお話とは別の現象ですね。

宣伝になりますが、こちらにざっくり書いてます↓

MATLABで行うFFT ~時系列データ生成から片側スペクトルの正規化まで~

[ryu]

お二方丁寧な回答ありがとうございます。

やっと理解することができました。

Hernia Babyさんのおしゃっていただいたようにまずやりたいことが違っておりました。

Atsushi Uenoさんのご回答のおかげで間違いに気づくことができました。

最後に１点質問させていただきたいのですが、フーリエ変換後の周波数分解能の関係で、ぴったり整数にはならないと思います。下の図は実際の周波数分解です。この場合、偶数波,奇数波と判断するにはそれぞれ四捨五入したものを扱ってもいいと思われますか。少し細かい話になってしまい申し訳ございませんが、ご意見いただけると嬉しいです。

Answer 1

[Hernia Baby]

やりたいこともわかりましたのでHarmonic ratioを求めるための手順を書いていこうかなと思います。

まず準備

clc,clear;
rng('default');
Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
t = t';

次数を1~50次まで用意します

基本周波数は 5 Hz なので 5*偶数次 Hz と 5*奇数次 Hz を作ります

n = 1:50;
idx = mod(n,2)==0;

f_base = 5;
even = f_base*n(idx);
odd  = f_base*n(~idx);

テキトーな係数を作ります

A_even = randi([0 2],size(even));
A_odd  = randi([0 3],size(odd));

正弦波どうしを足し合わせます

x_even = A_even.*sin(2*pi*even.*t);
x_odd  = A_odd.*sin(2*pi*odd.*t);
y = sum([x_even, x_odd],2) + rand(size(t));

figure
plot(t,y,'Color',[1 1 1]*.6)

平均パワーを最大振幅に変換して一度見てみましょう

[p,f] = pspectrum(y,Fs);
coefs = sqrt(2.*p);

figure
plot(f,coefs,'Color',[1 1 1]*.6)
xlim([0 250])

本当は基本周波数がわからないと思うので、最大を求めて偶奇で分けるのがいいのかもしれません

今回は各次元周波数に最も近い周波数を選んで比較しています

[f_even,~] = dsearchn(f,even');
[f_odd,~] = dsearchn(f,odd');

HR = sum(coefs(f_even))/sum(coefs(f_odd))

HR = 0.8946

最大を求めるのはタスクの局所極値を使用すれば楽にできます

ちなみに250Hzより高周波な部分はノイズです

% 局所的最大値の検出
maxIndices = islocalmax(coefs,"SamplePoints",f);

% 結果の表示
clf
plot(f,coefs,"Color",[1 1 1]*.6,"DisplayName","入力データ")
hold on

% 局所的最大値のプロット
plot(f(maxIndices),coefs(maxIndices),"o","Color",'r',...
    "MarkerFaceColor",'r',"DisplayName","局所的最大値")
title("極値の数: " + nnz(maxIndices))
hold off
legend
xlabel("f")

Answer 2

[ryu]

大変わかりやすい手順の説明ありがとうございました。

無事にコードを作成することができました。

今後ともよろしくお願い致します