MATLAB で学ぶ信号処理　伝達関数を逆FFTするとどうなる？

2021.09.08

はじめに

さっそく質問ですが、みなさん逆FFTって使っていますか？

音響関係の信号処理をしている方なら逆FFTは毎日のように使用している方が多いと思います。

しかし、振動関係では逆FFTを使用している方はあまりいないですよね。

本記事では、初心者向けにFFTと逆FFTの関係を説明し、伝達関数を逆FFTするとどうなるのかを説明致します。

 

FFTと逆FFTの関係

ここでは、FFTと逆FFT（IFFT:Inverse Fast Fourier Transform）の大まかな関係を説明します。

FFTをブロック線図形式で表すと図1のような関係で示すことができます。

図1　FFTのブロック線図

逆FFTのブロック線図を図2に示します。図1のFFTとは逆の変換をしていることがわかります。

図2　IFFTのブロック線図

 

MATLABでの検証

ここで、簡単にMATLABで検証してみましょう。

下記の波形1と波形2が一致していれば、図1と図2の関係は成立していることが確認できます。
波形1：任意の信号（ランダム波）
波形2：波形1をFFTした後に、逆FFTした信号

 

結果は図3のようになります。
波形1が青、波形2が赤で示してあります。上の図を拡大したものが中央の図です。波形1と波形2が一致していますね。
波形1と波形2の差分が一番下の図です。0になっていることがわかります。

図3　検証結果

 

clear all;clc;close all

N=2^14;
sig1=rand(1,N)-0.5;

SIG1=fft(sig1,N);
sig2=ifft(SIG1,N);

figure
subplot(311)
plot(sig1,'b')
hold on
plot(sig2,'r--')
xlim([1 N])

subplot(312)
plot(sig1,'b')
hold on
plot(sig2,'r--')
xlim([1 100])

subplot(313)
plot(sig1-sig2,'k')
ylim([-0.5 0.5])
xlim([1 N])

 

伝達関数を逆FFTすると？

結論から述べると、伝達関数を逆FFTするとインパルス応答になります。

まあ、ハンマリング試験の時にインパクトハンマーで叩いて応答を測定しているので、「インパルス応答になるでしょ！」って思っていた人もいるかもしれませんね。

ただ、伝達関数の測定方法はいろいろな方法がありますよね。振動関連だと下記があげられます。どんな測定方法で取得した伝達関数でも逆FFTするとインパルス応答に戻るんですよ。

  • インパクトハンマーによる打撃加振
  • 加振機を用いた加振試験
    →加振機への入力信号は様々（スイープ、ランダムなどなど）

 

MATLABでの検証

では、MATLABで検証してみましょう。
5自由度のバネマス系でFRFを作成し、そのFRFを逆FFTしてみます。

このとき、MATLABのような理論解では減衰行列Cを適切に設定しないと、逆FFTの検証がうまくいかないので注意してください。（うまくいかない理由については、ここで説明はしません。）

MATLABの検証結果が図4です。

図4　FRFを逆FFTして求めたインパルス応答

 

上記検証に用いたMATLABコードでわからない点がある方は、まずは下記リンクで勉強してみてください。FFTの結果がナイキスト周波数を基点に虚像となっている点についての理解が不足している可能性があります。

MATLABで学ぶ信号処理　 高速フーリエ変換（fast Fourier transform:fft）について 1
MATLABにはfft関数が用意されていて、簡単に誰でもfftができます。しかし、振動・騒音関係でMATLABのfftを使用する場合は注意が必要です。 騒音器や振動測定器で測定した時間波形をMATLABやPythonなどのプログラムで単純にfftをすると、騒音器や振動測定器で測定した”周波数分析結果”と一致しません。 今回は測定器の周波数分析結果と自作プログラムのFFT結果が異なる原因を説明したいと思います。また、騒音器や振動測定器と同様の結果が得られるようにFFTのプログラムを作成するにはどうすればいいかを紹介いたします。
mech-eng-uni.com
2021.03.20

 

ここで触れた内容はifftの入り口にすぎません。もっと深い使い方があるので、次からはもう少し深い内容について解説していきます。

ifftが使いこなせると、振動・騒音・音響関連で研究や仕事をしている方は、これまで解けなかった課題が簡単に解けるようになります。

MATLABコード

実行ファイル

clear all
close all
% % % % 初期設定
freq=[0.0001 1:1:2^15-1];
m_vec=ones(1,5)*0.1;
[M]=eval_Mmatrix(m_vec);
k_vec=ones(1,5)*10^4;
[K]=eval_Kmatrix(k_vec);
a=10^1;
b=10^-4;
C=a*M+b*K;
F1=zeros(length(m_vec),length(freq));
F1(1,:)=ones(1,length(freq));
% % % % 変位X（ノイズなし）
X1=eval_direct_x_2ndedition(M,K,C,F1,freq);
% % % % 理想的なFRF
FRF_x4_f1=X1(4,:)./F1(1,:);

% % % % % FRFの確認
% figure
% semilogx(freq, abs(FRF_x4_f1),'k')
% hold on
% semilogx(freq, real(FRF_x4_f1),'b')
% semilogx(freq, imag(FRF_x4_f1),'r')

% % % % FRFを鏡面
FRF=[FRF_x4_f1 0 conj(FRF_x4_f1(end:-1:2))];

sig=ifft(FRF,2^15*2);

figure
plot(real(sig))
xlabel('時間')
ylabel('変位')

 

functionファイル

function x=eval_direct_x_2ndedition(M,K,C,F,freq)
% % 2021.05.28 Fを修正
x=zeros(size(F,1),length(freq));
j=sqrt(-1);
for ii1=1:1:length(freq)
w=2*pi*freq(ii1);
x(:,ii1)=inv(K+j*w*C-w^2*M)*F(:,ii1);
end

end

function [K]=eval_Kmatrix(k_vec)
K=zeros(length(k_vec));
for ii1=1:length(k_vec)
if ii1==1 K(ii1,ii1)=k_vec(ii1);
else
K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);
end
end

end

function [M]=eval_Mmatrix(m_vec)
M=diag(m_vec);
end

コメント