Notes

signal = importdata("C:UsersS259-17DesktopabLORENZ.txt");
%signal = normalize(dane1(2000:5000));

%signal = (dane1(260000:263000))';
%signal = (dane1(230000:233000))';
%signal = (dane1(23000:25000))';

Fs = 1000; % czestotliwość próbkowania
t = 0:(1/Fs):(1 - 1/Fs);
comp1 = cos(2*pi*200*t).*(t>0.7);
comp2 = cos(2*pi*60*t).*(t>=0.1 & t<0.3);
trend = sin(2*pi*1/2*t);
rng default
wgnNoise = 0.4*randn(size(t));
signal = comp1 + comp2 + trend + wgnNoise;
signal = normalize(signal);

%signal = cos(2*pi*200*t);

%signal = cos(2*pi*5*t)+0.2*randn(size(t));

% for i = 1:1:1000
% fff = 0;
% if (mod(i, 5) == 0 )
% fff = randn(1,1);
% end
% signal(i) = cos(2*pi*5*t(i)) + fff;
% end
% signal = normalize(signal);

pom_co = 1/Fs; % czas pomiedzy probkami
%Fs = 1/pom_co; % czestotliwość próbkowania
N = length(signal);
ax = 1:N; ax = pom_co * ax;
f = (1:N) / (N * pom_co);

poz = 1;
set(gcf,'position',[10,10,600,1200]);

% filtr _____________________________________________
%x2 = movmedian(signal,10); signal = x2;
%x2 = movmean(signal,10); signal = x2;
% ___________________________________________________

subplot(6,1,1); plot(ax, signal);xlabel('Seconds');
subplot(6,1,2); histogram(signal);

% autokorelacja _____________________________________
[au la] = autocorr(signal,200);
subplot(6,1,3); plot(la, au);xlabel('Lag');

% fft _______________________________________________
ff = fft(signal);
subplot(6,1,4); plot(abs(ff)/N);
zoom = 1; % powiekszenie
subplot(6,1,5); plot(f(1:int16(N/(2*zoom))), 2*abs(ff(1:int16(N/(2*zoom)))/N));xlabel('Hz');
fflenght = length(ff);
r1 = 0.1 * fflenght;
ff(1:r1)= 0;
ff(fflenght - r1:fflenght)= 0;
ff1 = real(ifft(ff));
subplot(6,1,6); plot(ax, ff1);



% atraktor ____________________________________________

dim = 3;
figure;
subplot(2,1,1);
[~,lag] = phaseSpaceReconstruction(signal,[],dim);
lag = 33;
eps = (0.1*((max(signal)) - (min(signal))))^2;
[t, wyk, RQA] = recurent(signal, lag, dim, eps);
RQA_M(poz,:) = RQA(:);
poz = poz + 1;
eRange = 100;
lap = lyapunovExponent(signal,Fs,lag,dim,'ExpansionRange',eRange);
plot3(t(:,1),t(:,2),t(:,3));title("Dane lag = "+lag+" Lap = "+lap);
for ky = 0:500
CCC(ky+1) = nowa_e_d(wyk, 2, ky);
end

%figure;
subplot(2,1,2);
spy(wyk);view([-90 90]);
title("Dane"+" RR= "+ RQA(1) + " DET= "+ RQA(2) + " Lmax= "+ RQA(3) + " ENT= " + RQA(5)+ ...
" LAM= " + RQA(6)+ " ENT1= " + RQA(7)+ " VMAX= " + RQA(8)+ " ent_n= " + RQA(9));
% 1 2 3 4 5 6 7 8 9
%RR DET Lmax DIV ENT LAM ENT1 Vmax nowa_e(wyk, 1000, 2);
set(gcf,'position',[1800,10,600,1200]);
for kk = 2:lev+1
figure;
subplot(2,1,1);
[~,lag] = phaseSpaceReconstruction(mra(kk-1,200:N),[],dim);
ddane = normalize(mra(kk-1,200:N));
eps = (0.02*(max(ddane) - min(ddane)))^2;
[t, wyk, RQA] = recurent(ddane, lag, dim, eps/(4*kk));
RQA_M(poz,:)=RQA(:);
poz = poz + 1;
eRange = 100;
lap = lyapunovExponent(ddane,Fs,lag,dim,'ExpansionRange',eRange);
plot3(t(:,1),t(:,2),t(:,3));title("D" + (kk-1) + " lag = " + lag + " Lap = " + lap);
%figure;
subplot(2,1,2);
spy(wyk);view([-90 90]);title("D" + (kk-1)+" RR= "+ RQA(1) + " DET= "+ RQA(2) + ...
" Lmax= "+ RQA(3) + " ENT= " + RQA(5)+ " LAM= " + RQA(6)+ " ENT1= " + RQA(7)+ ...
" VMAX= " + RQA(8)+ " ent_n= " + RQA(9));
set(gcf,'position',[1800,10,600,1200]);
end
figure;
subplot(2,1,1);
[~,lag] = phaseSpaceReconstruction(mra(kk-1,:),[],dim);
eps = (0.02*((max(signal)) - (min(signal))))^2;
[t, wyk, RQA] = recurent(mra(lev+2-1,:), lag, dim,eps);
RQA_M(poz,:) = RQA(:);
poz = poz + 1;
eRange = 100;
lap = lyapunovExponent(mra(kk-1,:),Fs,lag,dim,'ExpansionRange',eRange);
plot3(t(:,1),t(:,2),t(:,3));title("A"+(kk-1)+ " lag = " + lag + " Lap = " + lap);
%figure;
subplot(2,1,2);
spy(wyk);view([-90 90]);title("A"+(kk-1)+" RR= "+ RQA(1) + " DET= "+ RQA(2) + ...
" Lmax= "+ RQA(3) + " ENT= " + RQA(5)+ " LAM= " + RQA(6)+ " ENT1= " + RQA(7)+ ...
" VMAX= " + RQA(8)+ " ent_n= " + RQA(9));
set(gcf,'position',[1800,10,600,1200]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t, wyk, RQA] = recurent(signal, tau, dim, eps)
N = length(signal);
t = zeros(N - tau * dim, dim);
for j = 1 : dim
for i = 1 : (N - tau * dim)
t(i,j) = signal( i + (j - 1) * tau);
end
end
wyk = zeros(N - tau * dim, N - tau * dim);
% wykres rekurencyjny _______________________________
for i = 1:(N-tau*dim)
x = t(i,:);
for j = 1:(N-tau*dim)
x1 = t(j,:);
d = (x-x1).^2;
s = sum(d,"all");
if s < eps
wyk(i,j)=1;
end
end
end
RR = sum(wyk,"All")/((N-tau*dim)*(N-tau*dim));

