9
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

MATLABで行うFFTその2 ~リーク誤差と窓関数~

Last updated at Posted at 2022-01-28

はじめに

今回は リーク誤差 と対策方法である の紹介をしていきます。

前回のあらすじ

前回の記事 では FFTで片側スペクトルの正規化 まで行いました。

しかし 狙った振幅よりも低い値 になってしまっていました。

今回はその原因と対策を書いていきます。

前回の復習

前回の復習から行いましょう。

Code
clear,clc,close all;
dt = 0.01;
L  = 512;
t = 0:dt:dt*(L-1);
Fs = 1/dt;

合成波について最大振幅は 10Hz で 0.5, 30Hz で 2 としています。

Code
y = .5*sin(2*pi*10*t) + 2*sin(2*pi*30*t);
plot(t,y)
xlim([0 t(end)])
xlabel '時間[sec]'
ylabel '信号'

figure_0.png

ここでFFTして正規化まで行います。

Code
% FFT
cy = fft(y);

% 周波数(サンプリング周波数の半分だけ)
f = Fs*(0:(L/2))/L;
f = f(1:end-1);

% 正規化
P = abs(cy(1:ceil(length(cy)/2)))./(length(y)/2);

% 図示
stem(f,P)
hold on

% 10Hz付近, 30Hz付近の行を検索
idx = knnsearch(f',[10;30]);
scatter(f(idx),P(idx),'o')
name = P(idx)';
text(f(idx),P(idx)+0.1,cellstr(num2str(name)));
ylim([0 1.7])
xlabel '周波数[Hz]'
ylabel 'P1(f)'
hold off

figure_1.png

数字でピークをみてみましょう

Code
Peak = [f(idx)', P(idx)'];
fprintf('ピーク周波数 %-.2f Hz \t 振幅 %-.2f\n\n',Peak(1,:),Peak(2,:))
Output
ピーク周波数 9.96 Hz 	 振幅 0.47

ピーク周波数 30.08 Hz  振幅 1.51

んー、ずれてますねえ(ニッコリ)

10Hzはともかく、30Hz付近は過小評価されています。

リーク誤差

このずれはリーク誤差によって生じています。

離散フーリエ変換は範囲を有限として、外側は範囲内の繰り返しとして仮定しています。

Code
subplot(1,2,1)
plot(t,y,'--')
hold on
scatter(t(1),y(1),'red','filled')
xlim([0 t(20)])
xlabel '時間[sec]'
ylabel '信号'
title 'はじめ'
subplot(1,2,2)
plot(t,y,'--')
hold on
scatter(t(end),y(end),'red','filled')
xlim([t(end-20) t(end)])
xlabel '時間[sec]'
ylabel '信号'
title 'おわり'

figure_2.png

一致してませんね…これがリーク誤差の原因です。

つまり最初と最後がぴったし同じで繰り返さないとダメなんですね。

そんなん実践で毎回使えるやついる~?

image_0.png

対策:窓関数

じゃあどうするのか?偉い人は考えました。

最初と最後を0に収束させればよくね?!

そこで出てくるのが窓関数ってやつです。

ここではとりあえず迷ったら使えばいいと噂のハニング窓を使います。

Code
figure
hw = hann(length(y))';
plot(t,y,':','Color',[.7 .7 .7],'LineWidth',0.1); hold on
plot(t,hw.*y,'r','LineWidth',0.8);
xlim([0 t(end)])
xlabel '時間[sec]'
ylabel '信号'
legend({'通常','窓付き'})
hold off

figure_2.png

こんな感じです

窓付きをFFTしますか。

ただここでも正規化が必要になります。

窓関数を使用した場合の正規化

窓関数を使った場合の正規化は以下のようになります。

$g\left(n\right)$:窓関数

$\textrm{nfft}$:ブロックサイズ

$\left|S\right|$:片側スペクトル

c\;=\;\frac{\Sigma \left(g\left(n\right)\right)}{\textrm{nfft}}
S_{\textrm{norm}} =\frac{\left|S\right|}{c\cdot \left(\textrm{nfft}/2\right)}

コードで書いてみましょう。

Code
% FFT
cy2 = fft(hw.*y);

% 正規化
c = sum(hw)/length(hw);
P = abs(cy2(1:ceil(length(cy2)/2)))./(length(y)/2)./c;

% 図示
stem(f,P)
hold on

% 10Hz付近, 30Hz付近の行を検索
idx = knnsearch(f',[10;30]);
scatter(f(idx),P(idx),'o')
name = P(idx)';
text(f(idx),P(idx)+0.1,cellstr(num2str(name)));
ylim([0 2.5])
xlabel '周波数[Hz]'
ylabel 'P1(f)'
hold off

figure_3.png

Code
Peak = [f(idx)', P(idx)'];
fprintf('ピーク周波数 %-.2f Hz \t 振幅 %-.2f\n\n',Peak(1,:),Peak(2,:))
Output
ピーク周波数 9.96 Hz 	 振幅 0.49

ピーク周波数 30.08 Hz  振幅 1.80

完ぺきではないですが、リーク誤差はある程度改善されました。

おわりに

次回はオーバーラップでもやろうかな。

9
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
9
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?