classdef BeatFeatures
    %FEATURES Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        % ssa features
        ssa_freq;
        ssa_var;
        ssa_energy;
        ssa_en_var;
        
        dist_max;
        ratio_max;
        
        mean_if;
        var_if;
        
        width_s1;
        width_s2;
        ampl_s1;
        ampl_s2;
        ratio_ampl;
        dist_s1_s2;
        
        stat_feat;
    end
    
    methods
        function obj = calc_ssa(obj, bsig, fs, L, energy_thresh, max_freq, N_max, min_width_max, s1end, s2beg, N_env_samples, spars_thresh, doPlot, num_beat)
            [sig, freq, f_var] = SSA_decomposition(bsig, L, fs, 0);
            
            v_sig = sig./repmat(max(sig),[size(sig,1),1]);
            [freq, ind] = sort(freq);
            f_var = f_var(ind);
            sig = sig(:,ind);
            
            % get rid of low energy content
            ind = obj.calc_high_energy_ind(sig,bsig,energy_thresh);
            freq(~ind) = 0;
            f_var(~ind) = 0;
            
            % hf energy mean and var
            obj.ssa_energy = mean(v_sig.^2).';
            obj.ssa_en_var = iqr(v_sig.^2).';
             
            obj.ssa_freq = freq;
            obj.ssa_var = f_var;
            
            [obj.dist_max, obj.ratio_max, obj.mean_if, obj.var_if, ~] = obj.get_features_from_decompositions(sig, ind, N_max, min_width_max, s1end, s2beg, fs);
            
            if doPlot
                p_sig = sig(:,ind);
                p_dist = obj.dist_max(:,ind);
                p_ratio = obj.ratio_max(:,ind);
                p_env = obj.env(:,ind);
                
                figure(20);
                hold on;
                clf;
                subplot(size(p_sig,2)+2,1,1), hold on, plot(bsig./max(bsig), 'k-'), hold off;
                title(['Decomposition of beat ', num2str(num_beat)]);
                subplot(size(p_sig,2)+2,1,2), hold on, plot(obj.w_inst_freq, 'k.-'), hold off;
                for i=1:size(p_sig,2)
                    subplot(size(p_sig,2)+2,1,i+2)
                    hold on;
                    plot(p_sig(:,i),'b-');
                    plot(p_dist(:,i), p_ratio(:,i), 'r*');
                    plot(obj.env_ind, p_env(:,i), 'g-');
                    hold off;
                end
                pause;
            end
        end
        
        function ind = calc_high_energy_ind(obj, sig, bsig, energy_thresh)
            e_sig = sig;
            for i = 2:size(sig,2)
                e_sig(:,i:end) = e_sig(:,i:end) + sig(:,1:end-i+1);
            end
            e_sig = mean(e_sig.^2)./mean(bsig.^2);
            
            ind = e_sig <= energy_thresh;
        end
        
        function [dist_to_max, ratio_max, mean_if, var_if, sparsity] = get_features_from_decompositions(obj, sig, h_en_ind, N_max, width_max, s1end, s2beg, fs)
            h_en_ind = find(h_en_ind == 1);
            
            dist_to_max = zeros(size(sig, 2),1);
            ratio_max = zeros(size(sig, 2),1);
            mean_if = zeros(size(sig,2),1);
            var_if = zeros(size(sig,2),1);
            sparsity = zeros(size(sig, 2),1);
            
            for i = 1:length(h_en_ind)
                % extract distance to N_max highest maxima
                decomp = sig(:,h_en_ind(i));
                decomp = decomp./max(decomp);
                
