diff --git a/SPARE.m b/SPARE.m
new file mode 100644
index 0000000..4b5e348
--- /dev/null
+++ b/SPARE.m
@@ -0,0 +1,301 @@
+function [new_noNorm, old, gt] = SPARE(ground_truth, raw_PPG, time_PPG, winLen, stride, fs_PPG, True_BPMs, plotting, log)
+
+%len = round(winLen*0.3086*fs_PPG);
+
+len = stride * fs_PPG;
+% Initialize
+old = zeros(1, len);
+gt = zeros(1, len);
+new_noNorm = zeros(1, len);
+times = time_PPG(1:len);
+
+count = 0;
+% Filter the data
+[dataFull, ~] = PPG_filter(raw_PPG, fs_PPG);
+
+signal_end = time_PPG(end);
+segment_start = time_PPG(1);
+segment_end = time_PPG(1) + winLen;
+
+
+while(segment_end < signal_end)
+    count = count+1;
+    %disp(count)
+    if mod(count,1) == 0
+        disp((1-(signal_end-segment_end)/(signal_end-time_PPG(1)))*100)
+    end
+    data = getSignalSegmentTEST(dataFull, time_PPG, segment_start, segment_end);
+    data_t = getSignalSegmentTEST(time_PPG, time_PPG, segment_start, segment_end);
+    data_ref = getSignalSegmentTEST(ground_truth, time_PPG, segment_start, segment_end);
+    data_noise = getSignalSegmentTEST(raw_PPG, time_PPG, segment_start, segment_end);
+
+    %% 
+    
+    if (~all(isnan(data)))
+        RR_mean = 1/True_BPMs(count)*60;
+        if (log)
+            disp('BPM from ECG:')
+            disp((1/RR_mean*60))
+        end
+        
+        PPG = data';
+        
+        % remove NaNs from signal
+        PPG(isnan(PPG)) = 0;
+
+        % devide the PPG in components
+        comps = SSA(PPG, 10, floor(length(PPG)*0.75/2) -1);
+        N = 60*fs_PPG;
+        b = zeros(N,1);
+        
+        % define frequency axis
+        %f = (0:fs_PPG/(N-1):fs_PPG)*60; 
+        
+        % do FFT per components
+        for i=1:length(comps(:,1))
+            [BPM, s] = SSR(comps(i,:)', fs_PPG, N);
+            b = b + s;
+        end
+        
+        b = b(40:500);
+        BPM = BPM(40:500);
+        
+        p = abs(b).^2;
+        BPMs = BPM';
+        
+        % find the index of all peaks
+        [~,index_P] = findpeaks(p);
+        index_Full = index_P;
+
+        % BPM of all peaks
+        BPM_peaks_full = BPMs(index_Full);
+
+        % compute second derivative
+        second = diff(diff(comps'))';
+        
+        % FFT per component of the second derivative 
+        a = zeros(N,1);
+        for i=1:length(second(:,1))
+            [~, Comps] = SSR(second(i,:)', fs_PPG, N);
+            a = a + Comps;
+        end
+        
+        a = a(40:500);
+        a = abs(a).^2;
+       
+        % find all peaks of the spectrum of the second derivative
+        [pks,index_P] = findpeaks(a);
+        
+        % BPM of all the peaks
+        BPM_peaks = BPMs(index_P);
+        
+        % compute distance between the biggests peaks and the HR from ECG
+        distance = abs(BPM_peaks-((1/RR_mean)*60));
+                
+        % remove the index that are less than 10% of the biggest one among
+        % the candidates
+        index_cand = find(distance<7);
+        pks_cands = pks(index_cand);
+        maximum = max(pks(index_cand));
+        
+        index_P_first = index_cand(pks_cands>=0.2*maximum);
+        
+        % candidates for being the representative of the fundamental of the signal
+        BPM_cand = BPM_peaks(index_P_first);
+        
+        % same but for the second harmonic
+        distance = abs(BPM_peaks-((1/RR_mean)*60*2));
+        BPM_cand2 = BPM_peaks(distance<14);
+        
+        % same but for the third harmonic
+        distance = abs(BPM_peaks-((1/RR_mean)*60*3));
+        BPM_cand3 = BPM_peaks(distance<20);
+        
+        % from the three sets of peaks (BPM_cand, BPM_cand2, BPM_cand3),
+        % find the triplet that staisfy the condition:
+        % (HR_cand in BPM_cand, HR_cand2 in BPM_cand2, HR_cand3 in BPM_cand3) s.t. HR_cand2/HR_cand = 2, HR_cand3/HR_cand = 3 
+        valMin = +inf;
+        valMin2 = +inf;
+        for i=1:numel(BPM_cand)
+            for j=1:numel(BPM_cand2)
+                for k=1:numel(BPM_cand3)
+                    val = abs(BPM_cand2(j)/BPM_cand(i)-2);
+                    val2 = abs(BPM_cand3(k)/BPM_cand(i)-3);
+                    if (val <= valMin && val2 <= valMin2)
+                        valMin = val;
+                        valMin2 = val2;
+                        HR_cand = BPM_cand(i);
+                        HR_cand2 = BPM_cand2(j);
+                        HR_cand3 = BPM_cand3(k);
+                    end
+                end
+            end
+        end
+        
+        % sanity check: HR_cand/HR_cand is at most 30% away from 2 and
+        % HR_cand do not differ much from HR from ECG
+        if (valMin<0.3 && abs(HR_cand-(1/RR_mean)*60) < 15)
+
+            HR_p = HR_cand;
+            HR_p2 = HR_cand2;
+            HR_p3 = HR_cand3;
+        else
+            % otherwise, use HR from ECG (extrema ratio)
+            HR_p = (1/RR_mean)*60;
+            HR_p2 = (1/RR_mean)*60*2;
+            HR_p3 = (1/RR_mean)*60*3;
+        end
+
+        % find index of the peaks computed from the spectrum of the orignal
+        % signal
+        [~, idx] = min(abs(BPM_peaks_full-(HR_p)));
+        [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx))));
+
+        if (log)
+            disp('HR peak from PPG:')
+
+            disp(BPMs(idx))
+        end
+        HR_2 = BPMs(idx);
+        BPM_fund = BPMs(idx);
+        p_fund = p(idx);
+
+        max_idx = idx;
+        while (max_idx<numel(p) && p(max_idx+1)<=p(max_idx)) 
+            max_idx = max_idx+1;
+        end
+        min_idx = idx;
+        while (min_idx-1>0 && p(min_idx-1)<=p(min_idx))
+            min_idx = min_idx-1;
+        end
+
+
+        f_min = BPMs(min_idx)/60;
+        f_max = BPMs(max_idx)/60;
+        distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean))*60 );
+        if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>10)
+            f_min = ((1/RR_mean)*60 -15)/60;
+            f_max = ((1/RR_mean)*60 +15)/60;
+
+        end
+        x = bandpass(PPG, [f_min f_max], fs_PPG);
+        HR = ((f_max+f_min)/2)*60;
+
+        [~, idx] = min(abs(BPM_peaks_full-(HR_p2)));
+        [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx))));
+        if (log)
+            disp('2th HR peak from PPG:')
+
+            disp(BPMs(idx))
+        end
+        BPM_second = BPMs(idx);
+        p_second = p(idx);
+
+
+
+        max_idx = idx;
+        while (max_idx<numel(p) && p(max_idx+1)<=p(max_idx) )
+            max_idx = max_idx+1;
+        end
+        min_idx = idx;
+        while (min_idx-1>0 && p(min_idx-1)<=p(min_idx))
+            min_idx = min_idx-1;
+        end
+
+        f_min = BPMs(min_idx)/60;
+        f_max = BPMs(max_idx)/60;
+        distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean*2))*60 );
+        if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>15)
+            f_min = ((1/RR_mean)*60*2 -15)/60;
+            f_max = ((1/RR_mean)*60*2 +15)/60;
+        end
+        y = bandpass(PPG, [f_min f_max], fs_PPG);
+
+
+        [~, idx] = min(abs(BPM_peaks_full-(HR_p3)));
+        [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx))));
+
+        BPM_third = BPMs(idx);
+        p_third = p(idx);
+
+        max_idx = idx;
+        while (max_idx<numel(p) && p(max_idx+1)<=p(max_idx))
+            max_idx = max_idx+1;
+        end
+        min_idx = idx;
+        while (min_idx-1>0 && p(min_idx-1)<=p(min_idx))
+            min_idx = min_idx-1;
+        end
+
+        f_min = BPMs(min_idx)/60;
+        f_max = BPMs(max_idx)/60;
+        distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean*3))*60 );
+
+        if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>25)
+            f_max = 1/RR_mean*3+30/60;
+            f_min = 1/RR_mean*3-30/60;
+        end
+        z = bandpass(PPG, [f_min f_max], fs_PPG);
+
+
+        segment_start = segment_start+len/fs_PPG;
+        segment_end = segment_start + winLen;
+
+        Hal = x+y+z;
+        %disp(kurtosis(raw_PPG))
+%         if (kurtosis(raw_PPG)<2.5)
+%             Hal = data_noise.';
+%         end
+
+        if (len*2 <= numel(Hal))
+            if (max(abs(data_noise(len:len*2))) > 2.5)
+                new_noNorm = [new_noNorm, Hal(len:len*2-1)];
+            else
+                new_noNorm = [new_noNorm, data_noise(len:len*2-1)];
+            end
+            old = [old, data_noise(len:len*2-1).'];
+            gt = [gt, data_ref(len:len*2-1)];
+            times = [times, data_t(len:len*2-1)];
+        else
+            new_noNorm = [new_noNorm, Hal(len:end)];
+            old = [old, data_noise(len:end).'];
+            gt = [gt, data_ref(len:end)];
+            times = [times, data_t(len:end)];
+        end
+
+        if (plotting)
+            try
+            figure
+            plot(subplot(2,1,1), Hal(len:len*2-1));    
+            hold on; plot(subplot(2,1,1), PPG(len:len*2-1));
+
+            legend('filtered', 'original')
+            %[~, p_ref] = mySpectrum(data_ref, fs_PPG);
+            plot(subplot(2,1,2), BPMs, p);
+            hold on;
+            plot(subplot(2,1,2), BPMs, p);
+            plot(subplot(2,1,2), BPM_fund,p_fund,'r*');
+            plot(subplot(2,1,2), BPM_second,p_second,'r*');
+            plot(subplot(2,1,2), BPM_third,p_third,'r*');
+            xlim([40, 350])
+            catch e
+                disp(e)
+            end
+        end
+        oldRR = RR_mean;
+    else
+        aiaia = 0;
+        segment_start = segment_start+len/fs_PPG;
+        segment_end = segment_start + winLen;
+    end
+     
+end
+end
+
+
+function data_segment = getSignalSegmentTEST(signal, signal_t, start, ending)
+
+index = find((start<signal_t) & (signal_t<ending));
+data_segment = signal(index);
+
+end
\ No newline at end of file
diff --git a/fftSPARE.m b/fftSPARE.m
new file mode 100644
index 0000000..b9e7f79
--- /dev/null
+++ b/fftSPARE.m
@@ -0,0 +1,298 @@
+function [new_noNorm, old, gt] = fftSPARE(ground_truth, raw_PPG, time_PPG, winLen, stride, fs_PPG, True_BPMs, plotting, log)
+
+len = stride * fs_PPG;
+% Initialize
+old = zeros(1, len);
+gt = zeros(1, len);
+new_noNorm = zeros(1, len);
+times = time_PPG(1:len);
+
+count = 0;
+% Filter the data
+[dataFull, ~] = PPG_filter(raw_PPG, fs_PPG);
+
+signal_end = time_PPG(end);
+segment_start = time_PPG(1);
+segment_end = time_PPG(1) + winLen;
+
+
+while(segment_end < signal_end)
+    count = count+1;
+    %disp(count)
+    
+    data = getSignalSegmentTEST(dataFull, time_PPG, segment_start, segment_end);
+    data_t = getSignalSegmentTEST(time_PPG, time_PPG, segment_start, segment_end);
+    data_ref = getSignalSegmentTEST(ground_truth, time_PPG, segment_start, segment_end);
+    data_noise = getSignalSegmentTEST(raw_PPG, time_PPG, segment_start, segment_end);
+
+    %% 
+    
+    if (~all(isnan(data)))
+        RR_mean = 1/True_BPMs(count)*60;
+        if (log)
+            disp('BPM from ECG:')
+            disp((1/RR_mean*60))
+        end
+        
+        PPG = data';
+        
+        % remove NaNs from signal
+        PPG(isnan(PPG)) = 0;
+
+        % devide the PPG in components
+        comps = SSA(PPG, 10, floor(length(PPG)*0.75/2) -1);
+        N = 60*fs_PPG;
+        b = zeros(N,1);
+        
+        % define frequency axis
+        f = (0:fs_PPG/(N-1):fs_PPG)*60; 
+        
+        % do FFT per components
+        for i=1:length(comps(:,1))
+            Comps = fft(comps(i,:)', N);
+            b = b + Comps;
+        end
+        
+        % just interested in the amplitude
+        P2 = abs(b);
+        % limit the analysis up to 400 BPM
+        L = 400;
+        P1 = P2(1:L);
+        P1(2:end-1) = 2*P1(2:end-1);
+        p = P1;
+        BPMs = f(1:L)';
+        
+        % find the index of all peaks
+        [~,index_P] = findpeaks(p);
+        index_Full = index_P;
+
+        % BPM of all peaks
+        BPM_peaks_full = BPMs(index_Full);
+
+        % compute second derivative
+        second = diff(diff(comps'))';
+        
+        % FFT per component of the second derivative 
+        a = zeros(N,1);
+        for i=1:length(second(:,1))
+            Comps = fft(second(i,:)', N);
+            a = a + Comps;
+        end
+        P2 = abs(a);
+        P1 = P2(1:L);
+        P1(2:end-1) = 2*P1(2:end-1);
+        
+        % find all peaks of the spectrum of the second derivative
+        a = P1;
+        % find all peaks of the spectrum of the second derivative
+        [pks,index_P] = findpeaks(a);
+        
+        % BPM of all the peaks
+        BPM_peaks = BPMs(index_P);
+        
+        % compute distance between the biggests peaks and the HR from ECG
+        distance = abs(BPM_peaks-((1/RR_mean)*60));
+                
+        % remove the index that are less than 10% of the biggest one among
+        % the candidates
+        index_cand = find(distance<7);
+        pks_cands = pks(index_cand);
+        maximum = max(pks(index_cand));
+        
+        index_P_first = index_cand(pks_cands>=0.2*maximum);
+        
+        % candidates for being the representative of the fundamental of the signal
+        BPM_cand = BPM_peaks(index_P_first);
+        
+        % same but for the second harmonic
+        distance = abs(BPM_peaks-((1/RR_mean)*60*2));
+        BPM_cand2 = BPM_peaks(distance<14);
+        
+        % same but for the third harmonic
+        distance = abs(BPM_peaks-((1/RR_mean)*60*3));
+        BPM_cand3 = BPM_peaks(distance<20);
+        
+        % from the three sets of peaks (BPM_cand, BPM_cand2, BPM_cand3),
+        % find the triplet that staisfy the condition:
+        % (HR_cand in BPM_cand, HR_cand2 in BPM_cand2, HR_cand3 in BPM_cand3) s.t. HR_cand2/HR_cand = 2, HR_cand3/HR_cand = 3 
+        valMin = +inf;
+        valMin2 = +inf;
+        for i=1:numel(BPM_cand)
+            for j=1:numel(BPM_cand2)
+                for k=1:numel(BPM_cand3)
+                    val = abs(BPM_cand2(j)/BPM_cand(i)-2);
+                    val2 = abs(BPM_cand3(k)/BPM_cand(i)-3);
+                    if (val <= valMin && val2 <= valMin2)
+                        valMin = val;
+                        valMin2 = val2;
+                        HR_cand = BPM_cand(i);
+                        HR_cand2 = BPM_cand2(j);
+                        HR_cand3 = BPM_cand3(k);
+                    end
+                end
+            end
+        end
+        
+        % sanity check: HR_cand/HR_cand is at most 30% away from 2 and
+        % HR_cand do not differ much from HR from ECG
+        if (valMin<0.3 && abs(HR_cand-(1/RR_mean)*60) < 15)
+
+            HR_p = HR_cand;
+            HR_p2 = HR_cand2;
+            HR_p3 = HR_cand3;
+        else
+            % otherwise, use HR from ECG (extrema ratio)
+            HR_p = (1/RR_mean)*60;
+            HR_p2 = (1/RR_mean)*60*2;
+            HR_p3 = (1/RR_mean)*60*3;
+        end
+
+        % find index of the peaks computed from the spectrum of the orignal
+        % signal
+        [~, idx] = min(abs(BPM_peaks_full-(HR_p)));
+        [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx))));
+
+        if (log)
+            disp('HR peak from PPG:')
+
+            disp(BPMs(idx))
+        end
+        HR_2 = BPMs(idx);
+        BPM_fund = BPMs(idx);
+        p_fund = p(idx);
+
+        max_idx = idx;
+        while (max_idx<numel(p) && p(max_idx+1)<=p(max_idx)) 
+            max_idx = max_idx+1;
+        end
+        min_idx = idx;
+        while (min_idx-1>0 && p(min_idx-1)<=p(min_idx))
+            min_idx = min_idx-1;
+        end
+
+
+        f_min = BPMs(min_idx)/60;
+        f_max = BPMs(max_idx)/60;
+        distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean))*60 );
+        if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>10)
+            f_min = ((1/RR_mean)*60 -15)/60;
+            f_max = ((1/RR_mean)*60 +15)/60;
+
+        end
+        x = bandpass(PPG, [f_min f_max], fs_PPG);
+        HR = ((f_max+f_min)/2)*60;
+
+        [~, idx] = min(abs(BPM_peaks_full-(HR_p2)));
+        [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx))));
+        if (log)
+            disp('2th HR peak from PPG:')
+
+            disp(BPMs(idx))
+        end
+        BPM_second = BPMs(idx);
+        p_second = p(idx);
+
+
+
+        max_idx = idx;
+        while (max_idx<numel(p) && p(max_idx+1)<=p(max_idx) )
+            max_idx = max_idx+1;
+        end
+        min_idx = idx;
+        while (min_idx-1>0 && p(min_idx-1)<=p(min_idx))
+            min_idx = min_idx-1;
+        end
+
+        f_min = BPMs(min_idx)/60;
+        f_max = BPMs(max_idx)/60;
+        distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean*2))*60 );
+        if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>15)
+            f_min = ((1/RR_mean)*60*2 -15)/60;
+            f_max = ((1/RR_mean)*60*2 +15)/60;
+        end
+        y = bandpass(PPG, [f_min f_max], fs_PPG);
+
+
+        [~, idx] = min(abs(BPM_peaks_full-(HR_p3)));
+        [~, idx] = min(abs(BPMs-(BPM_peaks_full(idx))));
+
+        BPM_third = BPMs(idx);
+        p_third = p(idx);
+
+        max_idx = idx;
+        while (max_idx<numel(p) && p(max_idx+1)<=p(max_idx))
+            max_idx = max_idx+1;
+        end
+        min_idx = idx;
+        while (min_idx-1>0 && p(min_idx-1)<=p(min_idx))
+            min_idx = min_idx-1;
+        end
+
+        f_min = BPMs(min_idx)/60;
+        f_max = BPMs(max_idx)/60;
+        distance2HR = abs( ((f_max+f_min)/2 - (1/RR_mean*3))*60 );
+
+        if ((f_max-f_min)<0.02 || (f_max-f_min>1) || distance2HR>25)
+            f_max = 1/RR_mean*3+30/60;
+            f_min = 1/RR_mean*3-30/60;
+        end
+        z = bandpass(PPG, [f_min f_max], fs_PPG);
+
+
+        segment_start = segment_start+len/fs_PPG;
+        segment_end = segment_start + winLen;
+
+        Hal = x+y+z;
+        %disp(kurtosis(raw_PPG))
+        if (kurtosis(raw_PPG)<2.5)
+            Hal = data_noise.';
+        end
+
+        if (len*2 <= numel(Hal))
+            if (max(abs(data_noise(len:len*2))) > 2.5)
+                new_noNorm = [new_noNorm, Hal(len:len*2-1)];
+            else
+                new_noNorm = [new_noNorm, data_noise(len:len*2-1)];
+            end
+            old = [old, data_noise(len:len*2-1).'];
+            gt = [gt, data_ref(len:len*2-1)];
+            times = [times, data_t(len:len*2-1)];
+        else
+            new_noNorm = [new_noNorm, Hal(len:end)];
+            old = [old, data_noise(len:end).'];
+            gt = [gt, data_ref(len:end)];
+            times = [times, data_t(len:end)];
+        end
+
+        if (plotting)
+            figure
+            plot(subplot(2,1,1), Hal(len:len*2-1));    
+            hold on; plot(subplot(2,1,1), PPG(len:len*2-1));
+
+            legend('filtered', 'original')
+            [~, p_ref] = mySpectrum(data_ref, fs_PPG);
+            plot(subplot(2,1,2), BPMs, p);
+            hold on;
+            plot(subplot(2,1,2), BPMs, p_ref);
+            plot(subplot(2,1,2), BPM_fund,p_fund,'r*');
+            plot(subplot(2,1,2), BPM_second,p_second,'r*');
+            plot(subplot(2,1,2), BPM_third,p_third,'r*');
+            xlim([40, 350])
+        end
+        oldRR = RR_mean;
+    else
+        aiaia = 0;
+        segment_start = segment_start+len/fs_PPG;
+        segment_end = segment_start + winLen;
+    end
+     
+end
+end
+
+
+function data_segment = getSignalSegmentTEST(signal, signal_t, start, ending)
+
+index = find((start<signal_t) & (signal_t<ending));
+data_segment = signal(index);
+
+end
\ No newline at end of file
diff --git a/noiseGenerator.m b/noiseGenerator.m
new file mode 100644
index 0000000..a9c148c
--- /dev/null
+++ b/noiseGenerator.m
@@ -0,0 +1,30 @@
+function [noise] = noiseGenerator(kernel_size, execution_fraction, noise_size, sampling_frequency, std, max_freq, plotting)
+
+    stride = kernel_size/execution_fraction;
+    num_actions = floor(noise_size/stride);
+    noise = zeros(1, noise_size);
+    dt = 1/sampling_frequency;
+    time = linspace(-((kernel_size+1)/2-1)*dt, (kernel_size+1)/2*dt, kernel_size);
+    kernel = sinc(time/kernel_size).*sinc(time);
+    
+    if plotting
+        figure
+        plot(kernel)
+    end
+    actions = randn(1,num_actions)*std;
+    actions = lowpass(actions, max_freq, sampling_frequency/stride);
+
+
+    for i=0:num_actions-1
+        if i*stride+kernel_size < numel(noise)
+            noise(i*stride+1:i*stride+kernel_size) = noise(i*stride+1:i*stride+kernel_size)...
+                                               + actions(i+1) * kernel;
+        else 
+            noise(i*stride+1:end) = noise(i*stride+1:end)...
+                                               + actions(i+1) * kernel(1:numel(noise(i*stride+1:end)));
+        end
+    end
+
+    
+
+end
\ No newline at end of file