function y = distort(x, window, overlap, f, separation, order, c, fs, ftype, synthfn)
%y = distort(x, window, overlap, f, separation, order, c, fs, ftype)
%takes in an input x, distorts it, output is y
%window is a column vector that windows the output of each synthesised block
%  for example, use a triangle window and an overlap of 2 to have the
%  amplitudes of each frequency bin linearly interpolate
%the size of each window determines the length of the filter block
%overlap is the overlapping factor between each window
%f is a vector of separation frequencies for the filters
%separation is a number between 0 and 1 that determines how far each filter
%is separated from each other.  At separation = 1 they are aligned at their
%3dB points (or specified cutoff for cheby).
%order is the order of the filter, must be a multiple of 2
%c is an ascending vector for zeroing bins if we think they are getting
%power from a nearby bin.
%fs is the sampling frequency
%ftype is a string:
%  'b' denotes a butterworth filter
%  'cR' denotes a type 1 chebychev filter with ripple R dBs.
%synthfn is the synthesis function to apply to each frame.  the distort
%functionality should be incluced in synthfn
%  (may be omitted for now)
%it should take in the magnitude of frequencies vector and return the
%distorted magnitude of frequencies vector.
%The code is based on my shiftpitch function

%following is a list of design parameters that can be changed:
%the operation can be modified by changing the window, window length,
%overlap, distortion function, f (frequency boundary list), seperation,
%order of the butterworth filters, and c (dynamic cutoff coefficients).

%we need to define the frequency bins that we are going to use.
%each bin will be a bandpass filter between 2 consecutive elements of f.
%f here is in actual frequency, we will normalize.
%f is expected to be in increasing order.
%the number of filters, and hence the number of bins, is length(f) - 1.

%drop D on a guitar is 73.4 Hz so go a little lower than this
%high E 24th fret is 1318.4, to account for harmonics should probably go as
%high as 10KHz, but can experiment with this later

%this is a relatively conservative spacing
%f = linspace(200, 1000, 1000);
%f = [900 920];

%I want to be able to change the seperation of the filters.  If we use the
%frequencies as bounds the 3dB points will align.  While this is half the
%power, the peak voltage will be sqrt(power) = .7.  This means that even if
%there is no frequency in the bin, the synthesized peak will be quite large
%relative to the peak of the bin where the frequency turly does exist.  One
%attempt to fix this is to create more spacing.  (another would be to use
%higher order filters to achieve sharper cutoff).
%I define seperation to be a fraction less than 1.  The 3dB width of the
%bin will be the seperation between 2 frequency bounds times this fraction.
%It will still be centered around the same frequency.

%seperation = 1;

%butterworth order- for a bandpass this must be a multiple of 2
%order = 4;

%c is a vector that lists the relative cutoff of the next closest
%elements.  The current bin will be zeroed if any of the surrounding
%bins are greater than the cutoff.
%ex: c = [1 2 3].  If the surrounding bins are greater than c, the bins
%outside that are greater than 2c, or the bins outside that are greater
%than 3c, this bin will be zeroed.
%This will be make much more complicated when distortion is introduced and
%the spacing of the bins is in no way constant.
%c = [2 3 3];

%real code begins here
%fs = 44100;

%figure out what type of filter to make, not bothering to check the string:


numbins = length(f)-1;

f1 = f(1:(length(f)-1));
f2 = f(2:length(f));
freq = f1+(f2 - f1)/2; %center frequencies needed for synthesis later

%normalize these values so that we can make our filters.  Butter expects
%values between 0 and 1 where 1 corresponds to fs/2.  So when f = fs/2 we
%should get an fnorm of 1.
fnorm = 2* f / fs;

%now make the filter bank
f1norm = fnorm(1:(length(f)-1)); %left norms
f2norm = fnorm(2:length(f)); %right norms
fcnorm = 2*freq/fs; %normalized center frequencies

%these are the new bin 3dB cutoffs
fleft = fcnorm - (f2norm-f1norm) * separation / 2;
fright = fcnorm + (f2norm-f1norm) * separation / 2;

%need to make this generic if we want filters to be a different order
A = zeros(numbins, order + 1);
B = zeros(numbins, order + 1);

