classdef TapeFeatures
    %TAPE_FEATURES Features of one tape
    %   This class structures all the features of a tape including a vector
    %   of features for every beat of type BeatFeatures
    
    properties
       % beat features
       fm_b;
       
       sum_bf;
       pess_bf;
       perc_an_bf;
       perc_n_bf;
       
       spi;
       hist_spi;
       hist_edges;
       mean_spi;
       mmean_spi;
       mean_spi_win;
       d_mean_spi_win;
       max_hist_ind;
       
       Gt;
       Gtabs;
       Gtsumabs;
       Gtspar;
       Gtmean;
       Gtvar;
       Gtmmax;
       
       wav_feature;
       
       qrs_en_coeff_mean;
       qrs_en_coeff_var;
       qrs_en_coeff_syst_mean;
       qrs_en_coeff_syst_var;
       qrs_en_coeff_diast_mean;
       qrs_en_coeff_diast_var;
       
       inst_freq_syst_mean;           
       inst_freq_syst_var;
       inst_freq_diast_mean;
       inst_freq_diast_var;
      
       var_dist_s1_s2
    end
    
    methods
        function obj = calc_spi_feature(obj, pcg, fs, N, bin_width, N_mean_win, doPlot)
            [~, obj.spi, ~] = Spectral_Purity_Sliding_Window(pcg, fs, 0, N);
            
            obj.mean_spi = mean(obj.spi);
            
            obj.mean_spi_win = conv(abs(obj.spi), 1/N_mean_win*ones(N_mean_win, 1), 'valid');
            obj.mmean_spi = median(obj.mean_spi_win);
            
%             figure(40);
%             clf;
%             h(1) = subplot(2,1,1); 
%             plot(pcg./max(pcg));
%             h(2) = subplot(2,1,2); 
%             hold on;
%             plot([zeros(N_mean_win-1,1); obj.mean_spi_win]);
%             plot(obj.mean_spi*ones(length(obj.spi),1), 'r--');
%             plot([zeros(N_mean_win,1); obj.d_mean_spi_win*10 + obj.mean_spi]);
%             
%             linkaxes(h,'x');
%             pause;
            
            [obj.hist_spi, obj.hist_edges] = histcounts(obj.spi, 'BinWidth', bin_width, 'Normalization', 'probability');
            
%             [~, obj.max_hist_ind] = max(obj.hist_spi);
%             obj.max_hist_ind = (obj.max_hist_ind-1)*bin_width;
            obj.max_hist_ind = max(obj.hist_spi)/sum(obj.hist_spi);

            if doPlot
                figure(40);
                clf;
                subplot(2,1,1); 
                plot(pcg./max(pcg));
                subplot(2,1,2); 
                hold on;
                plot(obj.spi);
                plot(ones(length(obj.spi),1), 'r--');

                figure(41);
                clf;
                histogram(obj.spi, obj.hist_edges);
                
                pause;
            end
        end
        
        function obj = calc_denoising_and_spect_feat(obj, pcg, fs, doPlot)
            if doPlot
                figure(51)
                sgram(pcg,fs)
                %caxis([-60 0])
                % Undersampling to be more accurate in low frequencies!
                disp('subsampling by 3')
            end
            
            datas=decimate(pcg,3);
            
            if doPlot
                figure(52)
                sgram(datas,fs/3,'tfr',0.1)

                disp('Highpass filter')
            end
            
            [b,a] = butter(2,0.1,'high'); %22000/32*0.1=68
            yb=filtfilt(b,a,datas);
            
            % Gabor thresholding
            if doPlot
                disp('Gabor thresholding');
            end
       
            N=2000; 
            N=min(N,length(yb));
            
            gm = pgauss(N,0.1); %try 0.1
            M = floor(N/8);
            
            G=dgtreal(yb,gm,20,M);
            
            if doPlot
                figure(53)
                imagesc(log(abs(G)))
                caxis([-30 0])
            end
            
            % set to zero the small coefficients
            Gt=G;
            Gt(abs(Gt)<0.1*max(max(abs(Gt))))=0; %   0.1 for battement1
            % smallest time points to zero
            Gt(:,sum(abs(Gt).^2)<0.1)=0;
            
            obj.Gt=Gt;
            obj.Gtabs=abs(Gt);
            obj.Gtsumabs=sum(abs(Gt).^2);
            
            temp=sum(abs(Gt).^2); 
            temp(temp==0)= []; 
            temp=(temp-min(temp))/(max(temp)-min(temp));
            
            obj.Gtspar=length(temp)/length(abs(Gt));
            obj.Gtmean=mean(temp);
            obj.Gtvar=std(temp);
            obj.Gtmmax=mean(mean(abs(Gt)));
            
            if doPlot
                figure(54); plot(sum(abs(Gt).^2));
                figure(56)
                imagesc(log(abs(Gt)))
                % compute the inverse transform with dual window
                disp('computing inverse Gabor transform')
            end
            
            L=dgtlength(N,20,M);
            ga=gabdual(gm,20,M,L);
            ybitmp=idgtreal(Gt,ga,20,M);
            ybi=ybitmp(1:length(yb));
            
            if doPlot
                figure(57); clf; hold off;    
                plot(yb)
                hold on 
                plot(ybi)
                
                pause;
            end
        end
        
        function obj = calc_wavelet_feature(obj, pcg)
            deubachis = 'db8';
            level = 8;
            levelnum = 2^level;
            Wavelet_Nodes = 1:levelnum;

            T = wpdec(pcg ,level,deubachis);    % tree
            E = wenergy(T);                     % energy
            feats = (E(: , Wavelet_Nodes));
            obj.wav_feature = feats([3 4 6 21 24 27 29 30 31 32 81 256]).';
        end
        
        function obj = calc_qrs_en_coeff_feat(obj, pcg, fs, seg_states) 
            fsig = pcg;
            obj.qrs_en_coeff_mean = mean(fsig);
            obj.qrs_en_coeff_var = var(fsig);
            
            obj.qrs_en_coeff_syst_mean = mean(fsig(seg_states == 2));
            obj.qrs_en_coeff_syst_var = var(fsig(seg_states == 2));
            obj.qrs_en_coeff_diast_mean = mean(fsig(seg_states == 4));
            obj.qrs_en_coeff_diast_var = var(fsig(seg_states == 4));
            
