% 信号処理分野のための GNU Octave/MATLAB/Scilab 入門
%  解答例 by Yoshiki NAGATANI 20180830 - https://ultrasonics.jp/nagatani/matlab/
% 
%  課題21：信号 u（7.1.節）にハン窓（ハニング窓）などを掛けてFFTし，
%          矩形窓（窓なし）と結果を比較する。

clear;

% 元の信号
u_original = sin([0:0.5:10]);

% Hann窓【 0.5*(1-cos(2πn/N) 】
hann_win = 0.5 * ( 1.0 - cos( 2.0 * pi * [0:size(u_original,2)-1]/size(u_original,2) ) );

% 窓掛け
u_with_win = u_original .* hann_win;

% FFT
u_original_fft = fft(u_original);
u_with_win_fft = fft(u_with_win);

% dBへの変換（オリジナル波形の振幅スペクトルの最大値で正規化）
u_original_fft_dB = 20.0 * log10( abs(u_original_fft) / max(abs(u_original_fft)) );
u_with_win_fft_dB = 20.0 * log10( abs(u_with_win_fft) / max(abs(u_original_fft)) );

% 波形の表示
subplot(3,1,1);
plot(u_original, 'bo-');
hold on;
plot(u_with_win, 'ro-');
plot(hann_win, 'k--');
grid on;
xlabel('サンプル [点]'); ylabel('振幅 [arb. unit]');
legend('オリジナル波形', 'Hann窓を掛けた波形', 'Window');

% 振幅スペクトルの表示（リニアスケール）
subplot(3,1,2);
plot([0:size(u_original_fft,2)-1], abs(u_original_fft), 'bo-');
hold on;
plot([0:size(u_with_win_fft,2)-1], abs(u_with_win_fft), 'ro-');
grid on;
xlabel('周波数 [点]'); ylabel('振幅スペクトル [arb. unit]');
legend('オリジナル波形', 'Hann窓を掛けた波形');

% 振幅スペクトルの表示（ログスケール）
subplot(3,1,3);
plot([0:size(u_original_fft_dB,2)-1], u_original_fft_dB, 'bo-');
hold on;
plot([0:size(u_with_win_fft_dB,2)-1], u_with_win_fft_dB, 'ro-');
grid on;
xlabel('周波数 [点]'); ylabel('正規化した振幅スペクトル [dB]');
legend('オリジナル波形', 'Hann窓を掛けた波形');
