% training_params.{fs, band_boundaries, window_size, overlapping_windows}
% assume 16 bands
%   01 - 16: band_mean_t - mean power of each band over time
%   17 - 32: band_std_t - std power of each band over time
%   33 - 34: mpb - mean most powerful band of 4 , most powerful band of 16           
%   35:      mpb_std_t - std of most powerful band of 16
%   36:      band_width_50 - width of smallest set of bands containing 50 of power
%   37:      nsb_mean_t - mean of number of bands with greater than 25% power
%   38:      nsb_std_t - std of number of bands with greater than 25% power
%   39:      power_std - std of total power over time
%   40:      power_peaks - number of power peaks (below 30% to above 70%)
%   41:      periodicity - regularity of power peaks
%   42:      power range
%   43:      mean power
%   44 - 47: band_mean_t - reduced to 4 bands
%   48 - 51: band_std_t - reduced to 4 bands
%   52 - 138: time-dependent frequency information at different resolutions
function [features, window_length, fft_data, fft_pow] = extract_local_features_window(data, params)

    % fft_data(num_bands, num_times)    
    fft_data = compute_windowed_fft(data, params.fs, params.band_boundaries, ...
        params.window_size, params.overlapping_windows);
    
    fft_pow = sum(fft_data, 1);
    for i = 1:size(fft_data, 2)
        if fft_pow(i) ~= 0
            fft_data(:, i) = fft_data(:, i) / fft_pow(i);            
        end        
    end 
    window_length = size(fft_data, 2);
    
    features = zeros(138, 1);
    
    % mean and std power
    features(1:16) = mean(fft_data, 2);
    features(17:32) = std(fft_data, 0, 2);
    
    % most powerful band (perhaps should change mean to mode)
    [tmp, mpb_16] = max(fft_data, [], 1);
    mpb_4 = mpb_16/4;
    features(33) = ceil(mean(mpb_4));
    features(34) = round(mean(mpb_16));
    features(35) = std(mpb_16);
    
    % band width
    band_pow = sum(fft_data, 2);
    total_pow = sum(band_pow);
    found = 0;
    for width = 1:16
        for b = 1:16-width+1
            if sum(band_pow(b:b+width-1)/total_pow) > 0.5
                found = 1;                
                break;
            end
        end
        if found
            features(36) = width;
            break;
        end
    end
    
    % bands over 25% power
    num_bands_over_25 = sum(fft_data > 0.25, 1);
    features(37) = mean(num_bands_over_25);
    features(38) = std(num_bands_over_25);
    
    % std power
    features(39) = std(fft_pow);
    
    % power peaks
    min_power = min(fft_pow);
    max_power = max(fft_pow);
    power_70 = 0.7*(max_power-min_power)+min_power;
    power_30 = 0.3*(max_power-min_power)+min_power;    
    had_low = 1;
    had_high = 0;        
    for t = 1:length(fft_pow)
        curr_state = (fft_pow(t) > power_30) + (fft_pow(t) > power_70);
        if curr_state < 1
            had_low = 1;
            had_high = 0;
        elseif curr_state > 1
            had_high = 1;
            if had_low
                had_low = 0;
                features(40) = features(40) + 1;
                peak_loc(features(40)) = t;
            end
        end
    end        
    
    if features(40) > 1        
        features(41) = mean(peak_loc(2:end)-peak_loc(1:end-1));
        features(41) = mean(abs((peak_loc(2:end)-peak_loc(1:end-1))-features(41)));        
    else
        features(41) = length(fft_pow)/4;
    end
    features(42) = max_power-min_power;
    features(43) = mean(fft_pow);
    
    % coarser resolution frequencies
    for f = 44:47
        i1 = (f-44)*4+1;
        i2 = i1+3;
        features(f) = mean(features(i1:i2));
        features(f+4) = std(sum(fft_data(i1:i2, :), 1), 0 , 2);
    end
    
    features(52:52+87-1) = compute_local_features_window(fft_data, fft_pow);    
    
    
    