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

12 views (last 30 days)
ryu
ryu on 31 Mar 2022
Answered: ryu on 5 Apr 2022
加速度などのデータをフーリエ変換した後に,sin成分とcos成分を取得する方法を教えてください。
  6 Comments
Hernia Baby
Hernia Baby on 3 Apr 2022
あー、なるほど!理解しました!
まずやりたいことに関する回答ですが、
ここでいうOdd/Even Harmonicsってたぶん違うと思います。
基本周波数があって、それの偶数倍か奇数倍のことを言ってると思います。
cos/sinは実部虚部に該当しますが、今回求めたいものとは関係ないです。
ちなみに左側右側に関係するのはナイキスト周波数になります。
複素共役の関係上、反対の周波数側にも同じようなスペクトルが生成されます。
なのでサンプリング周波数の半分しか見れませんよって感じのお話です。
それ以上みると折り返し雑音(エイリアシング)が生じます。
なので今回のお話とは別の現象ですね。
宣伝になりますが、こちらにざっくり書いてます↓
ryu
ryu on 3 Apr 2022
お二方丁寧な回答ありがとうございます。
やっと理解することができました。
Hernia Babyさんのおしゃっていただいたようにまずやりたいことが違っておりました。
Atsushi Uenoさんのご回答のおかげで間違いに気づくことができました。
最後に1点質問させていただきたいのですが、フーリエ変換後の周波数分解能の関係で、ぴったり整数にはならないと思います。下の図は実際の周波数分解です。この場合、偶数波,奇数波と判断するにはそれぞれ四捨五入したものを扱ってもいいと思われますか。少し細かい話になってしまい申し訳ございませんが、ご意見いただけると嬉しいです。

Sign in to comment.

Accepted Answer

Hernia Baby
Hernia Baby on 4 Apr 2022
Edited: Hernia Baby on 4 Apr 2022
やりたいこともわかりましたので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")

More Answers (1)

ryu
ryu on 5 Apr 2022
大変わかりやすい手順の説明ありがとうございました。
無事にコードを作成することができました。
今後ともよろしくお願い致します

Community Treasure Hunt

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

Start Hunting!