www.pudn.com > aft.rar > dct_lms.m

function [W, e, Lambda] = dct_lms(u, d, M, alpha, beta, gamma, verbose) 
% function [W, e, Lambda] = dct_lms(u, d, M, alpha, beta, gamma, verbose) 
% 
% dct_lms.m - use DCT-LMS algorithm to estimate optimum weight vectors 
% for linear estimation 
% written for MATLAB 4.0 
% 
% Reference: Ch.10 of Haykin, _Adaptive Filter Theory_, 3rd edition 
% 
% 
% Input parameters: 
%                u       : vector of inputs (real scalars) 
%                d       : vector of desired outputs 
%                M       : final order of predictor 
%                alpha   : base step size for weight updates 
%                beta:   : remembering factor for sliding DCT coefficient updates 
%                gamma   : forgetting factor for estimated eigenvalue updates 
%                verbose : set to i for interactive processing 
% 
% 
% Output parameters: 
%                W       : row-wise matrix of Hermitian transposed weights 
%                          at each iteration 
%                e       : row vector of prediction errors at each time step 
%                Lambda  : row-wise matrix of estimates of process eigenvalues 
%                          at each iteration 
 
% length of maximum number of timesteps that can be predicted 
N = min(length(u),length(d)); 
 
% initialize weight matrix and associated parameters for LMS predictor 
W = zeros(M, N+1); 
Lambda = zeros(M, N); 
 
m = [0:(M-1)]'; 
k = ones(M, 1); k(1) = 1 / sqrt(2); 
W2M = exp(-j * pi / M); 
W2M2 = exp(-j * pi / 2 / M); 
W2Mm = exp(-j * pi * m / M); 
W2M2m = exp(-j * pi * m / 2 / M); 
F1 = W2M2m; 
F2 = conj (W2M2m); 
y = zeros (N, 1); 
e = zeros (N, 1); 
Lambda_n1 = zeros(M, 1); 
 
n=M; 
i = n-M+1:n; 
A1n1 = W2M.^(m*(n-i)) * u(i); 
A2n1 = W2M.^(m* (i-n+2*M-1)) * u(i); 
 
for n = M+1:N 
 
    A1(:, n) = W2M.^m    .*  A1n1 + u(n) * ones(M, 1) - u(n-M) * (-1).^m; 
    A2(:, n) = W2M.^(-m) .* (A2n1 + u(n) * ones(M, 1) - u(n-M) * (-1).^m); 
    C(:, n)  = 1/2 * k .* (-1).^m .* W2M2.^m .* (A1(:, n) + A2(:, n)); 
 
    A1n1 = A1 (: , n); 
    A2n1 = A2(:, n); 
 
    % predict next sample and compute error 
    y(n) = W(:, n).' * C(:, n); 
    e(n) = d(n) - y(n); 
 
    if (verbose ~= 0); 
       disp(['time step ', int2str(n), ': mag. pred. err. = ' , num2str (abs(e(n)))]) 
    end; 
 
% adapt eigenvalue estimate and weight vectors, adjusting for 
% offset of M+i in starting index for eigenvalue step size 
    Lambda(:, n) = gamma * Lambda_n1 + 1 / (n-M) * (C( :, n) .^2 - gamma * Lambda_n1); 
    Lambda_n1    = Lambda (: , n); 
    W(:, n+1)    = W(:, n) + alpha ./ Lambda(:, n) .* C(:, n) * e(n); 
 
end; % for n 
 
% discard last update to W since no more data to predict 
W = W(i:M, i:N); 
 
% end; % dct_lms