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

%Properties
wname = 'morl';
characteristic_frequency = 4;
freq_range = 1;

%Effect

%phase_delay_factor = 0;
phase_delay_factor = 0.125;
frequency_factor = 1;

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

y1 = 0;%sin(2*pi*t);
y10 = 0;% 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);

%Frequencies corresponding to scales
frequencies = scal2frq(y_cwtft.scales,y_cwtft.wav,1/fsample);

%Filter activation for scales / frequencies
activation_values = 1- min(1,abs(frequencies - characteristic_frequency)/freq_range);

%Filter frequency dump values for scales / frequencies
dump_values =  frequency_factor + (1-frequency_factor)*(1-activation_values);

%Filter frequency delay values for scales / frequencies
delay_values = floor(fsample * phase_delay_factor * activation_values);

for i = 1:size(y_cwtft.scales,2)   
    %Modify CWT for given scale
    y_cwtft.cfs(i,:) = dump_values(i) * circshift(y_cwtft.cfs(i,:)',delay_values(i))';
end
 
%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(frequencies,dump_values,'r');
hold on;
plot(frequencies,activation_values,'b');
hold on;
legend('Frequency cutoff','Filter activation');
xlabel('Frequency (Hz)');
ylabel('Value');
title('Filter activation');

figure
plot(frequencies,delay_values,'g');
hold on;
legend('Delay (frames)');
xlabel('Frequency (Hz)');
ylabel('Value');
title('Delay');

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;
legend('Diff');
xlabel('Time');
ylabel('Value');
title('Signal diff');