cfdlms.m

《自适应滤波算法与实现》（第二版）源码

%CFDLMS Constrained Frequency-Domain 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
%      L - decimation/interpolation factor
%      Nw - order of the adaptive filter in the subbands
%      
%   'ofile.mat' - output file containing:
%      MSE - mean-square error


clear all		% clear memory
load ifile;		% read input variables

M=2*L;                  % number of bins

% Constrained Frequency-Domain LMS
for l=1:Nr
    disp(['Run: ' num2str(l)])
    nc=sqrt(Sn)*randn(1,N*M);     % noise at system output 
    xin=rand(1,N*M);
    xin=sqrt(12)*(xin-.5);        % input signal
    d=filter(B,A,xin);            % desired signal
    d=d+nc;                       % unknown system output

    y=zeros(M,N);
    xsb=zeros(M,N);
    dsb=zeros(L,N);
    Zy=zeros((L)*(Nw+1)-1,1);
    sig=zeros(M,1);
    ww=zeros(M,(Nw+1));
    wwc=zeros((Nw+1),(M));
    w=reshape(ww,Nw+1,M);
    uu=zeros(M,Nw+1);
    xin=[zeros(1,L) xin];
    for k=1:N
      % Serial-to-parallel converter
      aux1=[1:L]+(k-1)*L;
      aux2=[L+1:M]+(k-1)*L;
      x_p=[fliplr(xin(aux2)) fliplr(xin(aux1))]; xsb(:,k)=x_p';
      aux=[1:L]+(k-1)*L;
      d_p=fliplr(d(aux)); dsb(:,k)=d_p';
      ui=fft(xsb(:,k))/sqrt(M);
      uu=[ui , uu(:,1:end-1)];
      for m=1:M
          uy(m,:)=uu(m,:)*ww(m,:).';
      end;
      y(:,k)=ifft(uy)*sqrt(M);
	     
      e=(dsb(:,k))-(y(1:L,k));		% error samples
      out(aux)=y(1:L,k);		% adaptive filter output
      IL=[eye(L) ; zeros(M-L,L)];
      et=fft([e;zeros(M-L,1)])/sqrt(M); % F*IL*e;
      clear wwc;
      for m=1:M
	  sig(m)=(1-a)*sig(m)+a*abs(ui(m)).^2;
	  wwc(m,:)= u/(gamma+(Nw+1)*sig(m))*conj(uu(m,:))*(et(m)); % new increment for the coefficient vector for subband m (unconstrained)
      end;
      % Constraint
      waux=fft(wwc)/sqrt(M);
      wwc=sqrt(M)*ifft([waux(1:L,:); zeros(M-L,Nw+1)]); % new increment for the coefficient vector for subband m (constrained)
      ww=ww+wwc; 	% new coefficient vector for subband m (constrained)
      e_fda(:,k)=abs(e).^2;
    end;
    e_fd(l,:)=e_fda(:)';
end;


MSE=mean(abs(e_fd),1);
MSE=shiftdim(MSE,1);

save ofile MSE