はじめに
前回の記事で逆FFTについて説明しました。
前回の記事を応用すると、周波数次元上でフィルタリングをして、時刻歴のデータに戻すということができるかと思います。
簡単な話だと、図1のように周波数次元の信号（Signal）にフィルター（Filter）を乗じて、逆FFTすれば、フィルタリング処理した信号（filtered signal）が生成できますよね。
図1
しかし、時刻歴の信号をそのままフィルタリング処理する方が一般的です。もし、時刻歴信号をFFTして、フィルターを乗じてから逆FFTしてしまうと、窓のつなぎ目がちぐはぐになってしまい、窓を跨いだ信号間では振幅＆位相が所望の値になりません。もっと具体的に述べると、FFTの窓をサンプリング周波数分のデータ数（1s分のデータ数）で実施した場合、0～1s間は所望の信号処理ができていますが、窓のつなぎ目である1sを跨ぐ0.5～1.5s間ではちぐはぐな振幅＆位相関係になってしまいます。
この問題を解決できるのが、時刻歴信号の信号処理です。時刻歴信号のフィルタリング手法としては、FIRフィルタやIIRフィルタ、バターワースなどが良く知られていますね。
上記手法は古典的であり、1980年代のようなロースペックなPCでもフィルタリング処理できる手法として活躍していました。高速で演算処理が可能な現代においては、各周波数における振幅や位相を変更したいというニーズがありますが、これらに対応するには別の手法が必要です。
それが、逆FFTと畳込み積分です。ちなみに参考文献は下記です。
ここまでの説明を読んで、「あー！どうせあのことを説明するんでしょ」と予測できている方は、おそらく少数かと思います。ご存じなかった方は、目から鱗の手法かと思いますので最後までお付き合いください。これができると振動や音響の信号処理で出来ることが圧倒的に増えます。
逆FFTについて
前回の記事で逆FFTについて説明しました。結論としては、「FRFを逆FFTするとインパルス応答になるよ」ってことです。
上記を少し発展させると、各周波数で振幅＆位相を操作するためのフィルターを作成して、そのフィルターを逆FFTすると、フィルター関数のインパルス応答が求まります。
少し説明がわかりにくいと思いますので、MATLABの結果を示しながら説明します。図2が自作したフィルター関数です。振幅を周波数軸上でサイン波のように変化させて、位相をコサイン波のように変化させるようなフィルター関数です。
図2 フィルター関数
これを逆FFTすると図3のようになります。
図3 フィルター関数を逆FFTして算出したインパルス応答
インパルス応答のような信号が求まりましたね。
インパルス応答を用いたフィルタリング手法
冒頭で説明した、逆FFTと畳込み積分を用いたフィルタリング手法（インパルス応答を用いたフィルタリング手法）を図4に示します。
図4 インパルス応答を用いたフィルタリング手法
図4に示すように、フィルター関数やFRFを逆FFTしてインパルス応答を求めて、インパルス応答と処理したい信号とを畳み込み積分することで、フィルタリング処理した信号を生成することができます。
図5 フィルタリング処理した信号
はたして、このフィルタリング処理した信号が図2で示したフィルター関数と一致するような操作ができているのかが重要ですよね。
そこで、オリジナルの信号をレファレンスとし、フィルタリング処理した信号をレスポンスとしてFRFを求めます。このFRFが図2と一緒の結果になっていれば、所望のフィルタリング操作ができていることになります。
図6 フィルター関数とFRFの比較
結果が一致していますね。
従って、所望のフィルタリング処理ができていることがわかります。
これができるようになると、音場制御やノイズキャンセリング、制振などなどめちゃめちゃできることが増えます。
また、FRFを測定して、任意の入力波形（時刻歴）が作用したときの応答波形（時刻歴）を簡単に求めることができます。
MATLABコード
実行コード
clear all
close all
clc
% % % % 初期設定
Nq=1024;
freq=0:1:Nq;
% % % % フィルターを作ってみる
phy=linspace(0,10*pi,Nq+1);
% phy=linspace(0,0,Nq+1);
Ampli=0.3*sin(phy)+3;
Phase=pi*cos(phy);
Filter_Nq=Ampli.*exp(i*Phase);
% % % % FRFを鏡面部分の作成
Filter_fs=[Filter_Nq(1:end-1) 1 conj(Filter_Nq(end-1:-1:2))];
% % % フィルター関数の周波数特性の確認
figure
subplot(211)
plot(abs(Filter_fs))
ylim([0 4])
xlim([0 Nq])
title('振幅')
xlabel('周波数 [Hz]')
subplot(212)
plot(angle(Filter_fs))
xlim([0 Nq])
title('位相')
xlabel('周波数 [Hz]')
% % % フィルター関数のインパルス応答化
filter_imp=real( ifft(Filter_fs,Nq*2) );
filter_imp=[filter_imp(end/2:end) filter_imp(1:end/2-1)];
% % % インパルス応答の確認
figure
plot(linspace(0,1,Nq*2),filter_imp)
% % % オリジナルの信号の作成
rand_sig=10*( rand(1,Nq*10)-0.5 );
% % % FFT用の窓関数の作成
[win,fft_coe]=get_window(Nq*2,'hanning');
% % % インパルス応答を用いたフィルタリング
filtered_sig=conv(rand_sig,filter_imp);
% % % フィルタリング語の信号はデータ数がオリジナルよりも長くなってしまうので、オリジナルの信号のデータ長を比較しやすいように調整
p_fil=zeros(1,length(filter_imp));
p_fil(floor(length(p_fil)/2)+2)=1;
normal_sig=conv(rand_sig,p_fil);
figure
plot(filtered_sig,'r')
hold on
plot(normal_sig,'b')
% % % % フィルタリングができているかをFRFを求めて確認
% % 50回平均してみましょうか
N=50;
Normal_sig=zeros(N,Nq*2);
Filtered_sig=zeros(N,Nq*2);
% % % % FFT
for ii1=1:N
Normal_sig(ii1,:)=fft(win.*normal_sig(2000+Nq/32*(ii1-1):2000+Nq/32*(ii1-1)+Nq*2-1),Nq*2);
Filtered_sig(ii1,:)=fft(win.*filtered_sig(2000+Nq/32*(ii1-1):2000+Nq/32*(ii1-1)+Nq*2-1),Nq*2);
end
% % % % FRF
ratio=mean(Filtered_sig./Normal_sig,1);
figure
subplot(2,1,1)
plot(freq,abs(Filter_Nq))
hold on
plot(freq,abs(ratio(1:length(freq))),'r--')
ylim([0 4])
title('振幅')
subplot(2,1,2)
plot(freq,angle(Filter_Nq))
hold on
plot(freq,angle(ratio(1:length(freq))),'r--')
ylim([-pi pi])
title('位相')
functionファイル
function [win,fft_coe]=get_window(N,window)
t=linspace(0,1,N);
switch window
case 'rectangular'
win = ones(1,N);
fft_coe=[1 1];
case 'hanning'
win = 0.5*ones(1,N) - 0.5*cos(2.0*pi*t);
fft_coe=[2 8/3];
case 'hamming'
win = 0.54*ones(1,N) - 0.46*cos(2.0*pi*t);
fft_coe=[1/0.54 1/( (27^2*2+23^2)/(50^2*2) )];
case 'blackman'
win = 0.42*ones(1,N) - 0.5*cos(2.0*pi*t) + 0.08*cos(4.0*pi*t);
fft_coe=[1/0.42 1/( (21^2*2+25^2+2^2)/(25^2*2^3) )];
end
end
コメント