vsnlms.m

刚写完的VSNLMS程序

function [input_signal, error_signal, desired_signal, filter_output, impulse, filter_current, mse, db, db_avg]=VSNLMS(filter_size, p, ss_lower, input_file, iterations)
% Function to perform the VSNLMS algorithm on an input file.
% Inputs:   Filter order, step size, stepsize upper bound, stepsize lower bound, 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, mean sqaure error, attenuation (dB), average attenuation.

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

% generate impulse response and convolve to get echoed signal
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 variables
filter_current = zeros(filter_size,1);
input_vector = zeros(filter_size, 1);
stepsize_curr = zeros(filter_size, 1);
stepsize_prev = zeros(filter_size, 1);
gradient_curr = zeros(filter_size, 1);
gradient_prev = zeros(filter_size,1);

% perform VSNLMS algorithm for set number of iterations
for i=1:iterations
    i
    % get input sample
    input_vector(1)=input_signal(i);
    filter_output(i)=dot(filter_current, input_vector);
    error= desired_signal(i)-filter_output(i)
    gradient_curr=error*input_vector;
    % calcualte variable upper limit
    ss_upper = 1/(2*(dot(input_vector,input_vector)));
    % update stepsize vector
    stepsize_curr = stepsize_prev + p*gradient_curr.*gradient_prev;
    
    % Maintain step_size vector between upper and lower limits.
    for j=1:filter_size
        if stepsize_curr(j) > ss_upper;
            stepsize_curr(j) = ss_upper;
        end
        if stepsize_curr(j) < ss_lower
            stepsize_curr(j) = ss_lower;
        end            
    end
    % update filter tap weights 
    filter_current = filter_current + 2*diag(stepsize_curr)*gradient_curr;
    
    % shift values in input vector along ready for next iteration
    for j=filter_size:-1:2
        input_vector(j)=input_vector(j-1);    
    end
    
    %store values of error and cost function.
    error_signal(i)=error;
    cost(i)=error*error;
    
    %update variable vectors for next iteration
    gradient_prev=gradient_curr;
    stepsize_curr=stepsize_prev;
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);