fblms.m

This program simulates plant identification using frequency block least mean square (FBLMS) alogrith

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                              %
%  This program simulates plant identification using  frequency block least  mean square (FBLMS) alogrithm      %
%  reference: 《LMS算法的频域快速实现》                                                                         %
%  LMS is modified by XXX in XXX place, see details in XXX relevant document                                   %                                                                       %
%   on Nov 11,2007                                       %
%                                                                                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [w2,y,fblms_mse,N,K] = fblms_1(x,d,num,gamma,u,N);
%generator1: [w2,y,fblms_mse,N,K]=fblms_1(x,d,num,0.9,0.08,32);
%generator2: [w2,y,fblms_mse,N,K]=fblms_1(x,d,num,0.9,0.08,256);
%generator3: [w2,y,fblms_mse,N,K]=fblms_1(x,d,num,0.9,0.08,512);
%generator2: [w2,y,fblms_mse,N,K]=fblms_1(x,d,num,0.9,0.08,1024);
%generator2: [w2,y,fblms_mse,N,K]=fblms_1(x,d,num,0.9,0.08,2048);

%----------------------------------------------------%
% Input parameters:
%            u:   step size
%            N:   number of weights
%            x:   input signal
%            d:   output signal
%            y:   expectation output sequence of the AF
%            e:   the difference sequence
%            K:   experiment time
%            w:   weight
%            #:   edit up to user
%-----------------------------------------------------%

%----------initialization-----------------------------%
K = length(x);                                                %  the length of the input signal,also the iteration times
y = zeros(1,K);                                               %  define the ouput signal of the idetification plant
w1 = zeros(1,2*N);                                            %  define a middle variable of weight vector
e = zeros(1,K);                                               %  define the error signal
W = zeros(1,2*N);                                             %  define the fft transforming variable of w1
Xk = zeros(1,2*N);                                            %  define the fft transforming variable of the input blocks
P1 = ones(1,2*N)/100;                                         %  define the power estimation of input signal 
H=fix(K/N);                                                   %  the total number of the blocks
w2 = zeros(N,H);                                              %  define another middle variable of the weights to store each weight
x1=[zeros(1,N),x,zeros(1,N)];                                 %  extend the input signal in order to cut them to blocks
block=11;                                                     %  the size used to smooth the mse

%------------fblms algorithm begin--------------------%
for i=1:H                                                     %  divide the x into K/N blocks
    Xk(1:2*N) = fft(x1(((i-1)*N+1):(i+1)*N));                 %  carry out the fft processing 2N at a time 
    o1 = real(ifft(Xk.*W));                                   %  get Y(k) by doing ifft to X(k)*Y(k) but only save the last block
    y((i-1)*N+1:i*N) = o1((N+1):(2*N));                       %  discard the first N(only save the last block)(overlap-save method) 
    e((i-1)*N+1:i*N) = d((i-1)*N+1:i*N)-y((i-1)*N+1:i*N);     %  get the error in time domain 
    E = fft([zeros(1,N),e((i-1)*N+1:i*N)]);                   %  insert N zero and tansform it to frequency domain
    P1 = gamma*P1+(1-gamma)*((abs(Xk)).^2);                   %  compute the power estimation
    D =u*(1./P1);                                             %  get the step parameter vector
    o2 = real(ifft(D.*E.*conj(Xk)));                          %  get the gradient estimation
    V = fft([o2(1:N),zeros(1,N)]);                            %  discard the last N and tansform it
    W = W+2*V;                                                %  get the new weight vector in frequency domain
    w1 = real(ifft(W(1:2*N)));                                %  tansform the weight vector to time domain
    w2(:,i)=w1(1:N);                                          %  store each weight value 
end;       
%------------fblms algorithm end----------------------%

%------compute the mean square error------------------%       %  compute the mean square error 
e=e.^2;
E=[zeros(1,fix(block/2)),e,zeros(1,fix(block/2)+1)];
fblms_mse=zeros(1,K);
for i=1:K
    sum_E=0;                                                         
    for j=1:block
    sum_E=sum_E+E(j+i);
    end
    fblms_mse(i)=sum_E/block;
end

%-------plot the weight comparison--------------------%       %  plot the weight comparison
figure(1);
stem(1:N,num,'r');
hold on;grid on;
stem(1:N,w1(1:N),'b');
legend('orginal weights','anolog weights');
title('FBLMS algorithm weights comparison'); 
xlabel('weight number');
ylabel('weight value');

%-------plot the desire and practical signal----------%       %  plot the desire and practical signal
figure(2);
subplot(211);
plot(1:K,d,'b',(1:K),0,'r');
title('desired signal'); 
xlabel('iteration n');
ylabel('magnitude');
subplot(212);
plot(1:K,y,'b',(1:K),0,'r');
title('the practical output signal'); 
xlabel('iteration n');
ylabel('magnitude');

%----------plot the mse-------------------------------%       %  plot the mse
figure(3);
plot(1:K,fblms_mse);
xlabel('iteration n');
ylabel('magnitude');
title('mean square error');

%-------display the mismatch--------------------------%       %  display the mismatch
MIS=10*log10(norm(num-w1(1:N))/(norm(num)));
disp('FBLMS_MIS=');
disp(MIS);

%-------display the misadjustment---------------------%       %  display the misadjustment
i=1:fix(K/N);
power(i)=sum(x((i-1)*N+1:i*N).^2);
MISA(i)=u*power(i)/N;
logmisa=10*log10(MISA);
disp('FBLMS_MISADJUSTMENT=');
disp(MISA(1));
disp('FBLMS_LOG_MISADJUSTMENT=');
disp(logmisa(1));