Notes

clear;
close all;

Fs = 1e3; % czestotliwosc probkowania
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;
Noise = 0.4*randn(size(t));

signal = comp1 + comp2 + trend + Noise;
plot(signal);

pom_co = 0.001; % czas pomiedzy probkami
Fs = 1/pom_co; % czestotliwosc probkowania
t = 0:1/Fs:1-1/Fs;
N = length(signal);
ax = 1 : N; ax = ax * pom_co;
f = (1:N)/(N*pom_co);

set(gcf, 'position',[20,20,400,800])
subplot(6,1,1); plot(ax,signal); xlabel('Seconds');

subplot(6,1,2); histogram(signal);

signal = normalize(signal);

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

%fft
ff = fft(signal);
zoom = 15;
subplot(6,1,4); plot(f(1:int16(N/(2*zoom))), 2*abs(ff(1:int16(N/(2*zoom))))); xlabel("Hz");

%falki
sig = struct('val', signal, 'wavelet', {'morl1'}, 'period', pom_co);
cwtS1 = cwtft(sig);
scales = cwtS1.scales;
MorletFourierFactor = 4*pi / (6 + sqrt (2+6^2));
freq = 1 ./ (scales.*MorletFourierFactor);
subplot(6,1,5); contour(ax, freq, real(cwtS1.cfs));
xlabel('seconds'); ylabel('Hz'); axis([0 ax(end) 0 Fs/(2*zoom)]);

%fractal
[dhseoul, hseoul, cpseoul] = dwtleader(signal);
subplot(6,1,6); plot(hseoul, dhseoul, 'g--*');
xlabel('h'); ylabel('D(h)'); grid on;

%dekompozycja falkowa
lev = 6; %poziom dekompozycji
mra = modwtmra(modwt(signal, "db6", lev));
figure;
set(gcf, 'position', [650,50,500,900]);
subplot(lev+2,1,1);plot(ax,signal); title("Dane");
subplot(lev+2,1,2);plot(ax(200:N),mra(1,200:N)); title("D1")
for kk = 2:lev+1
subplot(lev+2,1,kk); plot(ax(200:N),mra(kk-1,200:N)); title("D"+(kk-1));
end
subplot(lev+2,1,kk+1);
plot(ax, mra(lev+1,:)); title("A"+(kk-1));
xlabel("Time (s)");


figure;
set(gcf, 'position', [1300,50,500,900]);
subplot(lev+2,1,1);plot(ax,signal); title("Dane");
for kk = 2:lev+1
subplot(lev+2,1,kk);
ff = fft(mra(kk-1,:));
zoom = 1*((kk-1));
plot(f(1:int16(N/(2*zoom))), 2*abs(ff(1:int16(N/(2*zoom))))/N);
xlabel("Hz"); title("D"+(kk-1));
end
subplot(lev+2,1,lev+2);
ff = fft(mra(lev+1,:));
zoom = 1*((lev+1));
plot(f(1:int16(N/(2*zoom))), 2*abs(ff(1:int16(N/(2*zoom))))/N);
xlabel("Hz"); title("A"+(kk-1));

%moc dekompozycji falkowej
[c,l]=wavedec(signal,lev,"db6");
[Ea,Ed] = wenergy(c,l);
bar([Ed,Ea]);

%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));
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]);

% atraktor2 ____________________________________________

dim = 3;
figure;
subplot(2,1,1);
[~,lag] = phaseSpaceReconstruction(signal,[],dim);
lag = 33;
[t, wyk1, RQA1] = recurent2(signal, lag, dim);
RQA_M(poz,:) = RQA1(:);
poz = poz + 1;
eRange = 100;
%lap = lyapunovExponent(signal,Fs,lag,dim,'ExpansionRange',eRange);
lap = 0;
plot3(t(:,1),t(:,2),t(:,3));title("Dane lag = "+lag+" Lap = "+lap);
%figure;
subplot(2,1,2);
spy(wyk1);view([-90 90]);title("Dane"+" RR= "+ RQA1(1) + " DET= "+ RQA1(2) + " Lmax= "+ RQA1(3) + " ENT= " + RQA1(5)+ " LAM= " + RQA1(6));
%txt = {'Plotted Data:','y = sin(x)'};
%text(400,500,txt);
set(gcf,'position',[1800,10,600,1200]);
for ky = 0:500
CCC1(ky+1) = nowa_e_d(wyk1, 2, ky);
end

% entropia ____________________________________________

figure
CMSE = ent_w(signal);
plot(CMSE);
CI = trapz(CMSE); %pole pod wykresem CMSE liczone metodą trapezów

% Falki jeszcze raz ____________________________________________

figure
zoom = 15;
contour(ax, freq, real(cwtS1.cfs));
xlabel('Seconds'); ylabel('Hz'); axis([0 ax(end) 0 Fs/(2*zoom)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CMSE = ent_w(signal)
szereg=signal';
for factor=1:50 %liczba skali czasowych
% tak dobrać factor, aby size(y,2)/factor >=300
p=1; x=1;
while p<(size(szereg,1)-factor);
A(x, factor)=sum(szereg(p:(p+factor-1)))/factor;
p=p+1; %k+factor, wtedy MSE; k=k+1 dla CMSE
x=x+1;
B(factor)=x-1; %tablica tymczasowa przechowująca długości uśrednionych szeregów
end
C(factor)=x-1; %C przechowaj długości uśrednionych szeregów
%SD(j - wiersz - nr ramki; factor - kolumna - poziom uśrednienia)
SD(factor)=std(A(1:(C(factor)), factor));
r(factor)=0.10*SD(factor); %0.2*SD(factor);
end
for factor=1:50
nd=0; nn=0;
for i=1:(C(factor)-2)
for p=(i+1):(C(factor)-2)
if (abs(A(i,factor)-A(p,factor)) < r(factor)) & (abs(A(i+1,factor)-A(p+1,factor))< r(factor))
nn = nn + 1;
if (abs(A(i+2,factor)-A(p+2,factor))< r(factor))
nd = nd + 1;
end
end
end
end
CMSE(factor)=-log2(nd/nn);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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