clear all
close all
clc

%##########################################################################
%##########################################################################
% General information:
% 
% Paper:      Fully Adaptive Savitzky-Golay Type Smoothers, submited in October 2018
% 
% Authors:    Maciej Niedzwiecki, Senior Member, IEEE, and Marcin Ciolek,
%             Member, IEEE
% 
% Conference: IEEE International Conference on Acoustics, Speech, and
%             Signal Processing 2019, Brighton, UK
% 
% Both authors are from
% Gdansk University of Technology 
% Faculty of Electronics, Telecommunications and Informatics 
% Department of Automatic Control
% Poland, 80-233 Gdansk, Narutowicza 11/12 street
% 
% Author of MATALB code: Marcin Ciolek
% E-mail: marcin.ciolek@pg.edu.pl
% Date of code verification: 10.2018
% MATLAB R2017b (9.3.0.713579), License Number: 68140
%##########################################################################
%##########################################################################


%##########################################################################
%##########################################################################
% Input
% 
% y - noisy signal
% so - clean signal
% 
% Output
% 
% FPA_smoothing - the denoised signal obtained based on FPE approach
% CV_smoothing  - the denoised signal obtained based on CV approach
% FPE_SNR - SNR for FPE approach
% CV_SNR  - SNR for CV approach
% SG_SNR  - SNR for 9 nonadaptive Savitzky-Golay type
% smoothers corresponding to different choices of parameters k (20, 40, 80) and n (1, 3, 5)
%##########################################################################
%##########################################################################

%##########################################################################
%##########################################################################

test_signal = 2;

% Chirp_noise_5dB
% Chirp_noise_10dB
% Chirp_noise_15dB
% Chirp_noise_20dB
% Chirp_noise_25dB


% ECG_noise_5dB
% ECG_noise_10dB
% ECG_noise_15dB
% ECG_noise_20dB
% ECG_noise_25dB

if test_signal == 1
    
    load('Input\\Chirp.mat')
    load('Input\\Chirp_noise_10dB.mat')
    
    so = so';
   
else
    
    load('Input\\ECG.mat')
    load('Input\\ECG_noise_10dB.mat')
    
    so    = y';
   
end


%##########################################
%########################################## Simulation settings

N     = length(so);

% the size of local decision window, we recommend m = [20,30] 
m     = 25;

% the maximum order, in parallel scheme algorithms of orders n=1,3,5,..., Nk are considered 
Nk    = 5;
n_tab = 1:2:Nk;

% in parallel scheme algorithms of parameter k: k1,k2,k3,..., k_max are considered  
k_tab = [20,40,80];

n_m   = length(n_tab);
n_k   = length(k_tab);

winodow = 2;
% 1 - rectangular window
% 2 - cosinusoidal window

base = 2;
% 1 - powers of time
% 2 - Fourier series

%##########################################
%##########################################
SNR_all_ave = zeros(n_k,n_m);
SNR_sel_ave = zeros(2,1);

ts       = 142+1;
tk       = N-142;

for num_realization = 1:100

    
disp(num_realization)    
y     = so + u(:,num_realization); 


%##########################################
%##########################################
% Parallel estimation scheme

Z          = zeros(n_k,n_m,N);
FPE_cell   = Z;
CVA_cell   = Z;
Noise_cell = Z;
s_hat_cell = Z;
SNR_all    = zeros(n_k,n_m);

for j = 1:n_k
    
    k  = k_tab(j);  
    
   [FPE,CVA,s_hat,res_noise] = SavitzkyGolayType(y,so,k,m,Nk,base,winodow);

    FPE_cell(j,:,:)   = FPE;
    CVA_cell(j,:,:)   = CVA;
    Noise_cell(j,:,:) = res_noise;
    s_hat_cell(j,:,:) = s_hat;
   
    for i=1:n_m      
           SNR_all(j,i) = snr(so(ts:tk),res_noise(i,ts:tk)');
    end
end

clear FPE CVA res_noise s_hat k j i 
%##########################################
%##########################################



%##########################################
%##########################################
% Parallel smoothing scheme: minimization of the FPE and CVA statistics
% for algorithms with different design parameters k and n


FPE_n      = zeros(N,1);
FPE_k      = zeros(N,1);

CV_n       = zeros(N,1);
CV_k       = zeros(N,1);

FPE_res_noise  = zeros(N,1);
CV_res_noise   = zeros(N,1);

FPE_smoothing  = zeros(N,1);
CV_smoothing   = zeros(N,1);

SNR_sel_all    = zeros(2,1);

for t = ts:tk
    
      FPE   =  FPE_cell(:,:,t);
      CVA   =  CVA_cell(:,:,t);
      Noise =  Noise_cell(:,:,t);
      S_hat =  s_hat_cell(:,:,t);
    
     % FPA based selection
     [FPE_res_noise(t),FPE_smoothing(t),FPE_n(t),FPE_k(t)] = order_bandwidth(FPE,Noise,S_hat,n_k,n_m);
     
     % CV based selection
     [CV_res_noise(t),CV_smoothing(t),CV_n(t),CV_k(t)]     = order_bandwidth(CVA,Noise,S_hat,n_k,n_m);
     
end


SNR_sel_all(1) = snr(so(ts:tk),FPE_res_noise(ts:tk)); 
SNR_sel_all(2) = snr(so(ts:tk),CV_res_noise(ts:tk)); 

SNR_all_ave    = SNR_all_ave + SNR_all;
SNR_sel_ave    = SNR_sel_ave + SNR_sel_all;


 
%##########################################################################
%##########################################################################
%##########################################################################

end



%Averaged SNR
SG_SNR = SNR_all_ave/100;
SNR_sel_ave = SNR_sel_ave/100;

FPE_SNR = SNR_sel_ave(1);
CV_SNR  = SNR_sel_ave(2);

my_figures(y,so,FPE_smoothing,CV_smoothing)

clear FPE CVA Noise S_hat FPE_cell CVA_cell Noise_cell s_hat_cell SNR_sel_all
clear SNR_sel_ave SNR_all_ave Z ts tk t SNR_sel_all CV_res_noise FPE_res_noise

clc
disp('END')

