Can anyone help me in understanding an algorithm and tell me where am i wrong in simulating the adaptive filter

Hi

I came across a paper titled 'Framework for automated earthquake event detection based on denoising By Adaptive Filter'.The thing which fancied me is that based on how they achieved convergence by adding/subtracting a small value called 'mu_extra'.Based on the paper,i also tried to simulate the same but i didnt get the output as expected.

i have also attached the paper along with it

kindly check and help me

clc;
clear all;
close all; 
%sys_desired = [86 -294 -287 -262 -120 140 438 641 613 276 -325 -1009 -1487 -1451 -680 856 2954 5206 7106 8192 8192 7106 5206 2954 856 -680 -1451 -1487 -1009 -325 276 613 641 438 140 -120 -262 -287 -294 86] * 2^(-15);
%sys_desired = [-3 0 19 32 19 0 -3]*2^(-15);
%a=importdata("G:\Project\BKHL.HHE.new.dat");
digfilt = designfilt('lowpassfir', 'PassbandFrequency', 10, 'StopbandFrequency',20,'SampleRate', 100);
for itr=1:100
    
    %% Defining input and initial Model Coefficients
    
    %input 
    x=randn(1,60000);
    %a=importdata("G:\Project\BKHL.HHE.new.dat");
    %x=a';
    
    %% EVSS LMS 
    model_coeff_evss = zeros(1,length(digfilt.Coefficients));
    
    %% Initial Values of Model Tap
    model_tap = zeros(1,length(digfilt.Coefficients));
    
    %% System Output where a 40 dB Noise floor is added
    noise_floor = 40;
    sys_opt = filter(digfilt.Coefficients,1,x)+awgn(x,noise_floor)-x;
    %sys_opt = awgn(x,noise_floor); 
    
    mu_min = 0.007;
    
    % input variance
    input_var = var(x);
    
    % upper bound = 1/(filter_length * input variance)
    mu_max = 1/(input_var*length(digfilt.Coefficients));
    %mu_max=0.025;
    
    %% Defining initial parameters for EVSS-LMS algorithm
    mu_EVSS(1)=0.025;
    mu_EXTRA=1.5*10^(-4);
    %% EVSS LMS ALGORITHM
    for i=1:length(x)
        
        % model tap values (shifting of tap values by one sample to right)
        model_tap1=[x(i) model_tap(1:end-1)];
        
        % model output
        model_out_evss(i) = model_tap1 * model_coeff_evss';
        
        % error
        e_evss(i) = sys_opt(i) - model_out_evss(i);
        
        %Updating the coefficients
        model_coeff_evss = model_coeff_evss + mu_EVSS(i) * e_evss(i) * model_tap;
        if mu_EVSS(i)>mu_min
            c=1;
        else mu_EVSS(i)<=mu_min;
            c=2^( 5);
        end
        mu_EVSS(i+1) = c * mu_EVSS(i) + mu_EXTRA * sign(e_evss(i));
    end
    %% Storing the e_square values after a whole run of VSS LMS algorithm
    
    err_EVSS(itr,:) = e_evss.^2;
    
    %% Printing the iteration number
    clc
    disp(char(strcat('iteration no : ',{' '}, num2str(itr) )))
end
figure;
plot(10*log10(mean(err_EVSS)), '-b');
title('Error Curve Plot for VSS LMS Algorithm'); xlabel('iterations');ylabel('MSE(dB)');legend('VSS LMS Algorithm')
grid on

for the above code I conceived from these paper,the mse is not converging.

why is it so.

if my code is wrong,kindly correct me.

Answer 1

[Mathieu NOE]

hello again

try this updated code

some minor bugs and missing mu min / max bounding code (see comments with % <= here)

clc;
clear all;
close all; 
%sys_desired = [86 -294 -287 -262 -120 140 438 641 613 276 -325 -1009 -1487 -1451 -680 856 2954 5206 7106 8192 8192 7106 5206 2954 856 -680 -1451 -1487 -1009 -325 276 613 641 438 140 -120 -262 -287 -294 86] * 2^(-15);
%sys_desired = [-3 0 19 32 19 0 -3]*2^(-15);
%a=importdata("G:\Project\BKHL.HHE.new.dat");
digfilt = designfilt('lowpassfir', 'PassbandFrequency', 10, 'StopbandFrequency',20,'SampleRate', 100);
sys_desired = digfilt.Coefficients;
length_fir = length(sys_desired);

for itr=1:100
    
    %% Defining input and initial Model Coefficients
    
    %input 
    x=randn(1,3000);
    %a=importdata("G:\Project\BKHL.HHE.new.dat");
    %x=a';
    
    %% EVSS LMS 
    model_coeff_evss = zeros(1,length_fir);
    
    %% Initial Values of Model Tap
    model_tap = zeros(1,length_fir);
    
    %% System Output where a 40 dB Noise floor is added
    noise_floor = 40;
    d = randn(size(x))*10^(-noise_floor/20);
    sys_opt = filter(sys_desired,1,x)+d;
    
    mu_min = 0.007;
    
    % input variance
    input_var = var(x);
    
    % upper bound = 1/(filter_length * input variance)
    mu_max = 1/(input_var*length_fir);
    
    %% Defining initial parameters for EVSS-LMS algorithm
    mu_EVSS(1) = 0.025;
    mu_EXTRA = 1.5*10^(-4);
    %% EVSS LMS ALGORITHM
    for i=1:length(x)
        
        % model tap values (shifting of tap values by one sample to right)
        model_tap=[x(i) model_tap(1:end-1)];
        
        % model output
        model_out_evss(i) = model_tap * model_coeff_evss';
        
        % error
        e_evss(i) = sys_opt(i) - model_out_evss(i);
        
        %Updating the coefficients
        model_coeff_evss = model_coeff_evss + mu_EVSS(i) * e_evss(i) * model_tap;
        
        if mu_EVSS(i)>mu_min
            c=1;
        else mu_EVSS(i)<=mu_min;
            c=2^(-5); % <= here
            mu_EVSS(i) = mu_min; % <= here (missing line)
        end
        
        if mu_EVSS(i)>mu_max
            mu_EVSS(i) = mu_max; % <= here (missing line)
        end
        
        mu_EVSS(i+1) = c * mu_EVSS(i) + mu_EXTRA * sign(e_evss(i));
        
    end
    %% Storing the e_square values after a whole run of VSS LMS algorithm
    
    err_EVSS(itr,:) = e_evss.^2;
    
    %% Printing the iteration number
    clc
    disp(char(strcat('iteration no : ',{' '}, num2str(itr) )))
end
figure;
plot(10*log10(mean(err_EVSS,1)),'-b');
title('VSS LMS Algorithms'); xlabel('iterations');ylabel('MSE(dB)');
grid on;

figure;
plot(sys_desired,'-*b');
hold on;
plot(model_coeff_evss, '-*m');
legend('FIR model','VSSLMS FIR ID');

[PARVATHY NAIR]

thankyou @Mathieu NOE

[Mathieu NOE]

as always, my pleasure !

[PARVATHY NAIR]

hi @Mathieu NOE

 %% System Output where a 40 dB Noise floor is added
    noise_floor = 40;
    d = randn(size(x))*10^(-noise_floor/20);
    sys_opt = filter(sys_desired,1,x)+d;
    

why 'sys_opt = filter(sys_desired,1,x)+awgn(x,noise_floor)-x ' is changed to the above code?

[Mathieu NOE]

hello

simply because I don't have access to the awgn function (belongs to Communication Toolbox I believe)

but my own code does exactly the same thing (adding white noise)

[PARVATHY NAIR]

thankyou @Mathieu NOE