L = zeros(length(wyk(:,1)),length(wyk(:,1)));
for i = 1:(N-tau*dim)
for j = 1:((N-tau*dim)-i)
if wyk(i+j-1,j) == 1
L(i,j)=1;
end
end
end

dL = zeros(1,length(wyk(:,1)));
for i = 1:length(wyk(:,1))
count = 0;
for j = 1 : length(wyk(1,:)) - i
if L(i,j) == 1
count = count + 1;
end
if L(i,j) == 0
if count > 0
dL(count) = dL(count) + 1;
count = 0;
end
end
end
end
dL = 2* dL;
sdl = 0;
for i = 2:length(dL)
sdl = sdl + dL(i)*i;
end
DET = sdl/(RR*((N-tau*dim)*(N-tau*dim)));
LS = dL(1);
Lmax = 1;
for i = length(dL):-1:1
if dL(i) > 0
Lmax = i;
break;
end
end
DIV = 1/Lmax;
pp(1) = 0;
for i = 1:length(dL)
pp(i) = dL(i)*i;
end
pp = pp./sdl;
ENT = 0;
for j = 2:length(wyk(1,:))
if pp(j) > 0
ENT = ENT + pp(j)*log(pp(j));
end
end
ENT = -ENT;

RQA(1) = RR;
RQA(2) = DET;
RQA(3) = Lmax;
RQA(4) = DIV;
RQA(5) = ENT;

vL = zeros(1,length(wyk(:,1)));
count = 0;
for i = 1:length(wyk(:,1))
count = 0;
for j = (i+1):length(wyk(:,1))
if wyk(i,j) == 1
count = count + 1;
end
if ((wyk(i,j) == 0) | (j == length(wyk(:,1))))
if count > 0
vL(count) = vL(count) + 1;
count = 0;
end
end
end
end
sdl = 0;
for i = 2:length(vL)
sdl = sdl + vL(i)*i;
end
LAM = sdl/(RR*((N-tau*dim)*(N-tau*dim)));
pp1(1) = 0;
for i = 1:length(vL)
pp1(i) = vL(i)*i;
end
pp1 = pp1./sdl;
ENT1 = 0;
for j = 2:length(wyk(1,:))
if pp1(j) > 0
ENT1 = ENT1 + pp1(j)*log(pp1(j));
end
end
ENT1 = -ENT1;
RQA(6) = LAM;
RQA(7) = ENT1;
LS = vL(1);
Vmax = 1;
for i = length(vL):-1:1
if vL(i) > 0
Vmax = i;
break;
end
end
RQA(8) = Vmax;
RQA(9) = nowa_e(wyk, 1000, 2);
RQA(10) = nowa_e_d(wyk, 2, 0);
end

function ent = nowa_e(R, MM, N)
%R = crp(dane,dim,tau,eps,'euclidean', 'silent');
[a b] = size( R );
cc = zeros(1, 2^(N^2));
for i = 1:MM
x = randi( a - N );
y = randi( b - N );
n = R(x : x + N - 1, y : y + N - 1);

nn = n( 1, : );
for j = 2 : N
nn = [nn n(j,:)];
end;
k1 = length(nn);
k = bit2int( nn', k1 );
cc( k + 1 ) = cc( k + 1 ) + 1;
end;
su = sum( cc );
ent = 0;
for i = 1 : 2^(N^2)
p = cc( i ) / su;
if p ~= 0
ent = ent + p * log( p );
end;
end;
ent = -ent;
end

function ent = nowa_e_d(R, N, tau)
%R = crp(dane,dim,tau,eps,'euclidean', 'silent');
[a b] = size( R );
cc = zeros(1, 2^(N^2));
for i = 1:(a - N - tau)
x = i; y = i + tau;
n = R(x : x + N - 1, y : y + N - 1);
nn = n( 1, : );
for j = 2 : N
nn = [nn n(j,:)];
end;
k1 = length(nn);
k = bit2int( nn', k1 );
cc( k + 1 ) = cc( k + 1 ) + 1;
end;
su = sum( cc );
ent = 0;
for i = 1 : 2^(N^2)
p = cc( i ) / su;
if p ~= 0
ent = ent + p * log( p );
end;
end;
ent = -ent;
end