function addNoise(filename, numCopies)
%
% function addNoise(filename, numCopies)
% 
% Adds Noise to 2D trajectories, like this: 
%
% 1. Adding noise by interpolating Gaussian peaks
% 2. Insert 2 time compressions-decompressions, 1 in the 1st half and 
%    one in the 2nd half
%
% ---Input:
% filename: file with 2D trajectory consisting of 2 columns
% numCopies: how many variations of the time-series
%
% ---Output:
% Writes 1_Noise.dat, 2_Noise.dat ... to disk 
% 
%================================================================
% EXAMPLE:
%
% addNoise('1.dat', 50); % creates 50 copies
% 
%================================================================
%
% run without parameters for a demo.....
%
%================================================================
%
% BEGIN COPYRIGHT NOTICE
%
%    addNoise code -- (c) 2002 Michalis Vlachos (http://www.cs.ucr.edu/~mvlachos)
%
%    This code is provided as is, with no guarantees except that 
%    bugs are almost surely present.  Published reports of research 
%    using this code (or a modified version) should cite the 
%    article that describes the algorithm: 
%
%      M. Vlachos, J. Lin, E. Keogh, D. Gunopulos: 
%      "A Wavelet-Based Anytime Algorithm for K-Means Clustering of Time Series", 
%      SDM 2003, San Francisco, Workshop on Clustering High Dimensionality Data 
%      and Its Applications  
%
%    Comments and bug reports are welcome.  Email to mvlachos@cs.ucr.edu 
%    I would also appreciate hearing about how you used this code, 
%    improvements that you have made to it.
%
%    You are free to modify, extend or distribute this code, as long 
%    as this copyright notice is included whole and unchanged.  
%
%    END COPYRIGHT NOTICE


if nargin == 0 
    demo; 
    return;
end 

data = load(filename);

[rows,cols] = size(data);
x=data(:,1); 
y=data(:,2); 


newtrajID = 1;
disp(['Making copies for ', filename]);
   
for i=1:numCopies 
      
    %-------Adding Noise---------
      
    noisyX = x + addSmoothNoise(x);
    noisyY = y + addSmoothNoise(y);

    %---add time shift
    %
    % pick randomly either time compression (0) or decompression (1)
    % from defined percent of traj length
      
    tx = addTimeShift(noisyX);
    ty = addTimeShift(noisyY);
      
    noisyX = []; noisyX = tx;
    noisyY = []; noisyY = ty;
      
    noisyX = noisyX'; % make 1 row
    noisyY = noisyY'; 
      
    %write to file
    newFile = strcat(num2str(newtrajID),'_Noise','.dat')
    fid = fopen(newFile,'w');		
    fprintf(fid,'%8.4f %8.4f\n',[noisyX;noisyY]);
    fclose(fid);
      
    %clear
    noisyX = [];
    noisyY = [];
            
    newtrajID = newtrajID + 1;
end
   
    %plot3(time,x,y,'r-');
    %xlabel('time'); ylabel('x movement');zlabel('y movement');
	%grid on;
	%rotate3d on;
	%title(strcat('Plot of: ',num2str(trajnum)));




%======================================
% Interpolate Gaussian Peaks
%======================================

function smoothed = addSmoothNoise(data)

    n = length(data);
    skip = floor(n/20); % a peak every 20 points
    range = std(data)/5; % change this for more noise

    % noisy peaks
    randStart = round(rand(1,1)*skip);
    noisePos = [randStart:skip:n];
    noise = randn(1,length(noisePos))* range;


    % using random positions leads to great noise...
    % because spline doesn't fit well
    %
    %noisePos = round(unifrnd(skip,n-skip,1,noiseLen-2));
    %noisePos = sort(noisePos);
    %noisePos = [skip, noisePos, n-skip];

    % smooth them
    x = 1:n;

    smoothed = interp1(noisePos,noise,x,'spline');
    smoothed = smoothed';

    % beg. and end may contain a lot of noise
    % so put zero noise
    smoothed([[1:randStart] [noisePos(end) n]]) = 0;



%======================================
% Add time variations
%======================================

function newData = addTimeShift(data)

    % pick randomly either time compression (0) or decompression (1)
  	% from 1..6 percent of traj length
     
    n = length(data)
     
  	% lengths 2..4% of length
    posLen = round(unifrnd(2,4,1,1)); % from 2..4 , 1x1
    removeLen = floor(posLen*n/100)
      
    % pick starting position
    % returns random from [1+removeLen: n-removeLen]
    posStart = round(unifrnd(1+removeLen,n-removeLen-1,1,1))
      
    % 0 or 1
    method = round(rand);
      
    if method==0, % compression==decimation
	    len1_New = posStart - removeLen; 
    else % decompr==interpolation
        len1_New = posStart + removeLen; 
    end
   
    len2_New = n - len1_New

    data1 = resample(data([1:posStart]), len1_New, posStart);
    data2 = resample(data([posStart+1:n]), len2_New, n-posStart);
  
    newData = [data1; data2];
   

%======================================    
% Demo of addNoise 
%======================================

function demo()

    % random walk
    data = randn(100,2)*100;
    data = cumsum(data);
    save('test.dat','data','-ASCII');
    
    addNoise('test.dat', 2);
    
    % plot what happened
    noise = load('1_Noise.dat');
    hold on;
    plot3([1:length(data)], data(:,1), data(:,2), 'k');
    plot3([1:length(data)], noise(:,1), noise(:,2), 'r--');
    legend({'Original' 'Variation'});
    rotate3d on;
    grid on;
    view(-11,38);