function op = hybrid(ip, tsm_factor, win_len, P)
% A hybrid time-scale modification algorithm based upon A Hybrid Time-Frequency Domain Approach to Audio Time-Scale Modification, Dorran, David; Lawlor, Bob; Coyle, Eugene, Journal of The Audio Engineering Society, vol: 54, issue: 1/2, pages: 21 - 31
%
% David Dorran, Audio Research Group, Dublin Institute of Technology
% david.dorran@dit.ie
% http://eleceng.dit.ie/dorran
% http://eleceng.dit.ie/arg
%

% make sure input is mono and transpose if necessary 
[r, c] = size(ip);
if r > 1
    ip = ip';
end;    
[r, c] = size(ip);
if r > 1
    disp('Note :this version only works on mono signals');
    op = [];
    return
end;

% initialize the values of analysis_frame_length, analysis_window_offset, synthesis_frame_offset and length_of_overlap
hop_size = win_len/4;
max_phase_deviation = 0.33;
hansquared = hanning(win_len).^2;
amp_mod_wave = [hansquared(3*hop_size +1:4*hop_size)+hansquared(2*hop_size +1:3*hop_size)+hansquared(1*hop_size +1:2*hop_size); hansquared(3*hop_size +1:4*hop_size)+hansquared(2*hop_size +1:3*hop_size); hansquared(3*hop_size +1:4*hop_size)]*2/3;
max_ip = max(ip);

desired_tsm_len = round(length(ip)*tsm_factor);
stationary_length = P*1.5;%win_len/4; % found this to a reasonable value for the stationary length.
Lmax = win_len/2;

analysis_window_offset = round((stationary_length - P)/abs(tsm_factor - 1)); % this equation was derived.
synthesis_window_offset = round(analysis_window_offset * tsm_factor);
analysis_frame_length = round(Lmax + synthesis_window_offset);
number_of_analysis_frames = floor((length(ip)- analysis_frame_length)/analysis_window_offset);

if number_of_analysis_frames < 2 %not much time-scaling being done just return the input
    op = ip;
    return;
end;

%the next two lines just ensure that the last frame finishes at the very end of the signal (not essential)
zpad = zeros(1, (number_of_analysis_frames*analysis_window_offset) + analysis_frame_length - length(ip) + win_len*2);
ip = [ip zpad];

%initialize the output
op = zeros(1, desired_tsm_len);

%initialize the first output frame
op(1 : analysis_frame_length+win_len) = ip(1 : analysis_frame_length+win_len);

min_overlap = round(Lmax - P); %ensure that there is some minimum overlap
count = 0;

