%OLSBLMS Open-Loop Subband LMS algorithm%%   'ifile.mat' - input file containing:%      Nr - members of ensemble%      N - iterations%      Sx - standard deviation of input%      B - coefficient vector of plant (numerator)%      A - coefficient vector of plant (denominator)%      Sn - standard deviation of measurement noise%      u - convergence factor%      gamma - small constant to prevent the updating factor from getting too large% %      a - small factor for the updating equation for the signal energy in the subbands%      M - number of subbands%      hk - analysis filterbank, with the filter coefficients in the rows%      fk - synthesis filterbank, with the filter coefficients in the rows%      Nw - order of the adaptive filter in the subbands%      %   'ofile.mat' - output file containing:%      MSE - mean-square error%      MSEsub - mean-square error in the subbandsclear all		% clear memoryload ifile;		% read input variablesL=M;                    % interpolation/decimation factorfor l=1:Nr    disp(['Run: ' num2str(l)])    nc=sqrt(Sn)*randn(1,N*M);     % noise at system output     xin=(rand(1,N*M)-.5);    xin=sqrt(Sx)*sqrt(12)*xin;    % input signal    d=filter(B,A,xin);            % desired signal    d=d+nc;                       % unknown system output    % Analysis of desired and input signals    for k=1:M	xaux=filter(hk(k,:),1,d);	dsb(k,:)=xaux(find(mod((1:length(xaux))-1,L)==0)); % desired signal split in subbands	xaux=filter(hk(k,:),1,xin);	xsb(k,:)=xaux(find(mod((1:length(xaux))-1,L)==0)); % input signal split in subbands    end;      for m=1:M	w_ol(m,:)=zeros(1,Nw+1);      % initial coefficient vector for                                      % subband m	x_ol(m,:)=[zeros(1,Nw+1)];    % initial input vector for subband m	sig_ol(m)=0;    end    for k=1:N      if mod(k,200)==0 disp(['Iteration: ' num2str(k)]), end;      for m=1:M 	  % Open-loop	  % --------------------------------------------------	  x_ol(m,:)=[xsb(m,k) x_ol(m,1:Nw)]; % new input vector for subband m	  xaux=shiftdim(x_ol(m,:),1);	  waux=shiftdim(w_ol(m,:),1);	  e_ol(m,l,k)=dsb(m,k)-waux'*xaux;   % error at subband m	  sig_ol(m)=(1-a)*sig_ol(m)+a*(xsb(m,k))^2;	  unlms_ol=2*u/(gamma+(Nw+1)*sig_ol(m));	  w_ol(m,:)=w_ol(m,:)+unlms_ol*e_ol(m,l,k)*shiftdim(xaux,-1); % new coefficient vector for subband m      end;    end;    wk_ol(l,:,:)=w_ol;end;for m=1:M    mse_ol(m,:)=mean(e_ol(m,:,:).^2,2);    % MSE at subband mend;MSE=mean(mse_ol,1);                 % overall MSEMSEsub = mse_ol;save ofile MSE MSEsub;