## Analysis of FxLMS based Spline Adaptive Filtering Algorithm

Version 1.0.4 (58.2 KB) by
The SAF consists of a linear network of adaptive weights in cascade with an adaptive nonlinear network. now we will discuss Hammerstein.
Updated 12 May 2022
In this code, we will identify a nonlinear system using the most recent adaptive Spline filter based filter. These type of filters belong to a class of linear-in-the-parameters nonlinear adaptive filters.
Refer the Notes from
Shivendra Nandan, “Concept of Adaptive Filtering & Spline Adaptive Filtering Algorithm,” DOI: 10.13140/RG.2.2.17453.97767
In which, using Spline filter codes, in this experiment, we have used this filter in a nonlinear system identification scenario, where the nonlinearity is introduced by the some source. For more details, please refer to the technical report. All the parameters used in this code are in coordinance with the parameters used in the above example.
The code is nicely scripted. Hope this helps!!!!!
%MATLAB CODES
%***********************Code 1 ******************************************
% This function evaluates the output of spline nonlinearity
% s is the nonlinearity input s[n]
% af is the nonlinearity structure
% x is the nonlinearity output x[n]
function [ x, af] = ActFunc( s, af)
af.s = s; % Input of the n o n l i n e a r i t y
switch ( af. aftype)
case -1 % Signed Sigmoidal
x = (2*af.Gain/(1+exp(-s*af.Slope))-af.Gain) ;
case 0 % Linear
x = s*af.Slope; % Defaut value
case 1 % Unsigned Sigmoidal
x = af.Gain/(1+exp(-s*af.Slope)) ;
case 2 % Gaussian
x = af.Gain*exp(-s^2/af.Slope ) ;
case 3 % Polynomial
x = 0 ;
for j=1:af.Pord
x = x + af.Q(j)*s.^j; % Sum of monomials
af.g(j) = s.^j;
end
case 20 % Quadratic s p l i n e
np = af.lut_len; % Number of c o n t r o l points
Su = s/af.DeltaX + ( np-1)/2; % F i r s t part of Eq . ( 7 . b )
uIndex = floor(Su) ; % Span index i in Eq . ( 7 . b )
u = Su - uIndex; % Local a b s c i s s a u in Eq . ( 7 . a )
if uIndex<1 % The index must s t a r t from 1
uIndex = 1 ;
end
if uIndex>(np-2)
uIndex = np - 2 ; % The index cannot exceed np - 2
end
af.g = [1 u u^2]*af.C; % F i r s t part of Eq . ( 5 ) : u^T C
x = af.g*af.Q(uIndex : uIndex+2) ; % Eq . ( 5 ) : u^T C q_i
af.uIndex = uIndex; % For d e r i v a t i v e computation
af.uSpline = u; % For d e r i v a t i v e computation
otherwise % Cubic s p l i n e
np = af.lut_len; % Number of c o n t r o l points
Su = s/af.DeltaX + ( np-1)/2; % F i r s t part of Eq . ( 7 . b )
uIndex = floor(Su); % Span index i in Eq . ( 7 . b )
u = Su - uIndex; % Local a b s c i s s a u in Eq . ( 7 . a )
if uIndex<1 % The index must s t a r t from 1
uIndex = 1 ;
end
if uIndex>(np-3) % The index cannot exceed np - 3
uIndex = np - 3;
end
af.g = [u^3 u^2 u 1]*af.C; % F i r s t part of Eq . ( 5 ) : u^T C
x = af.g*af.Q(uIndex : uIndex+3); % Eq . ( 5 ) : u^T C q_i
af.uIndex = uIndex; % For d e r i v a t i v e computation
af.uSpline = u; % For d e r i v a t i v e computation
end
end
%***********************Code 2 ******************************************
% This function implements the adaptation of a Hammerstein SAF (HSAF) structure by using the LMS adaptive algorithm
% Where
% F is the adaptive -lter structure;
% x is the input signal sample x[n];
% d is the desired signal sample d[n];
% y is the output signal sample y[n];
% e is the error signal sample e[n].
function [F,y,e] = AF_LMS_HSPL_F(F,x,d)
M = F.M; % Length of the l i n e a r f i l t e r
F.xw(2:M) = F.xw(1:M-1) ; % S h i f t the input delay-l i n e
[F.xw(1),F.af] = ActFunc(x,F.af); % Load a new input i n t o the delay line
for j=2:M % Constructing the matrix U
if F.af.uIndex == F.af.indexes(j)
F.af.gM(j,:) = F.af.gM(j-1,:);
else
F.af.gM(j,:) = zeros(1,4);
end
end
F.af.gM(1,:) = F.af.g; % Load a new vector in matrix U
F.af.indexes(2:M) = F.af.indexes(1:M-1);
F.af.indexes(1) = F.af.uIndex;
y = F.xw.'*F.w; % F i l t e r output
e = d - y; % Error estimation
% LMS weights and c o n t r o l points update ---------------------------------
F.w = F.w + F.mu*F.xw.*e; % LMS in Eq . ( 1 0 )
if F.af.aftype > 1 % Kind a c t . f . -1 0 1 2 4 5 *
e_av = F.mQ*e;
ii = F.af.uIndex:F.af.uIndex + F.af.P; % P = Spline order
F.af.Q(ii) = F.af.Q(ii) + e_av*( F.w'*F.af.gM).'; % LMS in Eq . ( 1 1 )
end
end
%***********************Code 3 ******************************************
% This function creates and initializes nonlinear function implemented by splines
% where:
% AF is the spline nonlinear structure;
% a-nit is the type of initialization (linear, Gaussian, random,...). In particular:
% -1: signed sigmoid;
% 0: linear function;
% 1: unsigned sigmoid;
% 2: Gaussian;
% 3: random.
% aftype is the type of spline basis (B, Catmull-Rom, Hermite,...). In particular:
% -1: -xed signed sigmoid;
% 0: -xed linear function;
% 1: -xed unsigned sigmoid;
% 2: -xed Gaussian;
% 3: fexible polynomial function;
% 4: fexible Catmull-Rom spline;
% 5: fexible B-spline;
% 6: -exible Bernstein spline;
% 7: -exible parametric --spline;
% 8: -exible Hermite spline;
% 9: -exible Beziér spline;
% DeltaX is the space between knots -x;
% Gain is the x-axis range limits;
% Slope is the initial slope;
% M is the length of the linear -lter M;
% Pord is the polynomial order in case of polynomial nonlinearity Pord
function AF = create_activation_function( afinit, aftype, DeltaX, Gain,Slope, M, Pord)
% Spline a c t i v a t i o n function d e f i n i t i o n and i n i t i a l i z a t i o n --------------
if nargin==0 , help create_activation_function; return;end
if nargin < 7
Pord = 3 ;
if nargin <6
M=1;
if nargin <5
Slope= 2.16 ;
if nargin <4
Gain= 1.1 ;
if nargin <3
DeltaX= 0.4 ;
if nargin <2
aftype=-1;
if nargin <1
afinit=-1;
end
end
end
end
end
end
end
% -----------------------------------------------------------------------
% Check fo n o n l i n e a r i t y type
if afinit == -1
Slope = Slope*2*(1/Gain);
end
if afinit == 1
Slope = Slope*4*(1/Gain);
end
if (( afinit<-1) || ( afinit>3))
fprintf(' Activation function error type\n ');
aftype = -1;
end
% LUT parameters --------------------------------------------------------
Table_Length = G1_TabFuncLen(DeltaX,Gain,Slope,afinit);
lut_len = Table_Length; % Length of LUT
s = 0.0; % Linear combiner output
x = 0.0; % Nonlinearity output
uIndex = 1 ; % Span index
uSpline = 0 ; % Local a b s c i s s a
% Polynomial n o n l i n e a r i t y -----------------------------------------------
if (aftype == 3)
Q = zeros(Pord,1); % Polynomial n o n l i n e a r i t y
Q(1) = 1 ;
g = zeros(1,Pord); % Dot u vector
gM = zeros(Pord,M); % U matrix ( f o r Hammerstein f i l t e r ) in [ 2 ]
else
Q = zeros(lut_len,1) ; % Look-Up Table n o n l i n e a r i t y
end
% I f NOT s p l i n e ---------------------------------------------------------
if (aftype < 4)
C = 0 ;
end
% Catmul-Rom s p l i n e n o n l l i n e a r i t y ---------------------------------------
if (aftype == 4)
C = 0.5*[-1 3 -3 1 ;
2 -5 4 -1;
-1 0 1 0; % Page 775 , top of column 1
0 2 0 0];
end
% B s p l i n e n o n l l i n e a r i t y ------------------------------------------------
if (aftype == 5)
C = (1/6)*[-1 3 -3 1;
3 -6 3 0; % Page 774 , bottom of column 2
-3 0 3 0;
1 4 1 0 ];
end
% Bernstein polynomial n o n l i n e a r i t y -------------------------------------
if (aftype == 6)
C = [-1 3 -3 1;
3 -6 3 0; % Bernstein polynomials
-3 3 0 0;
1 0 0 0];
end
% Parametric s p l i n e n o n l i n e a r i t y ----------------------------------------
if (aftype == 7)
tau = 0.5; % tau = 0 . 5 ==> CR-s p l i n e
C = [-tau 2-tau tau-2 tau;
2*tau tau-3 3-2*tau -tau; % Parametric s p l i n e
-tau 0 tau 0;
0 1 0 0];
end
% Hermite s p l i n e n o n l i n e a r i t y -------------------------------------------
if (aftype == 8)
C = [2 -2 1 1;
-3 3 -2 -1; % Hermite polynomials
0 0 1 0;
1 0 0 0];
end
% Bezier spline no n l l i n e a r i t y
if (aftype == 9)
C = (1/6)*[1 3 -3 1;
3 -6 3 0; % Bezier polynomials
-3 3 0 0;
1 0 0 0];
end
% Quadratic B-s p l i n e n o n l i n e a r i t y ---------------------------------------
if (aftype == 20)
C = [ 1 1 0;
-2 2 0;
1 -2 1];
end
% Spline order
P = length(C) - 1 ;
% Dummy arrays f o r l ea r n ing algorithms . See Eqs . ( 5 ) and ( 6 ) ------------
if (aftype > 3)
g = zeros(1,P+1); % The dot u vector
gM = zeros(M,P+1); % The U matrix ( f o r Hammerstein SAF) in [ 2 ]
end
indexes = ones(M,1);
% STRUCTURE DEFINITION
AF = struct('aftype',aftype,'afinit',afinit,'s',s,'x',x,'Q',Q,'lut_len',lut_len,'uIndex',uIndex,'uSpline',uSpline,'Slope',Slope,'Gain',Gain,'DeltaX',DeltaX,'P',P,'Pord',Pord,'C',C,'g',g,'gM',gM,'indexes',indexes);
% -----------------------------------------------------------------------
% For s p l i n e i n t e r p o l a t i o n ----------------------------------------------
if ( aftype>1 && aftype ~= 3)
LutSlope = (Table_Length - 1)/2.0; % New slope
X = -LutSlope*DeltaX;
for j=1 : Table_Length % Table_Length
AF.Q( j) = G1_FUNC( X, Gain, Slope, afinit) ;
X = X + DeltaX;
end
end
end
%***********************Code 4 ******************************************
% This function creates and initializes the structure implementing a WSAF architecture and adapted by the LMS adaptive algorithm.
% Where
% lms_af is the adaptive -lter structure;
% M is the length of the linear -lter M;
% mu is the step-size of the linear -lter µ;
% mQ is the step-size of the nonlinearity µQ;
% delta is the regularization parameter -;
% af is the spline nonlinearity structure.
function lms_af = create_SPL_lms_adaptive_filter_1( M, mu, mQ, delta, af)
% LMS SAF adaptive f i l t e r d e f i n i t i o n and i n i t i a l i z a t i o n -----------------
w = zeros(M,1); % Linear f i l t e r taps
w(floor(M/2)) = 1 ; % Adaptive f i l t e r weights i . c .
dI = delta; % Regularizing parameter
xd = zeros(M,1) ; % Buffer of the desired s i g n a l
xw = zeros(M,1) ; % Buffer of the f i l t e r s t a t u s
lms_af = struct('M',M,'mu',mu,'mQ',mQ,'dI',dI,'w',w,'xd',xd,'xw',xw,'af',af);
end
%***********************Code 5 ******************************************
% This function evaluates the output of a Hammerstein nonlinear structure implemented by a HSAF architecture
% where
% F is the adaptive -lter structure;
% x is the input signal sample x[n];
% y is the output signal sample y[n];
% s is the linear combiner input array sn.
function [F,y,s] = FW_HSPL_F(F,x)
M = F.M; % Length of the l i n e a r f i l t e r
F.xw(2:M) = F.xw(1:M-1) ; % S h i f t of the input delay-l i n e
F.xw(1) = x; % Load a new input i n t o the delay-l i n e
s = zeros (M,1); % Linear combiner input array
for i=1:M
[s(i),F.af] = ActFunc(F.xw(i),F.af); % Evaluating the n o n l i n e a r i t y output
end
y = s'.*F.w; % F i l t e r output
end
%***********************Code 6 ******************************************
% Nonlinear Function Implementation
function value = G1_FUNC( X, G, S, ty)
switch ty
case -1
value = 2*G/(1+exp(-X*S)) - G; % Signed Sigmoid
case 0
value = X*S; % Linear
case 1
value = G/(1+exp(-X*S)); % Unsigned Sigmoidal
case 2
value = G* exp(-(X*X)/5.0) ; % Gaussian
case 3
value = 2*G*(rand - 0.5); % Random
otherwise
value = 0 ;
end
end
%***********************Code 7 ******************************************
% For Calculation of Table Length
function Table_Length = G1_TabFuncLen( DX, G, S, ty)
ii = 0 ;
X = 0 ;
if (ty~=2 && ty~=3)
F = 0.0 ;
crtGain = G - 0.005*G; % Max a t c func value ---------------------
while ( F<crtGain)
F = G1_FUNC( X, G, S, ty) ;
X = X + DX;
ii = ii + 1 ;
end
elseif (ty==3)
ii = 11 ;
else
crtGain = 0.005*G; % Gaussian
F = G;
while ( F>crtGain )
ii = ii + 1 ;
F = FUNC( X, G, S, ty) ;
X = X + DX;
end
end
Table_Length = ii*2 + 1; % always odd
end
%***********************Code 8 ******************************************
% Main File of Hammerstein Spline Adaptive Filter (HSAF)
% Implements a convergence test of a Hammerstein spline adaptive -lter (HSAF).
clc
clear
close all
% -----------------------------------------------------------------------
%% Parameters s e t t i n g
% Input parameters ------------------------------------------------------
Lx = 30000; % Length of input s i g n a l
nRun = 10 ; % Number of runs
out_noise_level_dB = 60 ; % SNR
out_noise_level = 10^(-out_noise_level_dB/20) ; % Noise l e v e l
x = zeros(Lx,1) ; % x (Nx1) input s i g n a l array d e f i n i t i o n
% Colored s i g n a l generation ---------------------------------------------
a = 0.1;
b = sqrt(1-a^2) ;
%x = f i l t e r ( b , [1 -a ] , randn ( s i z e ( x ) ) ) ; % H( z ) = b/(1+a * z^-1)
disp ( ' . . . . . . ' ) ;
% Adaptive f i l t e r d e f i n i t i o n --------------------------------------------
M = 7 ; % Length of l i n e a r f i l t e r
mu0 = 0.1 ; % Learning r a t e f o r l i n e a r f i l t e r
mQ0 = 0.1 ; % Learning r a t e f o r c o n t r o l points
if Lx < 30000 % Batch f o r evaluating MSE
B = 100 ;
else
B = 4000 ;
end
% Spline a c t i v a t i o n function d e f i n i t i o n and i n i t i a l i z a t i o n --------------
afinit = 0 ; % I n i t a c t . func . -1 0 . . . (ONLY -1, bip . s i g . or 0 =linear)
aftype = 4 ; % Kind a c t . f . -1 0 1 2 4 5 ; (4 = CR-spline , 5 = B-s p l i n e )
Slope = 1 ; % Slope
DeltaX = 0.2 ; % Delta X
x_range = 2 ; % Range l i m i t
% Creating the n o n l i n e a r i t y ---------------------------------------------
af0 = create_activation_function( afinit, aftype, DeltaX, x_range, Slope,M); % Model
af1 = create_activation_function( afinit, aftype, DeltaX, x_range, Slope,M); % SAF
%% I n i t i a l i z a t i o n
% --- Target D e f i n i t i o n -------------------------------------------------
TH1 = create_SPL_lms_adaptive_filter_1(M,mu0,mQ0,1e-2,af0); % Target h Model LMS
% TARGET: Nonlinear memoryless function implemented by Spline interpolated LUT
Q0 = [ -2.20
-2.00
-1.80
-1.60
-1.40
-1.20
-1.00
-0.80
-0.91
-0.40
-0.20
0.05
0.0
-0.40
0.58
1.00
1.00
1.20
1.40
1.60
1.80
2.00
2.20
] ;
TH1.af.Q = Q0;
QL = length ( Q0) ; % Number of c o n t r o l points
% Linear f i l t e r ---------------------------------------------------------
TH1.w = [ 0.6 -0.4 0.25 -0.15 0.1 -0.05 0.001 ]; % MA system to be identified
% --- SAF d e f i n i t i o n ----------------------------------------------------
H1 = create_SPL_lms_adaptive_filter_1(M,mu0,mQ0,1e-2,af1); % HSAF LMS
% I n i t i a l i z e ------------------------------------------------------------
N = Lx + M + 1 ; % Total samples
for i = Lx+1:N
x(i) =0;
end
dn = zeros(N,1) ; % Noise desired output array
d = zeros(N,1) ; % Desired s i g n a l array
y = zeros(N,1) ; % Output array
e = zeros(Lx,1) ; % Error array
em = zeros(Lx,1) ; % Mean square e r r o r
varW = zeros(M,1) ; % Variance value of w
qm = zeros(QL,1) ; % Mean value Spline c o e f f
varQ = zeros(QL,1) ; % Variance value Spline c o e f f
%% Main loop ------------------------------------------------------------
disp (' Algorithm start . . . ') ;
t = clock ;
for n = 0 : nRun-1
fprintf( ' Test nr . %d/%d\n ' , n+1 , nRun) ;
x = filter( b, [1 -a] , randn(size(x))) ; % H( z ) = b/(1+a * z^-1)
dn = out_noise_level * randn(size(x)); % Noise
% SAF I .C. ----------------------------------------------------------
H1.w (:) = 0 ;
H1.w ( 1 ) = 0.1 ;
% Set Activation Func I .C. ------------------------------------------
H1.af.Q = af0.Q;
% HSAF Evaluation ---------------------------------------------------
for k = 1 : Lx
% Computing the desired output ----------------------------------
%[TH1,d(k),snk] = FW_HSPL_F(TH1,x(k)); % Hammerstein model
% Updating HSAF -------------------------------------------------
[H1,y(k),e(k)] = AF_LMS_HSPL_F(H1,x(k),d(k)+ dn(k)) ; % SAF LMS ( Eqs .(10) and (11))
end
em = em + (e.^2 ) ; % Squared e r r o r
wm = n+111;
% SAF run-time mean and variance estimation -------------------------
wm = (1/(n+1))*H1.w + (n/(n+1))*wm;
varW = varW + (n/(n+1))*((TH1.w - wm).^2) ;
qm = (1/(n+1))*H1.af.Q+(n/(n+1))*qm;
varQ = varQ + ( n/(n+1) )*((TH1.af.Q - qm).^2);
end
em = em/nRun; % MSE
H1.af.Q = qm;
%------------------------------------------------------------------------
% Average MSE evaluation
mse = mean(em(end-B-M-1:end-M-1) ) ; % Average MSE
%------------------------------------------------------------------------
fprintf('\n') ;
%% Results
% -----------------------------------------------------------------------
% P r i n t t a b l e of means and variances
% -----------------------------------------------------------------------
fprintf('\n') ;
fprintf( 'Number of iterations = %d\n ' ,nRun) ;
fprintf( ' Learning rates : muW = %5.3 f muQ = %5.3 f \n ' , mu0, mQ0) ;
fprintf( ' a = %4.2f b = %4.2f\n',a, b ) ;
fprintf( 'Number of filter weights = %d\n ' , M) ;
fprintf( 'AF type = %d\n ' ,aftype) ;
fprintf( ' DeltaX = %4.2 f \n ' ,DeltaX) ;
fprintf( 'SNR_dB = %4.2 f dB\n ' ,out_noise_level_dB) ;
fprintf( ' Steady-state MSE = %5.7f , equal to %5.3f dB\n ' ,mse,10*log10 (mse)) ;
fprintf('\n ') ;
fprintf('Mean and Variance Tables ----------------------------------\n ' ) ;
for i=1:QL
fprintf('i=%2d q0 =%5.2 f qm =%9.6f varQ = %10.3e \n ',i,TH1.af.Q( i),qm( i),varQ(i));
end
fprintf('\n');
fprintf('-----------------------------------------------------------\n');
for i=1:M
fprintf('i=%d w0 =%5.2f wm =%9.6f varW = %10.3e \n ',i,TH1.w(i) , wm( i) , varW( i) ) ;
end
fprintf('-----------------------------------------------------------\n ' ) ;
% -----------------------------------------------------------------------
% P l o t t i n g f i g u r e s
% -----------------------------------------------------------------------
% Plot Spline functions -------------------------------------------------
yLIM = 1.5 ;
xLIM = 3.0 ;
figure1 = figure('PaperSize',[20.98 29.68]);
box('on');
hold on;
hold('all');
ylim([-yLIM yLIM]);
xlim([-xLIM xLIM]);
grid on;
KK = 500;
yy1 = zeros(1,KK);
yy2 = zeros(1,KK);
xa1 = zeros(1,KK);
dx = 2*xLIM/KK;
xx = -xLIM;
for k = 1 : KK
yy1(k) = ActFunc(xx, TH1.af) ; % Model
yy2( k) = ActFunc(xx, H1.af) ; % Adapted
xa1( k) = xx;
xx = xx + dx;
end
xlabel('Input','FontSize',12,'FontWeight','demi');
ylabel('SAF output','FontSize',12,'FontWeight','demi');
title('Profile of model and adapted nonlinearity','FontSize',12,'FontWeight','demi');
plot( xa1,yy1,'-g','LineWidth',2);
plot( xa1,yy2,'--','LineWidth',2);
set( gca,'FontSize',10,'FontWeight','demi');
% Filter coefficients ---------------------------------------------------
figure2 = figure('PaperSize',[20.98 29.68]) ;
hold on;
plot( TH1.w, 'LineWidth' , 2 ) ;
plot( H1.w, 'm' , 'LineWidth' , 2 ) ;
xlabel('time {\itn} ' , 'FontSize' , 12 , 'FontWeight' , 'demi' ) ;
ylabel('Linear combiner coefficients' , 'FontSize' , 12 , 'FontWeight' , 'demi' ) ;
set(gca, 'FontSize' , 10 , 'FontWeight' , 'demi' ) ;
% MSE dB ----------------------------------------------------------------
figure3 = figure('PaperSize' , [10 15]);
box('on');
hold on;
hold('all');
ylim([-out_noise_level_dB-5 10]) ;
grid on;
edb = 10*log10( em );
[bb,aa] = butter(2,0.02) ;
plot(filter(bb,aa,edb), 'Color' ,[1 0 0], 'LineWidth' ,2);
noiseLevel(1 : length(edb)-1) = -out_noise_level_dB;
plot(noiseLevel, '--' , 'Color' ,[0 0 1], 'LineWidth' ,2);
title('Hammerstein SAF convergence test' , 'FontSize' ,12, 'FontWeight','demi');
xlabel('Samples','FontSize',12,'FontWeight','demi');
ylabel('MSE','FontSize',12,'FontWeight','demi');
legend ('MSE','NoiseLevel');
set(gca ,'FontSize' , 10 , 'FontWeight' , 'demi');
set( gcf , 'PaperSize' , [20.98 29.68]);
% -----------------------------------------------------------------------
fprintf( 'END Hammerstein Spline Adaptive Filter (HSAF) ---------------------------------------------\n ' ) ;

### Cite As

Shivendra Nandan (2024). Analysis of FxLMS based Spline Adaptive Filtering Algorithm (https://github.com/shivendranandan/Matlab_codes), GitHub. Retrieved .

##### MATLAB Release Compatibility
Created with R2022a
Compatible with any release
##### Platform Compatibility
Windows macOS Linux

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Version Published Release Notes
1.0.4

New version

1.0.3