function [finalOutput] = Shapee(file1, file2, outputFile, frameSize, overlap, omega)
%Shapee(file1, file2, outputFile, frameSize, overlap, omega)
%takes file1 and file2
%and combines them according to Christopher Penrose' Shapee
%algorithm. The combination is output as a wav file (outputFile).
%omega is the width (in bins) of the shaping region.
%Generally this should be kept close to the main lobe width
%of the window for best results
%MATLAB version Jez Wells Nov 2004

NumShapingRegions = floor((frameSize-1)/omega)-1;

frameCount = 1;

%find length of first file file and number of channels
params = wavread(file1,'size');
samples1 = params(1,1);
channels1 = params(1,2);
%do the same for the second file
params = wavread(file2,'size');
samples2 = params(1,1);
channels2 = params(1,2);
if channels1 ~= channels2
    error('Input files do not have the same number of channels');
end
if samples1 ~= samples2
    fprintf('Input files are not the same length.\nOutput file will be the same length as shortest input file.\n');
    samples = min([samples1,samples2]);
end
samples = min([samples1,samples2]);
%ensure our output matrix will be big enough to accomodate a whole number of
%frames by rounding up
outSamples = (ceil(samples/frameSize) + 1) * frameSize;
finalOutput = zeros(outSamples, channels1);

%generate window
Hann = sqrt(hann(frameSize));
% set ranges for shaping regions
lower = (omega*[0:NumShapingRegions])+1;
upper = lower+omega;

% read file
[input1, Fs1, bits1] = wavread(file1);
[input2, Fs2, bits2] = wavread(file2);
if Fs1~=Fs2
    error('Sample rates of input files are not the same');
end
if bits1~=bits2 && frameCount==1
    bits1
    bits2
    fprintf('Bit depth of files is not the same. Output bit depth will be the highest out of the two.\n');
    bits = max([bits1, bits2]);
end
bits = max([bits1, bits2]);
% get inputs to same length

    
while frameCount < samples
    outPointer = frameCount + (frameSize-1);
    if outPointer > samples
        outPointer = samples;
        siz = outPointer - frameCount + 1;
    end
    ip1 = input1(frameCount:outPointer,:);
    ip2 = input2(frameCount:outPointer,:);
    
    if length(ip1) < frameSize
        %calculate difference between input and frame length
        diff = frameSize - length(ip1);
        %append zeros
        oo = zeros(diff,2);
        ip1 = [rot90(ip1), rot90(oo)];
        %rotate back after appending zeros
        ip1 = rot90(ip1,-1);
    end
    if length(ip2) < frameSize
        %calculate difference between input and frame length
        diff = frameSize - length(ip2);
        %append zeros
        oo = zeros(diff,2);
        ip2 = [rot90(ip2), rot90(oo)];
        %rotate back after appending zeros
        ip2 = rot90(ip2,-1);
    end
    %place individual channels in separate vectors and apply window
    left1 = ip1(:,1).*Hann;  
    right1 = ip1(:,2).*Hann;
    left2 = ip2(:,1).*Hann;  
    right2 = ip2(:,2).*Hann;
    
    %perform fft
    leftFFT1 = fft(left1, frameSize);
    rightFFT1 = fft(right1, frameSize);
    leftFFT2 = fft(left2, frameSize);
    rightFFT2 = fft(right2, frameSize);
    %extract magnitudes and phases for each bin
    leftMag1 = abs(leftFFT1);
    rightMag1 = abs(rightFFT1);
    %leftPhase1 = angle(leftFFT1);
    %rightPhase1 = angle(rightFFT1);
    leftMag2 = abs(leftFFT2);
    rightMag2 = abs(rightFFT2);
    leftPhase2 = angle(leftFFT2);
    rightPhase2 = angle(rightFFT2);
    
    %do something below here
    %to create your own cross synthesizer
    
    cumAmpLeft = [0;cumsum(leftMag1)];
    cumAmpRight = [0;cumsum(rightMag1)];
    cumFreqLeft = [0;cumsum(leftMag2)];
    cumFreqRight = [0;cumsum(rightMag2)];
  
    sigmaAmpMagLeft = cumAmpLeft(upper)-cumAmpLeft(lower);
    sigmaAmpMagRight = cumAmpRight(upper)-cumAmpRight(lower);
    sigmaAmpFreqLeft = cumFreqLeft(upper)-cumFreqLeft(lower);
    sigmaAmpFreqRight = cumFreqRight(upper)-cumFreqRight(lower);
    
    sigmaAmpFreqLeft(find(sigmaAmpFreqLeft==0)) = realmin;
    sigmaAmpFreqRight(find(sigmaAmpFreqRight==0)) = realmin;
    
    MagMultiplierLeft = sigmaAmpMagLeft./sigmaAmpFreqLeft;
    MagMultiplierRight = sigmaAmpMagRight./sigmaAmpFreqRight; 
    
    for n = 1:NumShapingRegions+1
        leftMag1(lower(n):upper(n)-1) = leftMag2(lower(n):upper(n)-1)* MagMultiplierLeft(n);
        rightMag1(lower(n):upper(n)-1)  = rightMag2(lower(n):upper(n)-1) * MagMultiplierRight(n);
    end
    
    %phase for frequency reference is phase of second input file
        
    
    
    %when you're done modifying
    %resynthesise
    %convert mag and phase back to complex numbers
    %we'll re-use some of the input arrays to save on memory
    leftReal1 = leftMag1 .* cos(leftPhase2);
    leftImag1 = leftMag1 .* sin(leftPhase2);
    rightReal1 = rightMag1 .* cos(rightPhase2);
    rightImag1 = rightMag1 .* sin(rightPhase2);
    leftFFT1 = complex(leftReal1,leftImag1);
    rightFFT1 = complex(rightReal1, rightImag1);
    leftOutput = ifft(leftFFT1, 'symmetric');    %inverse FFT of data with complex conjugate symmetry is faster
    rightOutput = ifft(rightFFT1, 'symmetric');
    %window finalOutput to remove any clicks from discontinuities at frame boundaries
    output(:,1) = leftOutput.*Hann;
    output(:,2) = rightOutput.*Hann;
    %place finalOutput in single matrix to be written to stereo .wav file    
    finalOutput(frameCount:(frameCount + (frameSize-1)),:) =  finalOutput(frameCount:(frameCount + (frameSize-1)),:) + output(1:(frameSize),:);
    frameCount = frameCount + (frameSize/overlap);
end
m1 = max(abs(finalOutput(:,1)));
m2 = max(abs(finalOutput(:,2)));
m = max([m1,m2]);
m = 0.99/m; % avoids "output is clipped" error message in wavwrite
finalOutput(:,1) = finalOutput(:,1)* m;
finalOutput(:,2) = finalOutput(:,2)* m;
%write entire finalOutput to wav file
wavwrite(finalOutput,Fs1,bits,outputFile);
%play the finalOutput 
wavplay(finalOutput, Fs1);