% Sample MATLAB Code for Audio Analysis
% 2019.6.26 by T. Hirano
%
% 【準備】
% 1. 左上の「現在のフォルダー」のフォルダから緑色矢印が出ているアイコンを使って
%    このファイルがあるフォルダに設定する。
% 2. 「エディター」タブの「実行」ボタンを押すと実行される。
%    (または、下の「コマンド ウィンドウ」で.mを除いたこのファイル名
%    "audio_analyze"と入力してEnterを押す。)
% 3. 下記コードの行頭の"%%"はセクションという意味で、
%    「エディター」タブの「実行して次に進む」でステップ実行できる。
%

%% Read Audio File
%
% http://www.cs.tut.fi/sgn/arg/intro/basics/pia60.wav
%
[y,Fs] = audioread('pia60.wav');    % y: 波形, Fs: サンプリング周波数 (Hz)

dt=1/Fs;    % サンプリング間隔 (sec)
duration=(length(y)-1)*(1/Fs);  % waveファイルの再生時間(T)
time_lis=[0:1/Fs:duration];

plot(time_lis,y);

%% Play Audio File
sound(y,Fs);

% モノラルの場合はyは1次元配列
% ステレオの場合はyは2次元配列
% 2次元配列→1次元配列にするには
% y1=y(:,1)
% とすればよい。

%% FFT
% abs: 絶対値, anble: 偏角
% 「通信信号処理」2018年度テキスト、第8回, p.8参照
ft_y=fft(y);

df=1/duration;
freq_lis=[0:df:df*(length(ft_y)-1)];

figure; % 新しくグラフウィンドウを開く（これがないと前のグラフが更新される）
plot(freq_lis,abs(ft_y));

%% Cut Waveform
f0=263.7;       % 基本周波数 (Hz) （上のグラフから読み取って入力する）
period=1/f0;    % 周期 (sec)
n0=floor(length(y)/2);  % 何サンプル目から取り出すか？
n=round(3*period/dt);   % 何サンプル分を取り出すか？
figure;
plot(y([n0:n0+n]));

%% Copy One Period
f0=263.7;       % 基本周波数 (Hz)
period=1/f0;    % 周期 (sec)
n0=floor(length(y)/2);  % 何サンプル目から取り出すか？
n=round(period/dt);     % 何サンプル分を取り出すか？

y_one_period=y([n0:n0+n]);      % 1周期を切り出す
n_period=10;                    % 何周期分コピーするか？
y2=repmat(y_one_period,[n_period,1]);    % 1周期の配列をコピーして繰り返す
figure;
plot(y2);

%% Play Sound
n_period=1000;      % 何周期分コピーするか？
y3=repmat(y_one_period,[n_period,1]);    % 1周期の配列をコピーして繰り返す
sound(y3,Fs);

%% Write Wavefile
audiowrite('repat1000.wav',y3,Fs);    % 引数は(ファイル名, 波形, サンプリングレート)