[MATLAB]RIR_reshapingを用いた残響抑制のプログラム

clear;

%インパルス応答読み込み
filename = '%%%%%%%.wav';
[c,Fs] = audioread(filename);

%パラメータ設定
iter = 100000; %更新回数
mwu = 10^(-7); %更新の際のステップサイズ
Lh = round(length(c)*1+1); %フィルタ長
Lc = length(c);%インパルス応答の長さ
Lg = Lc + Lh - 1;%新しく作成するインパルス応答の長さ

%フィルタの初期設定
h = zeros(Lh,1); 
h(1) = 0.01;%1つ目の要素が0.01で残りの要素が全て0となるように設定

%畳み込み行列Cの設定
C = zeros(Lg+Lh-1,1);
j=1;
for i = Lg:-1:Lh
   C(i) = c(j); 
   j=j+1;
end

%set N1N2N3
[M,N1] = max(abs(c));%直接音が到達するサンプル数をN1に代入
N2 = round(0.004*Fs);%直接音部のサンプル数を算出しN2に代入
N3 = Lg-N1-N2;%残響部の長さをN3に代入
N0 = N1 + round(0.2*Fs);%マスキング曲線の計算に必要な値をN0に代入

%窓関数の設定
%直接音部の窓関数Wd
Wd = [zeros(N1,1);ones(N2,1);zeros(N3,1)];%直接音からN2サンプルだけ1に設定
    
%残響部の窓関数Wu
n = N1+N2+1:N3;
W0(n) = 10.^(3/log(N0/(N1+N2))*log(n/(N1+N2))+0.5);%マスキング限界
Wu = [zeros(N1+N2,1);W0.'];%残響部に関数W0がかかるように設定

%学習部分


for p = 1:iter
   cpad = [c;zeros(Lh-1,1)];%cの後ろに0をパディング
   hpad = [h;zeros(Lc-1,1)];%hの後ろに0をパディング
   g = ifft(fft(cpad).*fft(hpad));%フーリエ変換、逆変換を利用して畳み込み和を計算

   [Md,Id] = max(abs(g.*Wd));%直接音部の絶対値が最も大きいインデックスをI_dに代入
   [Mu,Iu] = max(abs(g.*Wu));%残響部の絶対値が最も大きいインデックスをI_uに代入

   C_I_d = C(Lg+1-Id:Lg+Lh-Id);%畳み込み行列CからId行を抜き出す
   C_I_u = C(Lg+1-Iu:Lg+Lh-Iu);%畳み込み行列CからIu行を抜き出す
   h = h - mwu*(C_Iu./g(Iu) - C_I_d./g(I_d));
   %hにgdの最大値が大きく、guの最小値が小さくなるような勾配を与えて更新
end
 
%結果をプロット
n = (1:Lc);
y_c = sp2dB(c(n));%元のインパルス応答の振幅特性
y_g = sp2dB(g(n));%再形成されたインパルス応答の振幅特性
log10W0(n) = -20*log10(W0(n));%マスキング限界を振幅特性の式に対応させる

plot(n,y_c,n,y_g,n,log10W0)
xlabel('Samples')
ylabel('Magnitude(dB)')

function y = sp2dB(x)
M = max(abs(x)); %入力された信号の最も大きい絶対値を取得
y = 20.*log10(X/M); %Mを基準として入力された信号をdBで表示