% 信号処理分野のための GNU Octave/MATLAB/Scilab 入門
%  解答例 by Yoshiki NAGATANI 20180830 - https://ultrasonics.jp/nagatani/matlab/
% 
%  課題23：スイープ音（連続的に正弦波の周波数が変化する）を作成して
%          聞いてみる。また，そのスペクトルも確認する。
%  課題24：課題23で作成した信号のスペクトログラムを表示する。

clear;

% 各パラメータ設定
fs = 8000;           % サンプリングレート [Hz]
T  = 0.5;            % 時間長 [s]

% 時間軸
dt = 1.0/fs;         % 時間刻み [s]
t  = [0:T*fs-1]/fs;  % 時間軸

% 周波数の生成
freq_start =   10.0;    % 開始周波数 [Hz]
freq_end   = 3500.0; % 終了周波数 [Hz]
sig_freq = (t/T) * (freq_end - freq_start) + freq_start;
subplot(4,1,1); plot(t, sig_freq, 'b.');
grid on; title('周波数');
xlabel('時刻 [s]'); ylabel('周波数 [Hz]');

% 周波数 → 位相の作成（ 2πf = dφ/dt → φ=∫2πfdt  ）
sig_phase = cumsum(sig_freq) * dt * (2.0 * pi);
subplot(4,1,2); plot(t, sig_phase, 'r.');
grid on; title('位相（時刻 0 [s] 時点の位相を 0 [rad] とする）');
xlabel('時刻 [s]'); ylabel('位相 [rad]');

% 位相 → 時間波形の作成（ 信号 = sin(位相[rad]) ）
sig_waveform = sin(sig_phase);
subplot(4,1,3); plot(t, sig_waveform, 'k-');
grid on; title('時間波形'); xlabel('時刻 [s]'); ylabel('振幅 [arb. unit]');

% スペクトログラム（課題24）
fft_win = 128;   % FFT ウインドウサイズ [点]（この解答例ではオーバーラップはナシ）
sig_specgram = zeros(floor(size(sig_waveform,2)/fft_win), fft_win); % 空の配列を用意
for i=1:floor(size(sig_waveform,2)/fft_win) % 時刻方向に移動
    window_start = 1+(i-1)*fft_win;            % 解析区間の先頭の座標
    window_end   = window_start + fft_win - 1; % 解析区間の末尾の座標
    sig_specgram(i,:) = abs(fft(sig_waveform(window_start:window_end))); % FFT
end
subplot(4,1,4); imagesc(sig_specgram(:,:)');
grid on; title('スペクトログラム（折り返し部分も表示）');
xlabel('FFT窓の通し番号（時刻に比例）'); ylabel('FFT結果の点数（周波数に比例）');
drawnow;
% ↑簡単のために image 関数に渡せる形の行列として保存しているため，
% グラフ表示はこのままでは適切ではない。（図の下に行くほど周波数が高い。
% また，値も [s] や [Hz] にはなっていない。この点の改善は各自の課題とする。）

% スペクトル（課題23）
sig_fft = fft(sig_waveform); % FFT
figure;
plot([0:size(sig_fft,2)-1]/T, abs(sig_fft), 'b.-');
grid on; title('スペクトル（波形全体をフーリエ変換した結果）');
xlabel('周波数 [Hz]'); ylabel('振幅スペクトル [arb. unit]');
drawnow;

% 音の再生
soundsc(sig_waveform, fs);
