[W]=dilateColumns(data,varargin)

function [W]=dilateColumns(data,varargin)
%
% [varargout]=dilateColumns(data,varargin)
%
% Dilates a time series input and resizes it to the orignal sample length
% input.  If a matrix is input, dilation is done along columns. % The up/down 
% sampling is done via resample.m, so low-pass filtering of the
% up/downsampling is automatically implemented. An anti-alaising, lowpass 
% filter is applied to the data, and the signal delay is compensated for via 
% resample.
%
% USAGE:
% [W]=dilateColumns(data); %default dilation is 2
% [W]=dilateColumns(data,dilNum); %a vector of dilation values
%
% INPUT:
% data:     A matrix of data vectors stored column-wise.
% dilNum:   A vector of dilation values.  Negative dilation values
%           compress (downsample) the time series, while positive dilation
%           values dilate (upsample) the time series.
%
% OUTPUT:
% W:    The original data with the dilated data.  The output matrix is
%       arranged in order so that dilated data are stored in columns
%       adjacent to the non-dilated input data.
%
%--------------------------------------------------------------------------
% % EXAMPLES
%--------------------------------------------------------------------------
% t             = [linspace(1/200,2000/200,2000)]';
% x             = sin(6*pi*t);
% g             = zeros(size(t));
% g(801:1200)   = [gausswin(400)]';
% data          = [g,g.*x];
%
% %Example 1: Compress and dilate a sine function
% W = dilateColumns(x,[-2,2,4]);
% 
% figure; hold on;
% for k = 1:size(W,2)
%     plot(W(:,k) + 2*k);
% end;
% 
% % %Example 2: Compress and dilate a wave-train
% W = dilateColumns(data,[-2,2,4]);
% 
% figure; hold on; h = zeros(size(W,2),1);
% for k = 1:size(W,2)
%     h(k)=plot(W(:,k) + 2*k);
% end;
% set(h([2,6]),'linewidth',2,'color','r');
%--------------------------------------------------------------------------

%Default value of dilations
dilNum = 2;

for k = 1:length(varargin)

    if(isnumeric(varargin{k}))

        dilNum = varargin{k};

    end;

end;

[M,N]       = size(data);
L           = length(dilNum);
P           = ones(L,1);
Q           = ones(L,1);
Q(dilNum<0) = abs(dilNum(dilNum<0));
P(dilNum>0) = dilNum(dilNum>0);

I           = [1:L+1:N*(L+1)]';
W           = zeros(M,N*L);
W(:,I)      = data;

for k=1:L

    temp = resample(W(:,I),P(k),Q(k));

    if(Q(k)>P(k))

        samp            = 1+M/(2*Q(k)):M-M/(2*Q(k));
        W(samp',I+k)    = temp;

    end;


    if(P(k)>=Q(k))

        i_mid           = ceil(median(1:P(k)*M));
        samp            = i_mid-M/2:i_mid+M/2-1;
        W(:,I+k)        = temp(samp,:);

    end;

end;

NumDown = length(dilNum(dilNum<0));

if(NumDown>0)

    for k=1:NumDown

        if(k==1), U = W; end;

        J           = I + k;
        U(:,I+k-1)  = W(:,J);

    end;

    U(:,I+NumDown)  = W(:,I);
    W               = U;

end;