Here is a script-file that simulates the problem:
function  = noise_cancellation()
close all; clear all; clc; set(0, 'DefaultFigureWindowStyle', 'normal');
rng('default'); h_eeg = @(signal_length) -0.2 + (0.4).*rand(1, signal_length);
h_noise_source = @(indices, f, fs) 1.2 .* sin(2.*pi.*f./fs.*indices);
h_noise_contaminate = @(indices, f, fs) 0.12 .* sin(2.*pi.*f./fs.*indices + pi/2);
eeg_signal = h_eeg(signal_length);
noise_signal = h_noise_source(1:signal_length, f, fs);
noise_added = h_noise_contaminate(1:signal_length, f, fs);
eeg_signal_contaminated = eeg_signal + noise_added;
%% Two-tap LMS Adaptive Filter
[ tapWeights, error ] = LMS_MRW( stepSize, tapNumber, reference, desired );
sqrt( mean(eeg_signal(500:end).^2) );
sqrt( mean(error(500:end).^2) );
plot(1:length(reference), error, 'r'); hold on; plot(eeg_signal, 'b');
grid on; xlim([1 numel(error)]);
title('Filtered Noisy EEG Signal and Clean EEG Signal');
xlabel('Index'); ylabel('Amplitude (unitless)'); legend('Filtered Noisy EEG Signal', 'Clean EEG Signal');
h = tapWeights(:, end);
plot(tapWeights'); grid on; xlim([1 numel(error)]);
title('Filter Tap-weight Convergence Curves');
xlabel('Index'); ylabel('Tap-weights (unitless)'); legend(['Tap 1 -> ' num2str(h(1))], ['Tap 2 -> ' num2str(h(2))]);
figure; spectrogram(eeg_signal, 128, 100, 256, fs); caxis([-122 -21]); colorbar; title('EEG Signal (desired)');
figure; spectrogram(eeg_signal_contaminated, 128, 100, 256, fs); caxis([-122 -21]); colorbar; title('Noise Contaminated EEG Signal (reference)');
figure; spectrogram(error, 128, 64, 256, fs); caxis([-122 -21]); colorbar; title('Recovered (target)');
function [tapWeights, error] = LMS_MRW(stepSize, tapNumber, reference, desired)
tapWeights = zeros(tapNumber, length(desired));
inputVector = zeros(tapNumber, 1);
for updateIndex = 1:1:length(reference)
inputVector(1) = reference(updateIndex);
filterOutput = inputVector'*tapWeights(:, updateIndex);
error(updateIndex) = desired(updateIndex) - filterOutput;
tapWeights(:, updateIndex + 1) = tapWeights(:, updateIndex) + stepSize*inputVector*error(updateIndex);
for inputVectorIndex = length(inputVector):-1:2
inputVector(inputVectorIndex) = inputVector(inputVectorIndex - 1);
The adaptive algorithm used is the LMS with a step-size of 0.0005. The step-size is noticeably smaller and perhaps an extreme value compared to the value of 0.01 stated in the example. However, by using such a small step-size, we can readily see the convergence of the tap-weights and the gradual reduction of the 60 Hz noise in the discrete-time domain.
The reduction of the 60 Hz noise can also be seen by examining the signals using spectrograms. Here are three spectrograms for the EEG signal prior to degradation by the 60 Hz noise, the contaminated EEG signal, and finally the recovered EEG signal with the 60 Hz noise removed by the adaptive filter.
In the third figure, you can readily see the 60 Hz noise being gradually removed by the adaptive filter.
Finally, if you followed the computation of the autocorrelation matrix, R, and the cross-correlation vector, h, you could compute the tap-weight vector, x. The simulation gives values x(1) = -0.0571 and x(2) = - 0.1158, which are close to the theoretical values of -0.0578 and -0.1156, respectively.
If you set the "desired" variable to "noise_added" and "reference" to "noise_signal", the EEG signal is not included, you can remove the 60 Hz tone entirely, as shown by the red curve in the time-domain plot. See the plot below.
Again you can see the 60 Hz noise being attenuated. The step-size was kept at 0.0005.
When the EEG signal is present, however, the 60 Hz noise is reduced, but not entirely eliminated. The third spectrogram shown above is based on a step-size of 0.0005, so the error due to the 60 Hz noise appears to be eliminated after the tap-weights of the filter have converged. However, if the step-size is made larger, say to 0.05, the error is readily seen in the spectrogram.