time = (0:0.01:6)';
x = sin(40./(time+0.01));

figure,plot(time, x)
title('Information Signal x')
xlabel('time')
ylabel('x')


n1 = randn(size(time));
figure,plot(time, n1)
title('Noise Source n1')
xlabel('time')
ylabel('n1')


domain = linspace(min(n1), max(n1), 20);
[xx, yy] = meshgrid(domain, domain);
zz = 4*sin(xx).*yy./(1+yy.^2);

figure,surf(xx, yy, zz);
xlabel('n1(k)');
ylabel('n1(k-1)');
zlabel('n2(k)');
title('Unknown Channel Characteristics That Generate Interference');
axis([min(domain) max(domain) min(domain)  max(domain) min(zz(:)) max(zz(:))]);
set(gca, 'box', 'on');

n1d0 = n1;                            % n1 delay 0
n1d1 = [0; n1d0(1:length(n1d0)-1)];   % n1 delay 1
n2 = 4*sin(n1d0).*n1d1./(1+n1d1.^2);  % interference
axis_limit = [min(time) max(time) min([n1;n2]) max([n1;n2])];
figure,subplot(2,1,1)
plot(time, n1);
ylabel('noise n1'); axis(axis_limit);
subplot(2,1,2)
plot(time, n2);
ylabel('interference n2'); axis(axis_limit);

m = x + n2;             % measured signal
figure,plot(time, m)
title('Measured Signal')
xlabel('time')
ylabel('m')

delayed_n1 = [0; n1(1:length(n1)-1)];
trn_data = [delayed_n1 n1 m];
mf_n = 2;
ss = 0.2;
in_fismat=genfis1(trn_data, mf_n);
out_fismat = anfis(trn_data, in_fismat, [nan nan ss]);
estimated_n2 = evalfis(trn_data(:, 1:2), out_fismat);
estimated_x = m - estimated_n2;

axis_limit =  [min(time) max(time) min([n2; estimated_n2]) max([n2; estimated_n2])];
figure,subplot(2,1,1)
plot(time, n2)
ylabel('n2 (unknown)'); axis(axis_limit);

subplot(2,1,2)
plot(time, estimated_n2)
ylabel('estimated_n2'); axis(axis_limit);

axis_limit =  [min(time) max(time) min([x; estimated_x]) max([x; estimated_x])];
subplot(2,1,1)
plot(time, x)
ylabel('x (unknown) '); axis(axis_limit);
subplot(2,1,2)
plot(time, estimated_x)
axis(axis_limit); ylabel('estimated\_x')