% File: v2sdat2.m
%
% Modified from 
% File:   v2sdat.m
% Date:   15-May-97
% Author: I. Chuang <ike@lanl.gov>
% 
% MATLAB5 FUCNTION: Convert a Varian data file set to a matlab sdat.mat file
% See NEW.DESIGN for the parameter specification.
% 
% 19-May-97 ILC: added multidimensional data handling
% 31-Aug-97 ILC: added avaltab and more multidim array handling
% 09-Oct-98 ILC: fixed avaltab for when no array exists
% 09-Oct-98 ILC: caused baseline correction to complex conj - had wrong hz dir
%                i.e. FREQ AXIS REVERSED with respect to earlier versions
% 05-Jan-99 LV: created parameter "version" to handle expts aborted halfway
%               (see body program)
%
% usage: v2sdat(dirname,phase,apod,nresamp)
% 
% where:
% 
% dirname	- name of the data directory (e.g. 'd8.fid')
% phase  	- phase of real and imaginary fft to use
%	    	  0-360 = normal phase angle (default=calib.phase)
%		  1000  = absolute magnitude
%		  1200  = log10(abs(dat))
%		  1400  = complex phase (in degrees)
%		  2000  = use complex value; do not plot
%		  3000  = plot FID (and return d = time domain data)
%		  may also be a vector [freqp pc0 pc1 ; freqx xc0 xc1 ; ...] in
%		     which case the program chooses the first row that
%		     fits in the frequency range of the current spectrum
%		     and uses it for zeroth and first order phase correction.
% apod		- apodization factor (0.00 = default)
%		  controls how much filtering to use when requested.  0=none.
% nresamp	- oversampling factor; controls size of the fft
%		  1 = fft buffer same as original = TD (default)
%		  n = fft buffer n * original 

function sd = v2sdat2j(fname,receiver,phase,apod,nresamp)

     rcvr='yynn';

if(nargin<5)
  nresamp = 1;
  if(nargin<4)
    apod = evalin('base','calib.gapod');
    if(nargin<3)
      phase = evalin('base','calib.gphase');
      phase(3)
      if(nargin<2)
        receiver=1;
      end
    end
  end
end


if ~ischar(fname)
  if(evalin('base','isfield(calib,''fn'')'))
    prefn = evalin('base','calib.fn');
  else
    prefn = 'exp';
  end
  fname = [prefn num2str(fname)];
end

fname0 = fname;
if ~(exist(fname)==7)
  fname = ['vnmrsys/' fname '/acqfil'];