prev_opt_overlap_length = Lmax;  
% Loop through all the analysis frames
for m = 1 : number_of_analysis_frames
    disp(strcat(num2str(m/number_of_analysis_frames*100), '% complete'))
    %grab the mth input frame
    ip_frame = ip(analysis_window_offset * m : (analysis_window_offset * m) + analysis_frame_length - 1);
    
    %grab the maximum overlapping segments from the input frame and the current output
    seg_1 = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax : round(synthesis_window_offset*(m-1))+analysis_frame_length -1);
    seg_2 = ip_frame(1: Lmax);
    
    %compute the correlation of these segments
    correlation = xcorr(seg_2, seg_1,'coeff');
    
    %Find the best point to overlap (opt_overlap_length) making sure not to exceed the maximum or go below the minimum overlap.
    correlation(length(correlation) - Lmax -1: length(correlation)) = -100;
    correlation(1: min_overlap) = -100;    
    [max_correlation, opt_overlap_length] = max(correlation);
    
    %this next block of code choses the maximum overlap from a number of possible overlaps - this helps to reduce the duration of the segment which is discarded/repeated which in turn helps ensure that the next_op_win is more likely to be similar to next_win
    orig_opt_overlap_length = opt_overlap_length;
    correlation_peaks = locate_peaks(correlation);
    new_correlation = zeros(1, length(correlation));
    new_correlation(correlation_peaks) = correlation(correlation_peaks);
    valid_overlaps = find(new_correlation > max_correlation*0.85);
    opt_overlap_length =  max(valid_overlaps);

    if(max_correlation < 0)
        opt_overlap_length = Lmax/2; %this line is a bug workaround for the case of zero inputs
    end;
    offset = Lmax - opt_overlap_length;

    %next_op_win = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - 1 - round(opt_overlap_length/2)-2*hop_size : round(synthesis_window_offset*(m-1))+analysis_frame_length -1 - round(opt_overlap_length/2)-2*hop_size + win_len-1).*hanning(win_len)';
    if(round(synthesis_window_offset*(m-1))+analysis_frame_length - 1 - round(opt_overlap_length/2)-2*hop_size < 1)
        next_op_win = op(1 : round(synthesis_window_offset*(m-1))+analysis_frame_length -1 - round(opt_overlap_length/2)-2*hop_size + win_len-1);
        next_op_win = [zeros(1, win_len - length(next_op_win)) next_op_win].*hanning(win_len)';
    else
        next_op_win = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - 1 - round(opt_overlap_length/2)-2*hop_size : round(synthesis_window_offset*(m-1))+analysis_frame_length -1 - round(opt_overlap_length/2)-2*hop_size + win_len-1).*hanning(win_len)';
    end
    
    phase_estimates = angle(fft(next_op_win));
  
    %initialize the next_win variable
    %next_win = ip(analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size : analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size + win_len- 1).*hanning(win_len)';
    if(analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size < 1)
        next_win = ip(1 : analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size + win_len- 1);
        next_win = [zeros(1, win_len - length(next_win)) next_win].*hanning(win_len)';
    else
        next_win = ip(analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size : analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size + win_len- 1).*hanning(win_len)';
    end
    [indices, num_zeros] = zero_crossings(correlation(min_overlap +1 : length(correlation) - Lmax-2));
    
    start_nrg_seg_1 = sum(abs(seg_1(1:Lmax/2)));
    end_nrg_seg_1 = sum(abs(seg_1(Lmax/2:length(seg_1))));
    
    %if the phase estimates don't look like they are any good at all then use the original phases.
    temp_var = 0;
    if ((num_zeros/length(correlation(min_overlap +1 : length(correlation) - Lmax-2)) > 0.2) & max_correlation < 0.4) | (num_zeros/length(correlation(min_overlap +1 : length(correlation) - Lmax-2)) == 0 | max_correlation < 0.3 | end_nrg_seg_1/start_nrg_seg_1 > 2)
       phase_estimates  = angle(fft(next_win)); 
       temp_var = 1;
    end 
        
    new_ip_frame = zeros(1, analysis_frame_length+win_len*2);
    dont_stop_looping = 1;
    cur_index = 1;
    while dont_stop_looping == 1
        orig_phase_estimates = phase_estimates;
        cur_win = next_win;
        cur_mags = abs(fft(cur_win));
        cur_phases = angle(fft(cur_win));
               
        %locate the peaks
        pk_indices = [];
        pk_indices  =  locate_peaks(cur_mags(1:win_len/2));
        
        %update phase of channels in region of influence
        start_point = 1; %initialize the starting point of the region of influence
        pk_indices = [pk_indices win_len/2 win_len/2];
        for k = 1: length(pk_indices) -1
            peak = pk_indices(k);
            
            %---+++ Following lines shift the phase towards the original phase
            if(cur_index ~= 1)
                phase_diff = princarg(phase_estimates(peak) - cur_phases(peak));
                
                if (abs(phase_diff) < max_phase_deviation)
                    phase_estimates(peak) = cur_phases(peak);
                else
                    if(phase_diff < 0)
                        phase_estimates(peak) = phase_estimates(peak) + max_phase_deviation;
                    else
                        phase_estimates(peak) = phase_estimates(peak) - max_phase_deviation;
                    end;
                end;
            end
            %---+++ previous lines shift the phase towards the original phase
            
            % first calculate angle of rotation
            angle_rotation = phase_estimates(peak) - cur_phases(peak);                
            
            next_peak = pk_indices(k+1);
            end_point = round((peak + next_peak)/2);
            ri_indices = [start_point : peak-1, (peak+1) : end_point]; %indices of the region of influence
            phase_estimates(ri_indices) = angle_rotation + cur_phases(ri_indices);
            
            start_point = end_point + 1;
        end;  
        phase_estimates = [phase_estimates(1:win_len/2 +1) fliplr(phase_estimates(2:win_len/2)*-1)];
        
        new_win = real(ifft(cur_mags.*exp(i*phase_estimates))).*hanning(win_len)'.*2/3;
        new_ip_frame((cur_index-1)*hop_size +1: (cur_index-1)*hop_size +win_len) = new_ip_frame((cur_index-1)*hop_size +1: (cur_index-1)*hop_size +win_len) + new_win;

        next_win = ip(analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size + hop_size*cur_index: analysis_window_offset * m + round(opt_overlap_length/2) -2*hop_size + hop_size*cur_index + win_len- 1).*hanning(win_len)';
        cur_index = cur_index +1;
        if(analysis_window_offset * m + opt_overlap_length + hop_size*cur_index- 1 > (analysis_window_offset * m) + analysis_frame_length - 1 + win_len)
            dont_stop_looping = 0;
        end
        phaseInc = angle(fft(next_win)) - cur_phases;
        phase_estimates = phase_estimates + phaseInc;
    end
    sum(abs(cur_phases - phase_estimates + phaseInc)) ;  
        
    % append mth analysis frame to the current synthesised output
    end_seg = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - round(opt_overlap_length/2)-2*hop_size:round(synthesis_window_offset*(m-1))+analysis_frame_length -1 - round(opt_overlap_length/2)+hop_size);    
    ov_seg = [end_seg.*amp_mod_wave'+new_ip_frame(1:hop_size*3) new_ip_frame(hop_size*3+1:length(new_ip_frame))];
    op(round(synthesis_window_offset*(m-1))+analysis_frame_length - round(opt_overlap_length/2) - 2*hop_size: round(synthesis_window_offset*(m-1))+analysis_frame_length - round(opt_overlap_length/2) - 2*hop_size + length(ov_seg) -1) = ov_seg;
    prev_opt_overlap_length = opt_overlap_length;
end; % end of for loop
% END

%--------------------------------------------------------------------------
function indices = locate_peaks(ip)
%function to find peaks
%a peak is any sample which is greater than its two nearest neighbours
	index = 1;
	num = 2;
    indices = [];
	for k = 1 + num : length(ip) - num
		seg = ip(k-num:k+num);
		[val, max_index] = max(seg);
		if max_index == num + 1
			indices(index) = k;
			index = index + 1;
		end;
	end;