classdef CalculationHandler
    %CalculationHandler Summary of this class goes here
    %   Detailed explanation goes here    
    properties
        config;
    end
    
    methods
        function c = CalculationHandler(conf)
            c.config = conf;
        end
        
        function [t] = load_tape(obj,file_name)
            t = load_tape_audio(file_name);
            if isvector(t)
                t = t(1);
            end
            % normalize
            [~,~,~,pcg]=Find_outliers(t.pcg,2);

            if(~isnan(pcg))
                t.pcg=pcg; 
            else
                display('found a wierd PCG');
            end
        end
        
        function [t] = segment_tape(obj, t, method, doPlot)
%             [~, t.fsig, ~] = QRS_Energy(t.pcg, obj.config.fs,1,0);
            
            % resample
            rs_pcg = resample(t.pcg, obj.config.seg_fs , obj.config.fs);
            
            [seg_ind, s1_end, s2_beg, s2_end, beat_len, ~] = pcg_segmentation(rs_pcg, [], obj.config.seg_fs, method, doPlot);
            t.seg_ind = floor(obj.config.fs/obj.config.seg_fs).*seg_ind;
            t.s1_end = floor(obj.config.fs/obj.config.seg_fs).*s1_end;
            t.s2_beg = floor(obj.config.fs/obj.config.seg_fs).*s2_beg;
            t.s2_end = floor(obj.config.fs/obj.config.seg_fs).*s2_end;
            t.beat_len = floor(obj.config.fs/obj.config.seg_fs).*beat_len;
            t.n_beats = length(t.seg_ind)-1;
            
            % generate assigned states
            assigned_states = any(and(bsxfun(@ge, 1:length(t.pcg), t.seg_ind(1:end-1)), bsxfun(@lt, 1:length(t.pcg), t.s1_end)));
            assigned_states = assigned_states + 2*any(and(bsxfun(@ge, 1:length(t.pcg), t.s1_end), bsxfun(@lt, 1:length(t.pcg), t.s2_beg)));
            assigned_states = assigned_states + 3*any(and(bsxfun(@ge, 1:length(t.pcg), t.s2_beg), bsxfun(@lt, 1:length(t.pcg), t.s2_end)));
            assigned_states = assigned_states + 4*any(and(bsxfun(@ge, 1:length(t.pcg), t.s2_end), bsxfun(@lt, 1:length(t.pcg), t.seg_ind(2:end))));
            t.assigned_states = assigned_states;
            
%             figure(11);
%             hold on;
%             plot(assigned_states);
%             plot(t.pcg);
%             hold off;
%             pause;
            
            if doPlot
                seg_logic = t.seg_ind;
                figure(10);
                clf;
%                 h(1) = subplot(2,1,1);
                hold on;
                plot(t.pcg/max(t.pcg));
                plot(t.fsig/max(t.fsig));
                stem(seg_logic, ones(size(seg_logic)), 'r-');
                legend('PCG signal', 'Beginning of beat');
                hold off;

%                 h(2) = subplot(2,1,2);
%                 hold on;
%                 plot(t.tm, t.ecg/max(t.ecg), 'k-');
%                 plot(t.tm, seg_logic, 'r-');
%                 hold off;

%                 linkaxes(h,'x');
                pause;
            end
        end
        
        function t_features = extract_features(obj, t, doPlot)
            if ~t.has_segmentation()
                error('Please perform segmentation before extracting features!');
            end
            
            fm_beats = zeros(t.n_beats, obj.config.N_important_ssa_decomp*8 + 74);
            
            seg_ind = t.seg_ind;
            
            for i = 1:t.n_beats
                b_feature = BeatFeatures;
                
                beat = t.pcg(seg_ind(i):seg_ind(i+1));

                b_feature = b_feature.calc_basic_features(beat, t.s1_end(i)-seg_ind(i), t.s2_beg(i)-seg_ind(i), t.s2_end(i)-seg_ind(i));
                b_feature = b_feature.calc_ssa(beat , obj.config.fs, obj.config.N_ssa_sig, obj.config.E_perc, obj.config.max_freq, obj.config.N_max, obj.config.min_width_max, t.s1_end(i)-seg_ind(i), t.s2_beg(i)-seg_ind(i), obj.config.N_env_samples, obj.config.spars_thresh, doPlot, i);
                b_feature = b_feature.calc_stat_features(beat, t.s1_end(i)-seg_ind(i), t.s2_beg(i)-seg_ind(i), t.s2_end(i)-seg_ind(i)); 
                fm_beats(i,:) = b_feature.get_feature_matrix(obj.config.N_important_ssa_decomp);
            end
            
%             persistent fm_beats_all;
%             if isempty(fm_beats_all)
%                 fm_beats_all = load('val_beat_v2_2.txt');
%             end
%             fm_beats = fm_beats_all(1:t.n_beats,1:(end-1));
%             fm_beats_all(1:t.n_beats,:) = [];
            
            % classify beats to calc Tape features based on beat
            % classification
            persistent BC_BAG;
            if (isempty(BC_BAG))
                load('BC_BAGC_subset2.mat');
            end
            
            f_ind =[11 26 35 36 37 40 41 44 54 104 107]; % beat features
            res = zeros(size(BC_BAG,2), size(fm_beats,1));
%             for i=1:size(BC_BAG,2)
%                  res(i,:)=predict(BC_BAG(i).ClassificationEnsemble, fm_beats(:,f_ind));              
%             end
            for i=1:size(BC_BAG,1)
                 res(i,:)=predict(BC_BAG{i,1}, fm_beats(:,f_ind));  
            end
           
            beat_res = -2*any(res == -1)+1;
            
            t_features = TapeFeatures;
            t_features.fm_b = fm_beats.';
            
            t_features.sum_bf = sum(beat_res);
            t_features.pess_bf = -2*any(beat_res == -1)+1;
            t_features.perc_an_bf = sum(beat_res == 1)./length(beat_res);
            t_features.perc_n_bf = sum(beat_res == -1)./length(beat_res);
            
            t_features = t_features.calc_instantanous_frequency(t.pcg, obj.config.fs, t.assigned_states);
            t_features = t_features.calc_spi_feature(t.pcg, obj.config.fs, obj.config.N_win_spi, obj.config.bin_width_spi, obj.config.N_mean_win, doPlot);
            t_features = t_features.calc_denoising_and_spect_feat(t.pcg, obj.config.fs, doPlot);
            t_features = t_features.calc_wavelet_feature(t.pcg);
%             t_features = t_features.calc_qrs_en_coeff_feat(t.fsig, obj.config.fs, t.assigned_states);
            t_features = t_features.calc_basic_features(t.s2_beg, t.s1_end);
        end
        
        function result = classify(obj, t_features)
            fm = t_features.get_tot_tape_feature_matrix(obj.config.N_important_ssa_decomp);
%             fm = fm([1 2 4 5 7 8 9 17 18 41 43 44 45 46 48]);
            fm = fm([2 4 8 9 18 42 43 44 45 46 47 48]); % new features (without Gabor and maxima)
            
            persistent SSC_ADB; 
            if (isempty(SSC_ADB))
                load('SSC_ADB_tb_incl_f.mat');
            end
            
            for i=1:size(SSC_ADB,2)
                A(i)=predict(SSC_ADB(i).ClassificationEnsemble, fm.');      
            end
            
            % pessimistic
            result = -2*any(A == -1)+1;
        end
    end
end