%                 [pks, ind] = findpeaks(decomp, 'MinPeakWidth', width_max, 'SortStr', 'descend', 'NPeaks', N_max);
% %                 [pks, ind] = findpeaks(decomp, 'MinPeakWidth', width_max, 'MinPeakHeight', 0.25, 'SortStr', 'descend', 'NPeaks', N_max);
% 
%                 if length(ind) < N_max
%                     ind = [ind(:); zeros(N_max - length(ind), 1)];
%                     pks = [pks(:); zeros(N_max - length(pks), 1)];
%                 end
%                 dist_to_max(h_en_ind(i)) = sum(ind(:) > s1end & ind(:) < s2beg);
%                 ratio_max(h_en_ind(i)) = sum(pks(ind(:) > s1end & ind(:) < s2beg));
                
                % extract IF of decomposition
                [~, ~, i_f] = IF_sig(decomp, fs);
                mean_if(h_en_ind(i)) = mean(i_f);
                var_if(h_en_ind(i)) = var(i_f);
            end
        end
        
        function obj = calc_basic_features(obj, bsig, s1end, s2beg, s2end)
        
            obj.width_s1 = s1end(1);
            obj.width_s2 = s2end(1) - s2beg(1) + 1;
            obj.ampl_s1 = max( bsig(1:s1end(1)) );
            obj.ampl_s2 = max( bsig(s2beg(1):s2end(1)) );
            obj.ratio_ampl = obj.ampl_s2/obj.ampl_s1;
            obj.dist_s1_s2 = s2beg(1) - s1end(1);
            
