%***************************************%
% リアルタイム音声処理(Audio Toolbox)   %
% バイナリマスキング                    %
%         by Arata Kawamura (2021.01.13)%
%***************************************%
% 入力デバイス(マイク)の設定
clear;
deviceReader = audioDeviceReader;    % 読み込み(マイク)オーディオデバイス情報を取得
deviceReader.NumChannels = 2;        % ステレオ入力にする
setup(deviceReader);                 % デバイスを使えるようにする
% 出力デバイス(スピーカ)の設定
deviceWriter = audioDeviceWriter(...
    'SampleRate', deviceReader.SampleRate);  % 書き出し(スピーカ)オーディオデバイス情報を取得
setup(deviceWriter,zeros(deviceReader.SamplesPerFrame,deviceReader.NumChannels)); %% デバイスを使えるようにする
%************************************%
%            基本項目設定            %
% Fs :サンプリング周波数(44100)      %
% out:出力バッファ (長さ1024)        %
%************************************%
Fs  = deviceReader.SampleRate;       % サンプリング周波数
Buf = deviceReader.SamplesPerFrame;  % デバイスから1回に読み込む信号長
out = zeros(Buf, 2);                 % 出力デバイスに送る信号
%************************************%
%            録音用設定              %
%************************************%
RecTime = 10;                        % 実行時間の設定 [秒]
l   = 0;                             % 録音用のフレームインデクス
in_wav  = zeros(RecTime*Fs,2);       % 録音用入力信号
out_wav = zeros(RecTime*Fs,2);       % 録音用出力信号
%************************************%
%       信号処理用変数設定           %
%************************************%
N = Buf;                             % 入出力信号の配列の長さ--> FFTの長さ(1024)
xin1 = zeros(N, 2);                  % 入力信号格納
x = zeros(N*2, 2);                   % FFT用の時間信号
y = zeros(N*2, 2);                   % IFFT用の時間信号
w = window('hann', N*2);             % 窓関数
%************************************%

%*********************************************************%
% メインループ                                            %
% マイク入力がリアルタイムでスピーカ出力される            %
% RecTime[秒]経過後，入力と出力がwavファイルに記録される  %
%*********************************************************%
tic;
while toc <= RecTime
    audioData = deviceReader();        % バッファの長さ(1024個)のマイク入力信号
    %plot(1:length(audioData),audioData);ylim([-0.2,0.2]);drawnow % 波形の表示
    xin = audioData;                   % 入力信号を格納
    x   = [xin1; xin] .* w;            % 窓関数を掛ける
    X   = fft(x);                      % FFT(長さ2N)
    %***************************%
    %    信号処理部             %
    %    X: 入力信号スペクトル  %
    %    Y: 出力信号スペクトル  %
    %***************************%

    G   = abs( X(:,1) ) - abs( X(:,2) );  % LRの差分をとる
    G(G>0)   = 1;                         % 負ならゲイン1 (左抽出)
    G(G<0)   = 0;                         % 正ならゲイン0 (右除去)
    Y   = G .* X;                         % バイナリマスキング

    %************************%
    %   信号処理部ここまで   %
    %************************%
    y = y + ifft(Y); % ハーフ・オーバラップでyの後半1フレームをつくる 
    out          = y(1:N, :);          % yの前半の半フレーム(Nサンプル)を出力
    y(1:N,:)      = y(N+1:end,:);      % yの後半の1フレーム(2N)を前半にシフト
    y(N+1:end,:) = 0;                  % yの後半の半フレームを0にする
    xin1   = xin;                      % 入力信号を次回FFTの前半として記録
    %******************************************************************%
     in_wav(N*l+1:N*l+N, :) = xin;     % 入力信号の記録(録音用)
    out_wav(N*l+1:N*l+N, :) = out;     % 出力信号の記録(録音用)
    l = l+1;                           % 録音用時刻の更新
    %******************************************************************%
    deviceWriter(out);                 % スピーカデバイスから out を出力
end

%% デバイスの解放
release(deviceReader);
release(deviceWriter);

%% 入力と出力をwavファイルに保存
audiowrite('input.wav', in_wav, Fs);
audiowrite('output.wav', out_wav, Fs);