%             figure(60);
%             hold on;
%             plot(fsig);
%             plot(seg_states);
%             pause;
        end
        
        function obj = calc_basic_features(obj, s2beg, s1end)
            obj.var_dist_s1_s2 = var(s2beg - s1end);
        end
        
        function obj = calc_instantanous_frequency(obj, pcg, fs, seg_states)
            [~, ~, ifreq] = IF_sig(pcg, fs);
            
            obj.inst_freq_syst_mean = median(ifreq(seg_states == 2));
            obj.inst_freq_syst_var = var(ifreq(seg_states == 2));
            obj.inst_freq_diast_mean = median(ifreq(seg_states == 4));
            obj.inst_freq_diast_var = var(ifreq(seg_states == 4));
            
%             figure(70);
%             clf;
%             hold on;
%             h1 = subplot(211);
%             hold on; 
%             plot(ifreq,'k');
%             plot([1,length(ifreq)],[obj.inst_freq_syst_mean, obj.inst_freq_syst_mean],'r');
%             plot([1,length(ifreq)],[obj.inst_freq_diast_mean, obj.inst_freq_diast_mean],'b');
%             plot([1,length(ifreq)],[obj.inst_freq_syst_mean-sqrt(obj.inst_freq_syst_var),obj.inst_freq_syst_mean-sqrt(obj.inst_freq_syst_var)],'g--');
%             plot([1,length(ifreq)],[obj.inst_freq_syst_mean+sqrt(obj.inst_freq_syst_var),obj.inst_freq_syst_mean+sqrt(obj.inst_freq_syst_var)],'g--');
%             
%             plot([1,length(ifreq)],[obj.inst_freq_diast_mean-sqrt(obj.inst_freq_diast_var),obj.inst_freq_diast_mean-sqrt(obj.inst_freq_diast_var)],'y--');
%             plot([1,length(ifreq)],[obj.inst_freq_diast_mean+sqrt(obj.inst_freq_diast_var),obj.inst_freq_diast_mean+sqrt(obj.inst_freq_diast_var)],'y--');
%             h2 = subplot(212); plot(pcg);
%             linkaxes([h1,h2],'x');
%             pause;
        end
        
        function [fm] = get_tot_tape_feature_matrix(obj,N_imp_ssa)
            var_num_max = var(obj.fm_b(1:N_imp_ssa,:), 0, 2);
            var_ratio_max = var(obj.fm_b(N_imp_ssa*1+1:N_imp_ssa*2,:), 0, 2);
            var_tot_max = var(sum(obj.fm_b(1:N_imp_ssa,:), 1));
            
            mean_num_max = mean(obj.fm_b(1:N_imp_ssa,:), 2);
            mean_ratio_max = mean(obj.fm_b(N_imp_ssa*1+1:N_imp_ssa*2,:), 2);
            mean_tot_max = mean(sum(obj.fm_b(1:N_imp_ssa,:), 1));
            
            fm = [obj.mmean_spi; obj.max_hist_ind; obj.Gtspar; obj.Gtmean; obj.Gtvar; obj.Gtmmax; ...
                obj.wav_feature; ...
                %obj.qrs_en_coeff_mean; obj.qrs_en_coeff_var; obj.qrs_en_coeff_syst_mean; obj.qrs_en_coeff_syst_var; obj.qrs_en_coeff_diast_mean; obj.qrs_en_coeff_diast_var;...
                mean_num_max; mean_ratio_max; mean_tot_max; var_num_max; var_ratio_max; var_tot_max; ...
                obj.var_dist_s1_s2; obj.inst_freq_syst_mean; obj.inst_freq_syst_var; obj.inst_freq_diast_mean; obj.inst_freq_diast_var;...
                obj.sum_bf; obj.pess_bf; obj.perc_an_bf; obj.perc_n_bf];
        end
    end
end

