% 信号処理分野のための GNU Octave/MATLAB/Scilab 入門
%  解答例 by Yoshiki NAGATANI 20180830 - https://ultrasonics.jp/nagatani/matlab/
% 
%  課題22：信号 u（7.1.節）をそのまま FFT した結果と，
%          信号の後ろにゼロ詰めしてから FFT した結果とを比較する。

clear;

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

% ゼロ詰め用のゼロ行列作成
zeropadding_ratio = 4;    % 4倍にゼロ詰め（すなわち 4-1=3 倍分のゼロを付け足す）
zeropadding = zeros(1,size(u_original,2)*(zeropadding_ratio-1));

% ゼロ詰め（ゼロ行列を後ろに結合）
u_zeropadding = [u_original zeropadding];

% FFT
u_original_fft = fft(u_original);
u_zeropadding_fft = fft(u_zeropadding);

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

% 波形の表示
subplot(3,1,1);
plot(u_zeropadding, 'ro-');
hold on;
plot(u_original, 'b*-');
grid on;
xlabel('サンプル [点]'); ylabel('振幅 [arb. unit]');
legend('ゼロ詰め後の波形', 'オリジナル波形');

% 振幅スペクトルの表示（リニアスケール）
subplot(3,1,2);
plot([0:size(u_zeropadding_fft,2)-1]/zeropadding_ratio, abs(u_zeropadding_fft), 'ro-');
hold on;
plot([0:size(u_original_fft,2)-1], abs(u_original_fft), 'b*-');
grid on;
xlabel('周波数 [点]'); ylabel('振幅スペクトル [arb. unit]');
legend('ゼロ詰め後の波形', 'オリジナル波形');

% 振幅スペクトルの表示（ログスケール(dB)）
subplot(3,1,3);
plot([0:size(u_zeropadding_fft_dB,2)-1]/zeropadding_ratio, u_zeropadding_fft_dB, 'ro-');
hold on;
plot([0:size(u_original_fft_dB,2)-1], u_original_fft_dB, 'b*-');
grid on;
xlabel('周波数 [点]'); ylabel('正規化した振幅スペクトル [dB]');
legend('ゼロ詰め後の波形', 'オリジナル波形');
