%%% VSS-MDNLMS simulation code
%%%
%%% Marius Rotaru, Felix Albu, Henri Coanda, “A Variable Step Size Modified Decorrelated NLMS Algorithm for Adaptive Feedback Cancellation in Hearing Aids",
%%% in Proc. of ISETC 2012, Timisoara, Romania, 15-16 November 2012, pp. 263-266 ISBN: 978-146731176-2 DOI: 10.1109/ISETC.2012.6408070
%%%
%%% Felix Albu
%%% Valahia University of Targoviste, Romania
%%% Email: felix.albu@valahia.ro
%%% http://falbu.50webs.com
%% clear buffers
clear all; clc;
close all;
%% Input Signal
SNR = 20;
%% Parameters for AFC
epsi = 1.e-3; %%% small constant for preventing divide zero
pu = 2; %%% estimated power using first order recursion
Mis_Avg=[];
%%% Compensation Gain
gain = [20, 30, 40, 50, 60];
a_gain = 0.5 * gain;
cgain = [10^(a_gain(1)/20),10^(a_gain(2)/20),10^(a_gain(3)/20), ...
10^(a_gain(4)/20),10^(a_gain(5)/20)];
max_gain = cgain;
if 1
%
%[xn, Fs] = wavread('test_file_TIMIT.wav');
[xn, Fs] = audioread('test_file_TIMIT.wav');
xn = 0.5 * xn; %%% weighting
vd = var(xn);
vn = vd / 10^(SNR/10);
vv = sqrt(vn) * randn(length(xn), 1);
xn = xn + vv;
else
Fs = 16000;
len = Fs*32;
coeff = [ 1 -0.9 0.8 -0.7 0.6 -0.5 0.4 -0.3 0.2 -0.1 0.05 ];
vn = randn(len, 1);
xn = filter(1, coeff, vn); %%% input AR noise
vd = var(xn);
vn = vd / 10^(SNR/10);
vv = sqrt(vn) * randn(len, 1); %%% background noise
xn = xn + vv; %%% Final input signal
end
xn = fix(xn * 2^15);
%% Input parameters
P = 10; %%% prediction order
M = 128; %%% FFT size
N = M / 2; %%% overlap size
L = 64; %%% filter order
N0 = 64;
f0=3;
e = zeros(N0,1);
mu = 0.004;
tmu = mu;
len = length(xn);
frame = fix(len / N);
%% Feedback Path
fid = fopen('coef2.f','r');
Hf = fread(fid, inf, 'float');
Hf = Hf(1:L);
fclose(fid);
normH = sqrt( Hf' * Hf );
flag_p = 0;
ch_cnt = 1;
ch_path = [ 0 ];
ch_path = round((ch_path * Fs) / N);
%% Clear Filter Bank buffers for DHA System
xp = zeros(N, 1);
ep = zeros(N, 1);
iep = zeros(N, 1);
uep = zeros(N, 1);
%% Sine window
swin = sin(pi*(0.5/(2*N):1/(2*N):1))';
%% Feedback Cancellation Buffers
err = zeros(N, 1);
eroare = zeros(N, 1);
des = zeros(N, 1);
Wt = zeros(L, 1);
Wf = zeros(L, 1);
x = zeros(1, L);
u = zeros(L, 1);
xx = zeros(L, 1);
xvec = zeros(L, 5) + epsi;
evec = zeros(L, 5) + epsi;
uvec = zeros(L, 5) + epsi;
%% Linear Prediction Parmeters & Vector
% Burg Algorithms
lam = 1-(1/160); %%% forgetting factor for linear predictor
dm = zeros(P, 1) + epsi;
nm = zeros(P, 1) + epsi;
km = zeros(P, 1) + epsi;
f = zeros(P, 1) + epsi; %%% adaptive forward linear predictor vector
b = zeros(P, 2) + epsi; %%% adaptive backward linear predictor vector
fx = zeros(P, 1) + epsi; %%% forward linear predictor vector for output signal
bx = zeros(P, 2) + epsi; %%% backward linear predictor vector for output signal
fe_n = zeros(P, 1) + epsi; %%% forward linear predictor vector for error signal for MD-NLMS
be_n = zeros(P, 2) + epsi; %%% backward linear predictor vector for error signal for MD-NLMS
%% Output Vectors
MisT = []; Mis = [];
yn = []; yx = [];
%% DHA system buffers and parameters
maxv = 2^16;
env = zeros(5, 1);
freq = [500, 1000, 2000, 4000, 8000];
s_f = M / (2 * Fs);
%%% Number of between center frequency each band
CenF = [ (freq(1) * s_f), ((freq(2)+freq(1)) * s_f), ...
((freq(3)+freq(2)) * s_f), ((freq(4)+freq(3)) * s_f), ...
((freq(5)+freq(4)) * s_f), (2 * freq(5) * s_f) ];
CenF(2:end) = CenF(2:end) - CenF(1:end-1) + 1;
%% Compensation Gain
A1 = linspace(1,cgain(1),CenF(1));
A2 = linspace(cgain(1),cgain(2),CenF(2));
A3 = linspace(cgain(2),cgain(3),CenF(3));
A4 = linspace(cgain(3),cgain(4),CenF(4));
A5 = linspace(cgain(4),cgain(5),CenF(5));
A6 = linspace(cgain(5),1,CenF(6)+1);
Ao = [ A1, A2(2:end), A3(2:end), A4(2:end), A5(2:end), A6(2:end) ];
Wgain = Ao';
IGain = Wgain;
%% Variable Step Size Params
vss = 1; %%% VSS Enable/Disable
lam1 = 0.992; %%% forgetting factor for vss power estimation
lam2 = 0.992; %%% forgetting factor for mu averaging
Pe = 0; %%% Power Estimation of Error Signal
Pd = 0; %%% Power Estimation of Desired Signal
mu_tmp1 = zeros(len,1);
mu_tmp4 = zeros(len,1);
mu_tmp2 = 0;
mu_tmp3 = 0;
mu_max = 0.005;
mu_min = 0.0005;
%% Main Algorithm Loop
for ii = 1: frame-1
%% Filter Bank Input
ec = err;
et = [ep; ec]; %% Append the new frame to the previous to make M frame
et = et .* swin; %% Windowing
ep = ec;
%% FFT
Et = fft(et) / M;
%% Compute Compensation Gain, Reference Signal
Xo = Et(1:N+1) .* Wgain;
%% Synthesis for Ouptut DHA
xo = ifft([Xo; conj(Xo(end-1:-1:2))]) * M;
xo = xo .* swin;
xc = xp + xo(1:N);
xp = xo(N+1:2*N);
%% Feedback path change
if ii == floor((frame-1)/2)
flag_p = 1;
ch_cnt = ch_cnt + 1;
end
if flag_p == 1
Hf = [ 0.02*randn(5,1);Hf(1:L-5)];
normH = sqrt( Hf' * Hf );
flag_p = 0;
end
%% Synthesis for inverse gain
iEt = Et(1:N+1) ./ IGain;
ieo = ifft([iEt; conj(iEt(end-1:-1:2))]) * M;
ieo = ieo .* swin;
iec = iep + ieo(1:N);
iep = ieo(N+1:2*N);
%% Synthesis for unit out
ueo = ifft(Et) * M;
ueo = ueo .* swin;
uec = uep + ueo(1:N);
uep = ueo(N+1:2*N);
%% stock buffers
xvec = [ xvec(:,2:end), xc ];
evec = [ evec(:,2:end), err ];
uvec = [ uvec(:,2:end), uec ];
%% Adaptive Feedback Cancellation
for n = N0: N
x = [ xc(n), x(1: end-1) ];
%% desired & error signal
des(n) = xn((ii-1)*N+n) + x * Hf;
err(n) = des(n) - x * Wf;
if 1
hv(n)=err(n);
err(n) = 2*tanh(0.5*err(n)); % limit the error signal
eroare(n)=abs(hv(n)-err(n));
aaa = (abs(hv(n)-err(n))>0.15);
%%%%%frequency shift
eh=real(hilbert(err(n:-1:n-N0+1)').*exp(-1i*2*pi*(f0/Fs).*[n:-1:n-N0+1]));
err(n) = eh(1);
end
%% saturate error signal
if abs(err(n)) > maxv
flag_lp = 0;
else
flag_lp = 1;
end
%% linear predictor
fxx = evec(n,2);
fer = evec(n,4);
ier = iec(n);
f(1) = fer; b(1, 1) = fer;
fx(1) = fxx; bx(1, 1) = fxx;
fe_n(1) = ier; be_n(1, 1) = ier;
if flag_lp
for m = 2: P
%%% Burg Lattice Algorithm
dm(m-1) = lam * dm(m-1) + (1-lam) * (f(m-1)^2 + b(m-1,2)^2);
nm(m-1) = lam * nm(m-1) + (1-lam) * (-2)*(f(m-1)*b(m-1,2));
km(m) = nm(m-1) / (dm(m-1));
%%% Adaptive Lattice Predictor
f(m) = f(m-1) + km(m)*b(m-1,2);
b(m,1) = b(m-1,2) + km(m)*f(m-1);
%%% Weighted Lattice Predictor for Updating Input Vector
fx(m) = fx(m-1) + km(m)*bx(m-1,2);
bx(m,1) = bx(m-1,2) + km(m)*fx(m-1);
%%% Lattice Predictor for Updating Error Vector of MD-NLMS
fe_n(m) = fe_n(m-1) + km(m) * be_n(m-1,2);
be_n(m,1) = be_n(m-1,2) + km(m) * fe_n(m-1);
end
b(:,2) = b(:,1);
bx(:,2) = bx(:,1);
be_n(:,2) = be_n(:,1);
else
for m = 2: P
%%% Weighted Lattice Predictor for Updating Input Vector
%%% PAP/MNLMS
fx(m) = fx(m-1) + km(m)*bx(m-1,2);
bx(m,1) = bx(m-1,2) + km(m)*fx(m-1);
%%% Lattice Predictor for Updating Error Vector of MNLMS
fe_n(m) = fe_n(m-1) + km(m) * be_n(m-1,2);
be_n(m,1) = be_n(m-1,2) + km(m) * fe_n(m-1);
end
bx(:,2) = bx(:,1);
be_n(:,2) = be_n(:,1);
end
%% pseudo input vector & pseudo error vector
u = [ fx(end); u(1:end-1) ];
es = fe_n(end);
%% compute power
pu = u' * u;
%% Update Gradient
dW = es / ( pu ) * u;
if (vss==1)
Pe = lam1 * Pe + (1 - lam1)*(err(n)^2);
Pd = lam1 * Pd + (1 - lam1)*(des(n)^2);
mu_tmp1((ii-1)*N+n) = mu_max*abs(1 - (Pd/Pe));
mu_tmp2 = lam2*mu_tmp2 + (1-lam2)*mu_tmp1((ii-1)*N+n);
%% limit mu between mu_max/mu_min
mu_tmp3 = min (mu_max,mu_tmp2);
mu = max (mu_min, mu_tmp3);
%% just for plotting purpose
mu_tmp4((ii-1)*N+n) = mu_tmp3;
end
%% Update Weight Vector
if ii > 5
Wt = Wt + mu * dW;
end
%% Compute Misalignment
wet = Hf - Wt;
Mis = [ Mis; [sqrt(wet'*wet)/normH]];
end
Wf = Wt;
%% Append Input Signal of Hearing Aids
yx = [ yx; [err,xc] ];
%% Apend
if mod(ii, 100) == 0
MisT = [ MisT; Mis ]; Mis = [];
yn = [ yn; yx ]; yx = [];
end
end
%% Apend Rest
MisT = [ MisT; Mis ];
yn = [ yn; yx ];
Mis_Avg = [Mis_Avg 20*log10(MisT)];
%% Plots
t1 = (1:length(MisT))/Fs;
t2 = (1:length(yn))/Fs;
%%% error and output signal
figure,
subplot(4,1,1), plot(t2, yn(:,1));
axis([0,t2(end),-2*max(abs(xn)),2*max(abs(xn))]);
title('Error Signal');
xlabel('Seconds');
subplot(4,1,2), plot(t2, yn(:,2));
axis([0,t2(end),-min(a_gain)*max(abs(xn)),min(a_gain)*max(abs(xn))]);
title('Output Signal');
xlabel('Seconds');
subplot(4,1,3), plot(t2, xn(1:length(yn)));
axis([0,t2(end), -2*max(abs(xn)),2*max(abs(xn))]);
title('Input Signal');
xlabel('Time in sec.');
subplot(4,1,4), plot(t2, mu_tmp4(1:length(yn)));
axis([0,t2(end), 0,0.01]);
title('Mu');
xlabel('Seconds');
%%% Misalignment
figure, plot(t1, 20*log10(MisT));
xlim([0,t1(end)]);
%axis([0 32 -35 5])
Mis_Avg = 20*log10(MisT);