<span style="font-size: small;"><div class="statcounter"><a title="counters for myspace" href="http://www.statcounter.com/myspace/" target="_blank"><img class="statcounter" src="http://c.statcounter.com/5263102/0/20dc4843/1/" alt="counters for myspace"></a></div></span>

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