end
pd = [fname '/procpar'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read parameters from proc

tacq       = vgetparm(pd,'at');
ns         = vgetparm(pd,'nt');
np         = vgetparm(pd,'np');		% # of pts in the fid
sfo        = vgetparm(pd,'sfrq') * 1e6;
swh        = vgetparm(pd,'sw');
[x,pp]     = vgetparm(pd,'pslabel');
[x,dt]     = vgetparm(pd,'date');
[x,aname]  = vgetparm(pd,'array');	% array parameter names

sfrq        = vgetparm(pd,'sfrq') * 1e6;
dfrq        = vgetparm(pd,'dfrq') * 1e6;
dfrq2       = vgetparm(pd,'dfrq2') * 1e6;

ns         = vgetparm(pd,'nt');

[array,nexpt,avaltab] = processarray(pd,aname);
dim = length(array);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine sfo and receiver situation

%sfo=dfrq2;

if (receiver > length(rcvr))
  error('Receiver number requested is higher than number of receivers used.\n')
end
actrcvrs=findstr('y',rcvr)               % e.g. gives [2 3] for rcvr='nyy'
numactrcvrs=length(actrcvrs)          
actrcvrindex=find(actrcvrs==receiver)    % so actrcvrs(actrcvrindex)=receiver

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read in FID
fp = fopen([fname '/fid'],'r','ieee-be');
dheader   = fread(fp,8,'int32');	% read main  header, 32 bytes
nfid = np/2;

%for cnt = 1:nexpt			% each expt's FID stored in a column
for cnt = 1:nexpt*numactrcvrs
  dblockhdr = fread(fp,7,'int32');	% read block header, 28 bytes
  draw      = fread(fp,np,'int32');	% fid data, np bytes

  if(np ~= length(draw))
    fprintf(1,'[v2sdat] oops! short read.  Expected np=%d, got %d\n',...
	np,length(draw));
  end

  if (mod(cnt,numactrcvrs) == mod(actrcvrindex,numactrcvrs)) 
                % if this fid belongs to this sd, then put it in fid
                % both expresssions are =0 or >0
    dd2 = reshape(draw,2,nfid);
    f = dd2(1,:) + i * dd2(2,:);
%   fid(:,cnt) = (f - mean(f))';		     % baseline corrected FID
    fid(:,fix((cnt-1)/numactrcvrs)+1) = (f - mean(f))';  
  end
end

fclose(fp);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apodize and fft to get spectrum

nspect = nfid*nresamp;

apod = matchfrq(sfo,apod,'apod');	% get matching entry for struct
sfo;
apod=0 % ******!!!!!!!!!!!
if(apod~=0)
  filter = exp(-(1:nspect)/nspect*apod)';
else
  filter = 1;
end

for k = 1:nexpt

%fid(1:nspect-10,k)=fid(11:nspect,k); fid(nspect-10:nspect)=0;
%fid=fid-sum(fid)/nspect; 
  spect(:,k) = fftshift(fft(filter.*fid(:,k),nspect));
%  spect(:,k) = fft(filter.*fid(:,k),nspect);
end
fdat = linspace(-swh/2,swh/2,nspect)';
nskip = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% phase spectrum

phase = matchfrq(sfo,phase,'phase')	% get matching entry for struct

for k = 1:nexpt
  spect(:,k) = dophase(spect(:,k),fdat,phase,sfo,swh);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create data structure

sver = 2;

clear sd
sd.sver    = sver;	% structure version
sd.fid     = fid;	% fid complex data
sd.nfid    = nfid;	% number of points in fid
sd.tacq    = tacq;	% acquisition time (time length of the FID)
sd.nskip   = nskip;	% number of points to skip for initial dead time
sd.ns      = ns;	% number of scans performed in acquiring the fid
sd.spect   = spect;	% FFT'ed fid frequency spectrum
sd.nspect  = nspect;	% number of points in the spectrum
sd.fdat    = fdat;	% frequency axis data
sd.sfo     = sfrq;	% spectrometer RF carrier frequency, Hz
sd.dfo     = dfrq;	% spectrometer RF carrier frequency, Hz
sd.swh     = swh;	% spectral bandwidth in hertz
sd.phase   = phase;	% [zeroth, first] order phase correction
sd.apod    = apod;	% amount of exponential apodization
sd.pp      = pp;	% pulse program name
sd.fname   = fname0;	% file format version number
sd.dt      = dt;	% date-time stamp for the experiment
sd.array   = array;	% list of structures giving array'ed parameters
sd.avaltab = avaltab;	% list of array element names and values for each spect
sd.dim     = dim;	% dimension of data (i.e., 0, 1, 2,...)
sd.nexpt   = nexpt;	% number of experiments = size(fid,1), from dim>1
sd.nresamp = nresamp;	% number of times the FID was resampled for the FFT
sd.peaks   = [];
sd.int     = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save data structure to same file but with '.mat' extension

ndot = max(findstr(fname0,'.'));
if(isempty(ndot))
  ndot = length(fname0);
else
  ndot = ndot-1;
end
outfname = [ fname0(1:ndot) '.mat' ];
save(outfname,'sd');

assignin('base','sd',sd);	% update sd structure in base environment

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot

%clf

%dssh

%d = real(spect);
%dmin = min(min(d));
%dmax = max(max(d));
%if(dmin<0)
%  theaxis = [min(fdat) max(fdat) 1.1*dmin 1.1*dmax];
%else
%  theaxis = [min(fdat) max(fdat) 0.9*dmin 1.1*dmax];
%end
%fs = sprintf('Frequency from %f in [Hz]',sfo);

%for cnt = 1:nexpt
%  subplot(nexpt,1,cnt);
%  plot(fdat,d(:,cnt));
%  grid on;
%  axis(theaxis);

%  xlabel(fs);
%  ylabel 'Signal';

%  s = sprintf('[%s] Expt %d, File: %s',date,cnt,[pwd '/' fname]);
%  title(s);

%  fontsize = 8;
%  s = sprintf('ns=%d, pulprog=%s, sw=%6.2f %s=%s',ns,pp,swh,...
%      avaltab{cnt}.name,avaltab{cnt}.val);
%  ylabel2(s,fontsize,1.04+0.015*rem(cnt-1,4));
%end

%imzoom on;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot helper function

function ylabel2(str,size,xpos)

if(nargin<3)
  xpos = 1.05;
  if(nargin<2)
    size = 8;
  end;
end

q = text('posi',[0,0],'vis','off','units','norm','string',str,...
    'color','k','rotation',90,'horiz','center','FontSize',size,...
    'position',[xpos,0.5],'vertical','middle','horiz','center','vis','on');

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% arrayed data helper function

function [array,nexpt,avaltab] = processarray(pd,aname)

global theanames

% get parameter names and number of parameters

if isempty(aname)
  array.name = '';
  array.val = 0;
  nexpt = 1;
  avaltab{1}.name = '';
  avaltab{1}.val = 0;
  return;
end

r = aname;
clear array


version=1   % <<<================= !!!!!!!
%version=2

if version==1

nexpt = 0;
cnt = 0;
while ~isempty(r)
  cnt = cnt + 1;
  [array(cnt).name,r] = strtok(r,' ,');
  [foo,pstr] = vgetparm(pd,array(cnt).name);
  array(cnt).val = str2num(pstr);
  s = sprintf('[%d] name=%s, val=%s',cnt,array(cnt).name,pstr);
  disp(s);
%  nexpt = nexpt + length(array(cnt).val);
end

end
if version==2 
nexpt=4;	% number of arrayed expts
array.name='aux';
array.val=[1:nexpt]
end


% make array entry descriptor strings

theanames = [];
arrayname(array,1,[],[]);
avaltab = theanames;
nexpt = length(avaltab);

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% array entry descriptor string making helper function

function atab = arrayname(array,level,name,val)

global theanames

if(level>length(array))		% base case
  n = length(theanames);
  theanames{n+1}.name = name;
  theanames{n+1}.val = val;
  return;
end

for k = array(level).val
  if(level==1)
    arrayname(array,level+1,[array(level).name],[num2str(k)]);
  else
    arrayname(array,level+1,[name ',' array(level).name],[val ',' num2str(k)]);
  end
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% phase spectrum helper function

function y = dophase(x,fdat,phase,sfo,swh)

if(length(phase)==2)	% do general linear phase correction

  pv = (phase(1) + fdat * phase(2)/swh)*pi/180;
  % 26-Jun-97 ILC: changed above line to include /swh
  y = x .* exp(i * pv);

elseif(size(phase,2)==3)	% general linear phase correction from a set

  for k = 1:size(phase,1)	% go through rows until one is valid
    pent = phase(k,:);
    if(abs(pent(1)-sfo)<2*swh) 
      pv = (pent(2) + fdat * pent(3))*pi/180;
      y = x .* exp(i * pv);
      return;
    end
  end

  s = sprintf('oops! no valid phase found in list = %6.4f',phase);
  disp(s);

elseif(phase <= 360)

  y = exp(i*phase * pi/180) * x;

elseif(phase==1400)

  y = unwrap(angle(x)) * 180/pi;

elseif(phase==1000)

  y = abs(x);

elseif(phase == 1200)

  y = log10(abs(x));

end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function: get structure entry corresponding to frequency

function x = matchfrq(frq,sdat,param)

if ~isstruct(sdat) x = sdat; return; end	% return if not structure

for s = sdat	% search for match
  if (s.frq==frq)
    x = getfield(s,param);
    return; 
  end
end

x = 0;
return;
