function [varargout] = projOnClustSpctrgrm(x,W,varargin)
% projOnClustSpctrgrm: A spectrogram type function for computation of
% the energy distribution of a 'long' time series using projOnClusters.m
% It gives a 'time'-distribution of the projection onto the convex cones
% described by W, for x.
%
% The set giving cone k is described by a vector w = W(:,k) and scalar rho:
%
% C_k = { x : x'w >= rho*|| x || };
%
% The polar cone is:
%
% C*_k = { y: y'x <= 0, x element of C_k };
%
% where w is the reference vector for the set C. The projector for this set
% is described in Henry Stark: Vector Space Projections, 1998, pg.109 - 113.
%
%--------------------------------------------------------------------------
% USAGE
%--------------------------------------------------------------------------
%
% [B,varargout] = projOnClustSpctrgrm(x,W,varargin)
%
% ...
%
% [B,nt,indx] = projOnClustSpctrgrm(x,W,rho,overlap,'plot');
%
% INPUT
% x: the signal to be analyzed. Must be a column vector.
% W: A matrix whose columns are the cone representers (see [].pdf for a
% description of these). W must have the same number of rows as x does.
% overlap: The amount of time (%) in overlap each window has. Default is 80.
% A tukeywindow is used to taper each time window. The window
% width is the same as the length of the columns of W.
% rho: A row vector of correlation coefficients corresponding to the columns
% of W.
% plotopt: A string. If plotopt = 'plot', then a spectrogram plot is made.
%
% OUTPUT
% B: The spectrogram matrix
% nt: The time window index. It gives the index of the time vector for
% each window used in computed B.
% indx: The index of the signal used to compute the projection of B from W.
%
% If no output is requested, a plot is made. Alternatively, a plot is
% produced if an input argument 'plot' is given.
%
% % Examples:
% %Example 1: Using 5 signals to obtain a rough estimate of a step function
% N = 1024;
% t = linspace(-2*pi,2*pi,N);
% x = sinc(t);
% temp = zeros(2048,1);
% temp(512:2048-513) = x;
% x = temp;
% W = zeros(64,floor(2*N/64));
% W(2:17,1) = gradient(hann(16));
% for k = 2:floor((2*N)/64)
% temp = sin(k*pi*[1:64]'/64);
% W(:,k) = temp.*W(:,1);
% end;
% [B,windx,clustindx] = projOnClustSpctrgrm(x,W,0.96*ones(1,floor(2*N/64)),50,'plot all');
% disp('A conic-spectrogram decomposition of a truncated sinc function')
% disp('with 50% window overlap and 96% correlation threshold')
% disp('Try plotting the pcolor(log10(B)) as well');
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%% FUNCTION BODY
%
%--------------------------------------------------------------------------
%% Get sizes of input variables
%--------------------------------------------------------------------------
[T,J] = size(W); %The number samples number of convex cones being projected onto
L = size(x,1); %The number of samples of x and the number of time series
%--------------------------------------------------------------------------
%% Default Values for input
%--------------------------------------------------------------------------
overlap = 50; %Convert this to a measure of shifting localization window
overlap = 1-overlap/100;
rho = 0.96*ones(1,J); %Default value for cone geometry
plotopt = 'none'; %Default plot option
%--------------------------------------------------------------------------
%% Check input options
%--------------------------------------------------------------------------
for k = 1:length(varargin)
if(isscalar(varargin{k}))
overlap = varargin{k}; %convert overlap into a shift amount for the localization
overlap = 1-overlap/100; %window used to compute spectrogram
continue;
end;
if(isnumeric(varargin{k}))
rho = varargin{k};
end;
if(ischar(varargin{k}))
plotopt = varargin{k};
end;
end;
%--------------------------------------------------------------------------
%% Zero-pad input if time series length is less than that of cones
%--------------------------------------------------------------------------
if( L < T );
z = zeros(T,1);
z(1:L) = x;
x = z;
x = maxalign(x);
L = size(x,1);
end;
%--------------------------------------------------------------------------
%% Compute required lengths of windows and time series
%--------------------------------------------------------------------------
R = ceil(L/T); %length ratio of time series to cone vectors
z = zeros(R*T,1); %length of zero padded time series
z(1:L) = x; %place time series x at beginning of zero time series
x = z; %replace x with zero padded version
L = size(x,1); %recompute time series length
tau = round(overlap*T); %get time-window length
nW = floor((L+1-T)/tau+1); %Get number of windows used
%--------------------------------------------------------------------------
%% Compute Window-function
%--------------------------------------------------------------------------
g = zeros(L,1); %Make window-function from Tukey Window
g(1:T) = tukeywin(T,0.5); %Use a flat-top tukey-window (parameter controlled)
xg = x.*g; %Multiply window function by time series to compute
%first column of conic-spectrogram
%--------------------------------------------------------------------------
%% Initiliaze all outputs and loop-variables
%--------------------------------------------------------------------------
B = zeros(J,nW); %The spectrogram matrix
nt = repmat(mean(1:tau),1,nW); %The index for the localizing window
indx = [1:J]'; %The index for the cluster waveform/cone
indx = repmat(indx,1,nW); %and it's matrix form
[u,e,b] = projOnClusters(xg(1:T,1),W,rho); %Intiliaze the conic-spectrogram matrix
B(:,1) = b; %from conic projection of localized time series
%--------------------------------------------------------------------------
%% Compute conic-spectrogram columns for all possible window shifts
%--------------------------------------------------------------------------
for k = 2:nW;
g = circshift(g,tau); %Shift windowing function to localize time series
xg = x.*g; %Recompute windowed-time series
I = (k-1)*tau:(k-1)*tau+T-1; %Get indices of localized time series
nt(k) = mean(I); %Get index for spectrogram window
xcut = xg(I); %Cut out windowed time series
Wtemp = xcorAlign([xcut,W]); %Align xcut with cone vectors
W = Wtemp(:,2:end); %Extract aligned cone vectors
xcut = Wtemp(:,1); %Extract xcut after alignment is done
[u,e,b] = projOnClusters(xcut,W,rho); %Compute conic decomposition and output energy from POCS
B(:,k) = b.^2; %Square cone projection norm to get signal energy in that cone
end;
%--------------------------------------------------------------------------
%% Assign requested outputs
%--------------------------------------------------------------------------
if(nargout>0)
for k = 1:nargout
if(nargout>=1), varargout{1}=B; end; %the conic-spectrogram matrix
if(nargout>=2), varargout{2}=nt; end; %the time series index
if(nargout>=3), varargout{3}=indx; end; %the cluster waveform index
end;
end;
%--------------------------------------------------------------------------
%% Plot Spectrogram if 'plot' is input
%--------------------------------------------------------------------------
if(strcmp(plotopt,'plot'))
pcolor(nt,indx,log10(B));
axis xy; axis tight; view(0,90);
shading('interp'); colormap('jet');
set(gcf,'color','w');
end;
%--------------------------------------------------------------------------
%% Plot Exhaustive Summary of Data if 'plot all' is input
%--------------------------------------------------------------------------
if(strcmp(plotopt,'plot all'))
subplot(3,3,[1,4]);
plotColumns(maxalign(W),{'k','r','b'});
set(gca,'Ytick',2:2:indx(end));
xlabel('Sample','fontsize',14,'fontweight','b');
ylabel('Cluster Waveform Index','fontsize',14,'fontweight','b');
set(gca,'fontsize',14,'fontweight','b');
subplot(3,3,[2,3,5,6]);
%pcolor(nt,indx,log10(B));
pcolor(nt,indx,(B));
set(gca,'Yaxislocation','right');
axis xy; axis tight; view(0,90);
set(gca,'Ytick',2:2:indx(end));
shading('interp'); colormap('jet');
xlabel('Window Index','fontsize',14,'fontweight','b');
ylabel('Cluster Waveform Index','fontsize',14,'fontweight','b');
set(gca,'fontsize',14,'fontweight','b');
subplot(3,3,[8,9]);
plot(x); xlim([min(nt),max(nt)]); %axis tight;
xlabel('Sample','fontsize',14,'fontweight','b');
ylabel('Analyzed Signal','fontsize',14,'fontweight','b');
set(gca,'fontsize',14,'fontweight','b');
subplot(3,3,[7]);
s=stem(trapz(B,2),'r'); set(s,'linewidth',2); axis tight;
xlabel('Cluster Index','fontsize',14,'fontweight','b');
ylabel('Energy Sum ','fontsize',14,'fontweight','b');
set(gca,'fontsize',14,'fontweight','b');
set(gcf,'color','w');
end;
