% 信号処理分野のための GNU Octave/MATLAB/Scilab 入門
%  解答例 by Yoshiki NAGATANI 20180830 - https://ultrasonics.jp/nagatani/matlab/
% 
%  課題26：インパルス，もしくは4.5.で作成したホワイトノイズに
%          LPF（低域通過フィルター），HPF（高域通過フィルター），
%          BPF（帯域通過フィルター）などを適用して，波形の変化を観察すると
%          ともに，聞き比べる。また，その振幅スペクトルを表示して確認する。

clear;

% 各パラメータ設定
fs = 8000;                 % サンプリングレート [Hz]
T  = 0.1;                  % 時間長 [s]
cutoff_freq_LPF_Hz = 630;  % LPFのカットオフ周波数 [Hz]
cutoff_freq_HPF_Hz = 500;  % HPFのカットオフ周波数 [Hz]

% 時間軸・周波数軸
dt = 1.0/fs;         % 時間刻み [s]
t  = [0:T*fs-1]/fs;  % 時間軸
df = 1/T;            % 周波数刻み [Hz]
f  = [0:T*fs-1]*df;  % 周波数軸


% オリジナル波形 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% インパルス
sig_original = zeros(1,floor(T*fs));
sig_original(floor(size(sig_original,2)/2)) = 1.0; 
% ホワイトノイズ（次の行のコメントアウトを削除）
% sig_original = rand(1,floor(T*fs)) - 0.5;

% FFT（この結果に対して各種フィルターを適用する）
sig_original_fft = fft(sig_original);

% 波形の表示
subplot(4,2,1);
plot(t, sig_original, 'k.-');
grid on; title('オリジナル波形');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]');
% 振幅スペクトルの表示
subplot(4,2,2);
plot(f, abs(sig_original_fft), 'k.-');
ylim([0 max(abs(sig_original_fft))*1.2]);
grid on; title('オリジナル波形の振幅スペクトル');
xlabel('周波数 [Hz]'); ylabel('振幅 [arb. unit]');

% LPF（低域通過フィルター） %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cutoff_freq_LPF_points_low  = round(cutoff_freq_LPF_Hz/df) + 1;
cutoff_freq_LPF_points_high = size(sig_original_fft,2) - round(cutoff_freq_LPF_Hz/df) + 1;

% LPF 適用
sig_LPF_fft = sig_original_fft;
sig_LPF_fft(cutoff_freq_LPF_points_low:cutoff_freq_LPF_points_high) = 0.0;

% iFFT
sig_LPF = real(ifft(sig_LPF_fft));

% 波形の表示
subplot(4,2,3);
plot(t, sig_LPF, 'r.-');
grid on; title('LPF波形');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]');
% 振幅スペクトルの表示
subplot(4,2,4);
plot(f, abs(fft(sig_LPF)), 'r.-');
grid on; title('LPF波形の振幅スペクトル');
xlabel('周波数 [Hz]'); ylabel('振幅 [arb. unit]');

% HPF（高域通過フィルター） %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cutoff_freq_HPF_points_low  = round(cutoff_freq_HPF_Hz/df) + 1;
cutoff_freq_HPF_points_high = size(sig_original_fft,2) - round(cutoff_freq_HPF_Hz/df) + 1;

% HPF 適用
sig_HPF_fft = sig_original_fft;
sig_HPF_fft(1:cutoff_freq_HPF_points_low) = 0.0;
sig_HPF_fft(cutoff_freq_HPF_points_high:end) = 0.0;

% iFFT
sig_HPF = real(ifft(sig_HPF_fft));

% 波形の表示
subplot(4,2,5);
plot(t, sig_HPF, 'g.-');
grid on; title('HPF波形');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]');
% 振幅スペクトルの表示
subplot(4,2,6);
plot(f, abs(fft(sig_HPF)), 'g.-');
grid on; title('HPF波形の振幅スペクトル');
xlabel('周波数 [Hz]'); ylabel('振幅 [arb. unit]');

% BPF（帯域通過フィルター） %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% BPF 適用（前述の LPF & HPF を両方適用）
sig_BPF_fft = sig_original_fft;
sig_BPF_fft(cutoff_freq_LPF_points_low:cutoff_freq_LPF_points_high) = 0.0;
sig_BPF_fft(1:cutoff_freq_HPF_points_low) = 0.0;
sig_BPF_fft(cutoff_freq_HPF_points_high:end) = 0.0;

% iFFT
sig_BPF = real(ifft(sig_BPF_fft));

% 波形の表示
subplot(4,2,7);
plot(t, sig_BPF, 'b.-');
grid on; title('BPF波形');
xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]');
% 振幅スペクトルの表示
subplot(4,2,8);
plot(f, abs(fft(sig_BPF)), 'b.-');
grid on; title('BPF波形の振幅スペクトル');
xlabel('周波数 [Hz]'); ylabel('振幅 [arb. unit]');

drawnow;

% 音の再生 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
soundsc(sig_original,fs); pause(T+0.5);
soundsc(sig_LPF,fs);      pause(T+0.5);
soundsc(sig_HPF,fs);      pause(T+0.5);
soundsc(sig_BPF,fs);