%             figure(61);
%             clf;
%             hold on;
%             plot(bsig);
%             plot([1,s1end(1)], [1,1], 'r-*');
%             plot([s2beg(1),s2end(1)], [1,1], 'r-*');
%             
%             pause;
        end
        
        function obj = calc_stat_features(obj, bsig, s1end, s2beg, s2end)
            bfeat1=sum(bsig(1:s1end).^2)/s1end;
            bfeat2=sum(bsig(1:s2beg).^2)/s2beg;
            bfeat3=sum(bsig(1:s2end).^2)/s2end;
            bfeat4=sum(bsig.^2)/length(bsig);           
            bfeat5=sum(bsig(s1end:s2beg).^2)/numel(s1end:s2beg);       
            bfeat6=sum(bsig(s1end:s2end).^2)/numel(s1end:s2end);       
            bfeat7=sum(bsig(s1end:end).^2)/numel(s1end:length(bsig));
            bfeat8=sum(bsig(s2beg:s2end).^2)/numel(s2beg:s2end);        
            bfeat9=sum(bsig(s2beg:end).^2)/numel(s2beg:length(bsig));
            bfeat10=sum(bsig(s2end:end).^2)/numel(s2end:length(bsig));
            %max min std mean median 
            bfeatm1=median(bsig(1:s1end));
            bfeatm2=median(bsig(1:s2beg));
            bfeatm3=median(bsig(1:s2end));
            bfeatm4=median(bsig);           
            bfeatm5=median(bsig(s1end:s2beg));       
            bfeatm6=median(bsig(s1end:s2end));       
            bfeatm7=median(bsig(s1end:end));
            bfeatm8=median(bsig(s2beg:s2end));
            bfeatm9=median(bsig(s2beg:end));
            bfeatm10=median(bsig(s2end:end));

            bfeatm11=mean(bsig(1:s1end));
            bfeatm12=mean(bsig(1:s2beg));
            bfeatm13=mean(bsig(1:s2end));
            bfeatm14=mean(bsig);           
            bfeatm15=mean(bsig(s1end:s2beg));       
            bfeatm16=mean(bsig(s1end:s2end));       
            bfeatm17=mean(bsig(s1end:end));
            bfeatm18=mean(bsig(s2beg:s2end));
            bfeatm19=mean(bsig(s2beg:end));
            bfeatm20=mean(bsig(s2end:end));

            bfeatm21=min(bsig(1:s1end));
            bfeatm22=min(bsig(1:s2beg));
            bfeatm23=min(bsig(1:s2end));
            bfeatm24=min(bsig);           
            bfeatm25=min(bsig(s1end:s2beg));       
            bfeatm26=min(bsig(s1end:s2end));    
            bfeatm27=min(bsig(s1end:end));
            bfeatm28=min(bsig(s2beg:s2end));
            bfeatm29=min(bsig(s2beg:end));
            bfeatm30=min(bsig(s2end:end));

            bfeatm31=max(bsig(1:s1end));
            bfeatm32=max(bsig(1:s2beg));
            bfeatm33=max(bsig(1:s2end));
            bfeatm34=max(bsig);           
            bfeatm35=max(bsig(s1end:s2beg));       
            bfeatm36=max(bsig(s1end:s2end));       
            bfeatm37=max(bsig(s1end:end));
            bfeatm38=max(bsig(s2beg:s2end));
            bfeatm39=max(bsig(s2beg:end));
            bfeatm40=max(bsig(s2end:end));

            bfeatm41=std(bsig(1:s1end));
            bfeatm42=std(bsig(1:s2beg));
            bfeatm43=std(bsig(1:s2end));
            bfeatm44=std(bsig);           
            bfeatm45=std(bsig(s1end:s2beg));       
            bfeatm46=std(bsig(s1end:s2end));       
            bfeatm47=std(bsig(s1end:end));
            bfeatm48=std(bsig(s2beg:s2end));
            bfeatm49=std(bsig(s2beg:end));
            bfeatm50=std(bsig(s2end:end));
                
            bfeat11=numel(1:s1end);
            bfeat12=numel(1:s2beg);
            bfeat13=numel(1:s2end);
            bfeat14=length(bsig); 
            bfeat15=numel(s1end:s2beg); 
            bfeat16=numel(s1end:s2end);
            bfeat17=numel(s1end:length(bsig));
            bfeat18=numel(s2beg:s2end);   
            bfeat19=numel(s2beg:length(bsig)); 
            bfeat20=numel(s2end:length(bsig));
            
            obj.stat_feat = [bfeat1 bfeat2 bfeat3 bfeat4 bfeat5 bfeat6 bfeat7 bfeat8 bfeat9 bfeat10 ...
                         bfeatm1 bfeatm2 bfeatm3 bfeatm4 bfeatm5 bfeatm6 bfeatm7 bfeatm8 bfeatm9 bfeatm10 ...
                         bfeatm11 bfeatm12 bfeatm13 bfeatm14 bfeatm15 bfeatm16 bfeatm17 bfeatm18 bfeatm19 bfeatm20 ...
                         bfeatm21 bfeatm22 bfeatm23 bfeatm24 bfeatm25 bfeatm26 bfeatm27 bfeatm28 bfeatm29 bfeatm30 ...
                         bfeatm31 bfeatm32 bfeatm33 bfeatm34 bfeatm35 bfeatm36 bfeatm37 bfeatm38 bfeatm39 bfeatm40 ...
                         bfeatm41 bfeatm42 bfeatm43 bfeatm44 bfeatm45 bfeatm46 bfeatm47 bfeatm48 bfeatm49 bfeatm50 ...
                         bfeat11 bfeat12 bfeat13 bfeat14 bfeat15 bfeat16 bfeat17 bfeat18 bfeat19 bfeat20]; 
        end
        
        function fm_b = get_feature_matrix(obj,N_imp_ssa)
            fm_b = [obj.dist_max(1:N_imp_ssa); obj.ratio_max(1:N_imp_ssa); obj.ssa_freq(1:N_imp_ssa); obj.ssa_var(1:N_imp_ssa); obj.ssa_energy(1:N_imp_ssa); obj.ssa_en_var(1:N_imp_ssa); ...
                obj.width_s1; obj.width_s2; obj.dist_s1_s2; obj.ratio_ampl; obj.mean_if(1:N_imp_ssa); obj.var_if(1:N_imp_ssa); obj.stat_feat.'];

              % for the entry we are interested only in tape features
%               fm_b = [obj.dist_max(1:N_imp_ssa); obj.ratio_max(1:N_imp_ssa); zeros(5*N_imp_ssa,1); ...
%                   0; obj.width_s1; obj.width_s2; obj.dist_s1_s2; zeros(15,1)];
              
        end
    end
end

