&amp;amp;amp;amp;amp;lt;div class="statcounter"&amp;amp;amp;amp;amp;gt;&amp;amp;amp;amp;amp;lt;a title="counters for myspace" href="http://www.statcounter.com/myspace/" target="_blank"&amp;amp;amp;amp;amp;gt;&amp;amp;amp;amp;amp;lt;img class="statcounter" src="http://c.statcounter.com/5263102/0/20dc4843/1/" alt="counters for myspace"&amp;amp;amp;amp;amp;gt;&amp;amp;amp;amp;amp;lt;/a&amp;amp;amp;amp;amp;gt;&amp;amp;amp;amp;amp;lt;/div&amp;amp;amp;amp;amp;gt;

%%% 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

%% 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);