function [spgr,cl] = spgr1a_mt1(dat,samp_rate,seg_pwr,seg_inc_fac,opt_str);
% function [spgr,cl] = spgr1a_mt1(dat,samp_rate,seg_pwr,seg_inc_fac,opt_str);
%
% Function to estimate log spectrogram of single time series
% using multitaper windows and periodogram based estimation procedure.
% Includes bias correction and variance calculation.
% 
% Input arguments
%  dat         Input data array.
%  samp_rate   Sampling rate (samples/sec).
%  seg_pwr     Segment length - specified as power of 2.
%  seg_inc_fac Increment factor between segments, 0<seg_inc_fac<=1, 0.5 = 50% overlap.
%  opt_str     Options string (Optional - string input).
%
%  Supported options are: 'r m t'
%     r  Rectification used on each segment.
%     m  Suppress mains in spectral estimate.
%     t  Linear de-trend applied to each segment.
%     M set multitaper window NW value (default is 2.5)
%        Mnw: uses 2*nw tapers, e.g. 'M2', 'M3.5'
%        range for nw is 2 to 5.
%
% Output arguments
%  spgr [freq x segments]: log10 spectrogram estimate
%  cl   Structure with analysis parameters.
%
% function [spgr,cl] = spgr1a_mt1(dat,samp_rate,seg_pwr,seg_inc_fac,opt_str);

% This version 19/3/12, D.M. Halliday.
% Copy of spgr1a_mt, uses pre-generated dpss sequences in ND_dpss_dat

% Check numbers of arguments

if (nargin<4)
  error(' Not enough input arguments');
end  

if (nargout<2)
  error(' Not enough output arguments');
end  

% Checks on input parameters
pts_tot=length(dat);                        % No of data points
seg_size=2^seg_pwr;                         % DFT segment length (T)
seg_inc=fix(seg_size*seg_inc_fac);          % Points in increment between segments
seg_tot=1+fix((pts_tot-seg_size)/seg_inc);  % Number of segments in spectrogram (including overlap)

%Maximum number of tapers seg_pwr limits for Multi-taper analysis
MT_tap_max=11;
MT_seg_pwr_min=6;
MT_seg_pwr_max=14;
NW_default=2.5;

% Process options.
trend_flag=0;  % Set defaults - options off.
mains_flag=0;
rect_flag=0;
MT_flag=1;     % Uses multitaper windows
NW=NW_default;
MT_tap=round(2*NW)-1;   % K ~ 2NW - 1
if (nargin<5)  % No options supplied.
  opt_str='';
end
options=deblank(opt_str);
while (any(options))              % Parse individual options from string.
  [opt,options]=strtok(options);
  optarg=opt(2:length(opt));      % Determine option argument.
  switch (opt(1))
    case 'r'             % Rectification option.        
      rect_flag=1;
    case 't'             % Linear de-trend option.
      trend_flag=1;
    case 'm'             % Mains suppression option.
      mains_flag=1;
    case 'M'             % MT taper option
      NW=str2num(optarg);
      i=round(2*NW)-1;   % K ~ 2NW - 1
	    if (i<1 | i>MT_tap_max)
	      error(['error in option argument -- M',optarg]);
	    end
			MT_tap=i;
    otherwise
      error (['Illegal option -- ',opt]);  % Illegal option.
  end
end

if (trend_flag)            % Index for fitting data with polynomial.
  trend_x=(1:seg_size)';
end

% Load dpss sequences
if (MT_flag)
	global NS_dpss
	if isempty(NS_dpss)
		load NS_dpss_dat
	end
	eval_str=['NS_dpss.pwr',num2str(seg_pwr),'.tap',num2str(MT_tap),';'];
	taper=eval(eval_str);

	eval_str=['NS_dpss.pwr',num2str(seg_pwr),'.eig',num2str(MT_tap),';'];
	MT_eig=eval(eval_str);
end

dat_mt=zeros(seg_size,MT_tap);
u2_norm=1/sqrt(seg_size);

% Storage for individual segments and individual periodograms
dat1=zeros(seg_size,1);
per1=zeros(seg_size,1);
% Starting index for each segment - initialise to first segment.
seg_start=1;

% Storage for Spectrogram - No DC values output, includes fN value
spgr=zeros(seg_size/2,seg_tot);

% Frequency resolution.
deltaf=samp_rate/seg_size;

% Normalising factor for periodogram: 1/(2 pi T)
t_fac=2*pi*seg_size;

% Loop round segments.
for ind=1:seg_tot

	dat1(:,1)=dat(seg_start:seg_start+seg_size-1); % Extract data
	seg_start=seg_start+seg_inc;                   % Update starting index

	dat1=dat1-mean(dat1);                          % Subtract mean.

	if rect_flag
    dat1=abs(dat1);                              % Rectification (Full wave).
  end

	if trend_flag                                  % Linear trend removal.
    p=polyfit(trend_x,dat1,1);                   % Fit 1st order polynomial.
    dat1=dat1-p(1)*trend_x-p(2);                 % Subtract from data.
  end  

	if (MT_flag)
		for ind_MT=1:MT_tap
			dat_mt(:,ind_MT)=dat1;
		end	
		fdat1=fft(dat_mt.*taper)/u2_norm;     % DFT of single segment, using multiple windows
		per1=mean(abs(fdat1.*fdat1)/t_fac,2); % Periodogram, average over multiple windows
	end
	% Multitaper average for each segment.
	spgr(:,ind)=per1(2:seg_size/2+1);

	if mains_flag      % Suppression of mains - smooth out using adjacent values. (May not work well with MT)
	  mains_ind=round(50.0/deltaf);
	  spgr(mains_ind,ind)=0.5*(spgr(mains_ind-2,ind)+sp(mains_ind+2,ind));
	  spgr(mains_ind-1,ind)=0.5*(spgr(mains_ind-2,ind)+sp(mains_ind-3,ind));
	  spgr(mains_ind+1,ind)=0.5*(spgr(mains_ind+2,ind)+sp(mains_ind+3,ind));
	end  
 	
end

% Bias correction. Bias in log10 domain
B=log10(exp(1))*(psi(MT_tap)-log(MT_tap));

% Output log 10 Spectrogram estimate with bias correction
spgr=log10(spgr)-B;

spg_c95=1.96*log10(exp(1))*sqrt(psi(1,MT_tap));
% N.B. Confidence interval for log plot of spectra is TWICE this value.

% Set values in cl structure.
cl.type=30;   % For spectrogram
cl.seg_size=seg_size;
cl.seg_tot=seg_tot;
cl.samp_rate=samp_rate;
cl.seg_inc_fac=seg_inc_fac;
cl.df=deltaf;
cl.freq=(1:seg_size/2)'*deltaf;  % Frequency in Hz.
cl.seg_inc_fac=seg_inc_fac;
cl.seg_inc=seg_inc;
cl.MT_flag=MT_flag;
cl.NW=NW;
cl.MT_tap=MT_tap;
cl.B=B;
cl.spg_c95=spg_c95;

% Set additional elements in cl structure.
cl.opt_str=opt_str;                   % Copy of options string.
cl.what='';														% Field for plot label.

% Display No of segments & resolution 
disp(['Segments: ',num2str(seg_tot),', Segment length: ',num2str(seg_size/samp_rate),' sec, (',num2str(seg_size),' pts), Increment Factor: ',num2str(seg_inc_fac),' (',num2str(seg_inc),' pts), Resolution: ',num2str(deltaf),' Hz.']);
