MATLABを用いて信号処理の基礎の基礎
最終更新：2022/10/07
信号を読んで, 可視化して, FFTして, フィルターして, 逆FFTして, 可視化する.
研究室でフーリエ変換を教えるために作った教育用のコード. MATLABの公式のFFT関数の説明に少し書き足した感じ.
コード
使う信号
MATLABにデフォルトで入っている gong.mat という信号を使う.
信号 y と サンプリングレート Fs がロードされる.
ちなみに, sound(y, Fs) で音を聞ける.
ゴングの音かな.
MATLAB
load gong.mat
t = (0:length(y)-1)/Fs;
d = y; % 特に意味は無いが, dという名前に変えた
plot(t, d); xlabel('time [sec]'); ylabel('amp [?]')
ylim([-1.5, 1.5])
axis tight
信号の切り抜き
$t_1 = 1$ [sec] から $t_2 = 2$ [sec] を抜き出して可視化してみる.
MATLAB
%%%% t1とt2を入力 %%%%
t1 = 1; t2 = 2;
%%%%%%%%%%%%%%%%%%%%%
isIn = (t1 < t) & (t < t2);
dp = d(isIn); % data partially のつもり
tp = t(isIn); % time partially のつもり
plot(tp, dp); xlabel('time [sec]'); ylabel('amp [?]')
FFTをかける
本家の解説 に基づいて作成.
MATLAB
L = length(dp);
Y = fft(dp);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
plot(f, P1)
% plot(f, P1) % logをとると見やすい. 必要ならコメントアウト.
xlabel('f (Hz)')
ylabel('|P1(f)|')
axis tight
band-pass フィルター
$f_1 - f_2$ [Hz] を残してフィルタリング ( band-pass フィルター ).
MATLAB
%%%% f1とf2を入力 %%%%
f1 = 500; f2 = 1500;
%%%%%%%%%%%%%%%%%%%%%
isOmit = ~((f1 < f) & (f < f2));
Y_filtered = Y;
% 左側をフィルタ
tmp = Y_filtered(1:length(isOmit));
tmp(isOmit) = 0;
Y_filtered(1:length(isOmit)) = tmp;
% 右側をフィルタ
isOmit = fliplr(isOmit);
tmp = Y_filtered(end-length(isOmit)+1:end);
tmp(isOmit) = 0;
Y_filtered(end-length(isOmit)+1:end) = tmp;
% フィルタ後を可視化
plot(abs(Y_filtered))
xlabel('f (Hz)')
ylabel('|P1(f)|')
axis tight
逆FFT
逆FFTして ( ifft ) , フィルタの効果をみる.
MATLAB
d_filtered = real(ifft(Y_filtered));
t = tiledlayout(2, 1); t.Padding = 'compact'; t.TileSpacing = 'compact';
nexttile
plot(tp, dp);
xlabel('time [sec]'); ylabel('amp [?]'); legend('before filter', box='off')
nexttile
plot(tp, d_filtered);
xlabel('time [sec]'); ylabel('amp [?]'); legend('after filter', box='off')
よわからないけど, 拡大するとなめらかになっているかも.
MATLAB
t = tiledlayout(2, 1); t.Padding = 'compact'; t.TileSpacing = 'compact';
nexttile
plot(tp, dp);
ylabel('amp [?]'); legend('before filter', box='off')
xlim([1.1, 1.2])
nexttile
plot(tp, d_filtered);
xlabel('time [sec]'); ylabel('amp [?]'); legend('after filter', box='off')
xlim([1.1, 1.2])