hello
well, it is a bit converging with your data , but it's not "as good" as with a truly stationnary and ergodic signals we prefer to use for system identification (white noise , random signal, which is the same)
why use seismic data for that task ???
also note that in the original code , the random noise input is updated at each iteration while with your signal we can only run 1 iteration; it makes the reading of the convergence much more difficult as there is no averaging effect
your signal has also a quantization problem : the signal is quite "staircase" as the max peak of signal is only coded on 10 levels
this is also a limiting factor of LMS algorithm performance
try to see if you can get a better signal with a analog preamp before signals gets digitized
With random signal , the LMS and VSSLMS identified models are perfectly overlaid (because we have the "best" input signal)
with your signal and one single iteration, there is still a convergence but not complete; you can try to use a longer excitation signal or increase the mu gain
this plot is obtained with the same convergence gains as in the original code :
now we can make the convergence much faster : i increased the gain for both algos by a factor 100
%learning rate for LMS algorithm
mu_LMS = 0.0004*100;
and
%% Defining initial parameters for VSS-LMS algorithm
mu_VSS(1)=1; %initial value of mu for VSS
alpha = 0.97;
gamma = 4.8e-4*100;
and we got almost the "perfect" results even though the signal is not really 100% appropriate for the tasks
and the convergence curve :
this is in the updated code below :
% Desired System (link: http://www.firsuite.net/FIR/AKSOY08_NORCHIP_G40)
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);
%% Ergodic Process
% Note that the uploaded Figure is obtained for itr=500
x = readmatrix('BKHL.HHZ.new.csv');
x = x./max(abs(x)); % normalize the data to +/-1 range (optionnal)
x = x';
for itr=1:10%:2
%% Defining input and initial Model Coefficients
%input
%x=randn(1,60000);
% Model for the LMS Algorithm
model_coeff = zeros(1,length(sys_desired));
% Model for the VSS-LMS Algorithm
model_coeff_vss = zeros(1,length(sys_desired));
%% Initial Values of Model Tap
model_tap = zeros(1,length(sys_desired));
%% 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)+awgn(x,noise_floor)-x;
sys_opt = filter(sys_desired,1,x)+d;
%% Lower and Upper Bounds of Learning Rate
%%%%%%%% These above informations can be accessed from this paper:
%%% R. H. Kwong and E. W. Johnston, "A variable step size LMS algorithm," in IEEE Transactions on Signal Processing, vol. 40, no. 7, pp. 1633-1642, July 1992.
%%%doi: 10.1109/78.143435
%%%%%%%%%%%%%%%%%%
%input variance
input_var = var(x);
% upper bound = 1/(filter_length * input variance)
mu_max = 1/(input_var*length(sys_desired));
%learning rate for LMS algorithm
mu_LMS = 0.0004*100;
% lower bound = learning rate for LMS algorithm
mu_min = mu_LMS;
%%%%%%%%%%%%%%% NOTE %%%%%%%%%%%%%%%%%%%%%%%%%
% If the difference between mu_max and mu_min is not sufficient
% Error Curves for both LMS and VSS LMS algorithms would be IDENTICAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Defining initial parameters for VSS-LMS algorithm
mu_VSS(1)=1; %initial value of mu for VSS
alpha = 0.97;
gamma = 4.8e-4*100;
%Note that these values were defined in the above paper
%% Execution of Both LMS and VSS-LMS algorithms
for i=1:length(x)
%% LMS Algorithm
% 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(i) = model_tap * model_coeff';
%error
e_LMS(i)=sys_opt(i)-model_out(i);
%Updating the coefficients
model_coeff = model_coeff + mu_LMS * e_LMS(i) * model_tap;
%% VSS LMS Algorithm
% model output
model_out_vss(i) = model_tap * model_coeff_vss';
% error
e_vss(i) = sys_opt(i) - model_out_vss(i);
%Updating the coefficients
model_coeff_vss = model_coeff_vss + mu_VSS(i) * e_vss(i) * model_tap;
%Updating the mu value using the VSS algorithm
mu_VSS(i+1) = alpha * mu_VSS(i) + gamma * e_vss(i) * e_vss(i) ;
%% Checking the constraints of mu as given in the paper
if (mu_VSS(i+1)>mu_max)
mu_VSS(i+1)=mu_max;%max
elseif(mu_VSS(i+1)<mu_min)
mu_VSS(i+1)= mu_min;
else
mu_VSS(i+1) = mu_VSS(i+1) ;
end
end
%% Storing the e_square values after a whole run of LMS and VSS LMS algorithms
err_LMS(itr,:) = e_LMS.^2;
err_VSS(itr,:) = e_vss.^2;
%% Printing the iteration number
clc
disp(char(strcat('iteration no : ',{' '}, num2str(itr) )))
end
%% Comparing the Error Curves
figure;
plot(10*log10(mean(err_LMS,1)),'-b');
hold on;
plot(10*log10(mean(err_VSS,1)), '-r');
title('Comparison of LMS and VSS LMS Algorithms'); xlabel('iterations');ylabel('MSE(dB)');legend('LMS Algorithm','VSS LMS Algorithm')
grid on;
figure;
plot(sys_desired,'-*b');
hold on;
plot(model_coeff, '-*r');
plot(model_coeff_vss, '-*m');
legend('FIR model','LMS FIR ID','VSSLMS FIR ID');