type = ftype(1);

if type ~= 'b'
    R = str2double(ftype(2:end));
end

for i = 1:numbins
    if type == 'b'
        [B(i,:), A(i,:)] = butter(order/2, [fleft(i), fright(i)]);
    else
        [B(i,:), A(i,:)] = cheby1(order/2,R, [fleft(i), fright(i)]);
    end
end

n = length(x);
winsize = length(window);
blocksize = (winsize / overlap);

%since this is intended to be a realtime effect, we will 'cut off' the end.
%numblocks is the number of blocks we perform a filter on.
%this discards any block that is less than blocksize, and then discards
%blocks so that we stop when window no longer fits across our data.
numblocks = floor(n / blocksize) - overlap + 1;
y = zeros(blocksize * (numblocks + overlap - 1), 1);

%filters are time averagers.  As such, the wave will need to ramp up from
%0.  To avoid this, we save the state of the filter from each frame and use
%it in the next.  Each state is a vector of length 2;
Zf = cell(numbins, 1);
for k = 1:length(numbins)
    %need a way to make generic in case we change order of filters
    %Zf{k} = [0,0];
    Zf{k} = zeros(1, order);
end

%as we are now pipelining the filtering, we must save the norm of each
%frame to use in future frames.  There will be a norm to save for each
%frequency, and we will need 'overlap' frames to compute a full window.
%This means we will need to store the last 3 blocks.  Keep a counter to
%note which slot should be written to next.

normsqr = zeros(numbins, overlap - 1);
next = 1;

%numblocks
for i = 1:numblocks
    start = (blocksize*(i-1)) + 1;
    stop = (blocksize*( i-1 + overlap ));
    a = x( start : stop ); %current sample to do fft on
    
    %first we need to get a frequency domain representation of our frame.
    %we will do this by using our filter bank to isolate each frequency bin.
    %we take the rms (sqrt(2) * 2-norm / sqrt(# samples)) as a peak value for synthesis later.
    
    peak = zeros(numbins, 1);%peak amplitude of the sine wave to be synthesized
    for j = 1:numbins
        [filt, Zf{j}] = filter(B(j,:),A(j,:),a,Zf{j});
        nsqr = (norm(filt))^2;
        %mix this norm with the norms of previous frames to acquire the
        %proper overlap.
        totalnorm = sqrt(sum(normsqr(j,:)) + nsqr);
        
        %this is our ring buffer of norm^2 values- only used if overlap > 1
        if overlap > 1
            normsqr(j,next) = nsqr;
        end
        
        %compute rms
        peak(j) = sqrt(2) * totalnorm / sqrt(winsize*overlap);%norm of last (overlap) frames
        
    end
    %increment ring buffer index
    next = mod(next,overlap - 1) + 1;
    
    %due to real world filters there will be power observed in bins
    %where there really is no power.  Thus we must set a threshold.  A hard
    %threshold will not work, as softer music will not sound at all then.
    %We use a relative threshold.  We must also take into account the
    %distance between bins.  By the nature of filters, real power from a bin will appear more in the
    %adjacent bin than one 10 bins away.
    
    %cutoff that takes bin distance into account, assuming some sort of equal spacing between bins:
    %recall c is a vector of thresholds
    cb = fliplr(c);
    ctotal = [cb c];
    if ~isempty(ctotal)
    for j = 1:numbins
        %bin to check against:
        %all done in 1 look so we can break as soon as we zero
        comp = [(j-length(c)):(j-1),(j+1):(j+length(c))];
        for k = 1:length(ctotal)
            if comp(k) > 0 && comp(k) <= numbins
                if peak(comp(k)) > comp(k) * peak(j)
                    peak(j) = 0;
                    break;
                end
            end
        end
    end
    end
    
    %now we must resynthesize our signal.  For distorting the signal this
    %will change to have many more waves synthesized.
    
    
    
    %first lets find the times that we are using for this frame
    period = 1/fs;
    tstart = start * period;
    tstop = stop * period;
    t = (tstart:period:tstop)'; %must be a column vector
    
    y(start:stop) = y(start:stop) + window .* synthfn(freq, peak, t);
end