adaptive_interference_cancellation_with_lms_algorithm.m

this is a code for adaptive interfrence cancellation on a certain signal using the least mean square

% Adaptive Interference Cancellation
clear all
close all
f_sampling=10e3;
number_of_points=2000;
time_instants=0:number_of_points-1;
% Interfering signal has fluctuating frequency
frequency_of_interfering_signal=50+5*cos(2*pi*rand+2*pi*10*time_instants/f_sampling)+0.1*(rand(size(time_instants))-0.5); %1x2000
% Signal plus interference plus noise
signal_plus_interference=1.1*sin(2*pi*20*time_instants/f_sampling)+...
                         sin(2*pi*100*time_instants/f_sampling)+...
                         0.7*sin(2*pi*150*time_instants/f_sampling)+...
                         0.8*sin(2*pi*200*time_instants/f_sampling)+...
                         2*sin(2*pi*frequency_of_interfering_signal.*time_instants/f_sampling)+...
                         0.2*randn(1,number_of_points);

% Input to LMS filter
input_to_LMS_filter=sin(2*pi*frequency_of_interfering_signal.*time_instants/f_sampling);

% LMS filter parameters
number_of_weights=20;
step_size=0.005;
weight_vector=zeros(1,number_of_weights); % Initialize with zeros

% Static notch filter
edge_frequencies=[1000 1500];
% Use fir1 instruction to design a bandstop filter that suppresses the
% frequency components from 45 to 55 Hz
% The filter order should be number_of_weights
% Use filter command to find the output of the filter
% Call output output_of_static_notch_filter
% ENTER YOUR CODE HERE
normalized_edge_frequencies = edge_frequencies/5000; 
num_of_notch_filter = fir1(20,normalized_edge_frequencies,'stop'); % fir filters have no poles thats why the output is a vector 1x20 --- always stable
den_of_notch_filter = 1; 
output_of_static_notch_filter = filter(num_of_notch_filter,den_of_notch_filter,signal_plus_interference);% filtered signal by notch filter

% Running LMS algorithm
% Initialize the error signal which is used in this configuration as the
% output we are interested in
error(1:number_of_weights-1)=signal_plus_interference(1:number_of_weights-1); 
for index=number_of_weights:number_of_points
    % Implement LMS algorithm
    % You should complete the error vector afkter this loop is finished
    % Do not use Matlab direct LMS instructions
    % ENTER YOUR CODE HERE
    output_of_adaptive_filter = conj(input_to_LMS_filter(index-19:index))*weight_vector.'  ; % outputs 1x1 element y(n) generally for a complex signal
    delta_weight = 2*step_size*input_to_LMS_filter(index-19:index)*(signal_plus_interference(index)-output_of_adaptive_filter); % weight update equation
    weight_vector = weight_vector + delta_weight; % recursive formula for modifying the weight vector 1x20
    error = [error (signal_plus_interference(index)-output_of_adaptive_filter)]; % concatenation between the error vector and the new error
    % from each iteration
end

% Compute the magnitude spectrum, i.e., abs(fft(.)), of the signal with intereference and the
% filtered signal with LMS and notch filter
fft_size=2^13;
frequency_axis=0:f_sampling/fft_size:f_sampling/20;
% The three variables should be called
% magnitude_spectrum_of_signal_plus_interference --> signal plus inference
% magnitude_spectrum_of_filtered_signal --> LMS processed signal
% magnitude_spectrum_of_output_of_static_notch_filter --> output of static
% notch filter
% ENTER YOUR CODE HERE
magnitude_spectrum_of_signal_plus_interference = abs(fft(signal_plus_interference,fft_size));
magnitude_spectrum_of_signal_plus_interference = magnitude_spectrum_of_signal_plus_interference(1:410); % length is taken according to the specified freq. axis
magnitude_spectrum_of_filtered_signal = abs(fft(error,fft_size));
magnitude_spectrum_of_filtered_signal = magnitude_spectrum_of_filtered_signal(1:410);
magnitude_spectrum_of_output_of_static_notch_filter = abs(fft(output_of_static_notch_filter,fft_size));
magnitude_spectrum_of_output_of_static_notch_filter = magnitude_spectrum_of_output_of_static_notch_filter(1:410);
% Time domain plots
subplot(2,1,1)
plot(time_instants,signal_plus_interference)
grid
title('Signal plus interference')
subplot(2,1,2)
plot(time_instants,error)
grid
xlabel('Time (samples)')
title('Output of LMS filter')
% Frequency domain plots
figure
subplot(3,1,1)
plot(frequency_axis,magnitude_spectrum_of_signal_plus_interference/max(magnitude_spectrum_of_signal_plus_interference))
grid
title('Signal plus interference')
subplot(3,1,2)
plot(frequency_axis,magnitude_spectrum_of_filtered_signal/max(magnitude_spectrum_of_filtered_signal))
grid
ylabel('Normalized magnitude spectrum')
title('Output of LMS filter')
subplot(3,1,3)
plot(frequency_axis,magnitude_spectrum_of_output_of_static_notch_filter/max(magnitude_spectrum_of_output_of_static_notch_filter))
grid
xlabel('Frequency (Hz)')
title('Output of static notch filter')
%%
% LPC
clear all
close all
load data_sequence
% Estimate the correlation function of signal x
% You need R(0) to R(15)
% You may use the xcorr command.
% YOUR CODE HERE
c = xcorr(x,15,'biased'); 

% Form matrices R and r of the Yule Walker equation
% you may use toeplitz command
% Use matrix inversion to get coefficients
% Call them lpc_coefficients
% YOUR CODE HERE
r = c(17:31);
R =toeplitz(c(16:30));
R_inverse = inv(R);
lpc_coefficients =R_inverse*r' ;
coefficients_from_Matlab=-lpc(x,15); 
coefficients_from_Matlab=coefficients_from_Matlab(2:16)';
if (norm(lpc_coefficients-coefficients_from_Matlab)<1e-6)
    display('Your result is the same as using Matlab lpc command')
else
    display('There is a difference between your coefficients and those produced by Matlab lpc command')
end