%Pseudo-Frequencies for CWT
frequencies = [4 2 1 ];

%Properties
wname = 'morl';

%Input
fsample = 500;
t=0:1/fsample:2*pi;

y1 = sin(2*pi*t);
y10 = sin(2*pi*t*2);
y20 = sin(2*pi*t*4);

y = y1+y10+y20;

%scales based on frequencies
scales = freq2scale(frequencies, wname, 1/fsample);

%FFT based CWT transform
y_cwtft = cwtft(y, 'plot', 'wavelet', wname, 'scales', scales );

%Obtain signal based on CWT
y_icwtft = icwtft(y_cwtft,'signal',y,'plot');

%Normalize output between -1 and 1
y_icwtft= (y_icwtft - min(y_icwtft)) / ( max(y_icwtft) - min(y_icwtft) );

%Bring back to input range
y_icwtft= min(y) + (max(y)-min(y))*y_icwtft;


figure
plot(t,y,'r');
legend('Input');
xlabel('Time');
ylabel('Value');
title('Input signal');



figure
plot(t,y,'r');
hold on;
plot(t,y_icwtft,'b');
hold on;
legend('Input','Output');
xlabel('Time');
ylabel('Value');
title('Signal transformation');

figure
plot(t,y-y_icwtft,'b');
hold on;
plot(t,abs(y-y_icwtft),'r');
hold on;
legend('Diff', 'Normalized Diff');
xlabel('Time');
ylabel('Value');
title({'Signal diff',sprintf('summary diff (%0.2f)',sum(abs(y-y_icwtft)))});