rls0.m

一种改进的RLS算法

function [input_signal, error_signal, desired_signal, filter_output, impulse, filter_prev, cost, mse, db,db_avg, fc0, fc1,fc2, fc3,fc4]=RLS(filter_size, lambda, input_file, iterations)
% Function to perform the RLS algorithm on an input file.
% Inputs:   Filter order, lambda (forgetting factor), input wav file, number of iterations.
% Outputs:  Input signal, error estimation signal (echo cancelled), desired signal (echoed signal), adaptive filter output, real impulse response
%           Estimation of impulse response, cost, mean sqaure error, attenuation (dB), average attenuation.

% Read in input file.
input_signal = wavread(input_file);

% Generate desired signal by creating impulse and convolving
impulse=zeros(filter_size,1);
for (i=1:5)
    impulse(((i-1)*filter_size/5)+1)=1/i;    
end
desired_signal = conv(input_signal, impulse);

% Initialise varaibles
filter_prev = zeros(filter_size,1);
input_vector = zeros(filter_size, 1);
psi_inv_prev = eye(filter_size);
intermediate= zeros(filter_size, 1);
gain = zeros(filter_size, 1);

% Perform RLS algorithm for set number of iterations.
for i=1:iterations
    i
    input_vector(1)=input_signal(i);
    intermediate = psi_inv_prev*input_vector;
    gain = (1/(lambda+dot(input_vector, intermediate)))*intermediate;
    filter_output(i)=dot(filter_prev, input_vector);
    error= desired_signal(i)-filter_output(i)
    filter_prev = filter_prev + gain*error;
    psi_inv_prev = (1/lambda)*(psi_inv_prev - gain*((input_vector')*psi_inv_prev));
    % Shift along input vector
    for j=filter_size:-1:2
        input_vector(j)=input_vector(j-1);    
    end
    error_signal(i)=error;
    cost(i)=error*error;
    if (i==1)
        fc0=filter_prev;
    end
    if (i==2500)
        fc1=filter_prev;
    end
    if (i==5000)
        fc2=filter_prev;
    end
    if (i==7500)
        fc3=filter_prev;
    end
    if (i==10000)
        fc4=filter_prev;
    end
   
end

% Find moving average of error squared.
for i=1:iterations-100
    mse(i)=mean(cost(i:i+100));
    i
end

%find moving avarage of db attenuation (averaged to smooth output).
for i=1:iterations-2500
    db(i)=-20*log10(mean(abs(desired_signal(i:i+2500)))'./mean(abs(error_signal(i:i+2500))));
    i
end

%find total avarage db attenuation
db_avg=mean(db);