function [brukerdata]=brukerimport_SPEN(DataType, inpath)
%   [brukerdata]=brukerimport(DataType, inpath)
%   Imports PFG-NMR data in Bruker format
%               DataType specifies wheter it is raw or Processed data.
%               DataType=0 is raw data, 1 processed, and 2 is processed
%               data for an array
%
%   Useage: Point to the fid/ser file that contains the raw data. The
%           imported data will be returned in the structure brukerdata.
%
%   Example:
%   See also: DOSYToolbox, dosy_mn, score_mn, decra_mn, mcr_mn, varianimport,
%             brukerimport, jeolimport, peakpick_mn, dosyplot_mn,
%             dosyresidual, dosyplot_gui, scoreplot_mn, decraplot_mn,
%             mcrplot_mn
%
%   This is a part of GNAT.
%   Copyright 2017  <Mathias Nilsson>
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%   mathias.nilsson@manchester.ac.uk

%% Changes
% MN      26 June    2017   Handling import and PSYCHEiDOSY experiment,
%                           given correct name for sequence
% MN      17 July    2017   Handling import and PSYCHEiDOSY experiment for
%                           both 2D and 3D versions
% GDP     15 August  2017   Importing 3D corrected for most of our experiments
% MN      23 October 2017   Added errordlg for import error
% MN      09 January 2018   Fixed and error calculating brukerdata.at for
%                           processed data

%% Initialisation
%Check the OS so the paths gets written the right way
if isunix==1
    OSslash='/';
else
    OSslash='\';
end
% FilePath=[FilePath OSslash];
%a=b
if nargin==2
    if DataType==0 %raw data
        %inpath should point to the directory containing the ser/fid file
        FilePath=[inpath OSslash];
        ProcPathRoot=[FilePath 'pdata' OSslash '1' OSslash];
    elseif DataType==1 %processed data
        %inpath should point to the directory containing the 1r/2rr file
        ProcPathRoot=[inpath OSslash];
        StrInd=strfind(inpath,'pdata');
        if isempty(StrInd)
            warndlg('Please select the directory containing the processed data (e.g 1r, 2rr)','Import Bruker Processed Data')
            uiwait(gcf)
        end
        
        FilePath=inpath(1:(StrInd-1));
        
    else
        error('Unknown DataType')
    end
    % [file FilePath]=uigetfile('*','Choose the Bruker data file containing the FID (*.ser or *.fid)');
else
    error('brukerimport takes exactly 2 arguments')
    %FilePath=uigetdir(pwd,'Choose the directory containing the NMR experiment (ser or fid file)');
end
%% Main
if FilePath
    %% set some default parameters
    brukerdata.sp1=[];
    brukerdata.ni=[];
    brukerdata.logfile={'brukerimport'};
    %% Check data dimension
    if exist([FilePath OSslash 'acqu4s'], 'file') == 2
        error('Data of order >3 is not supported')
    elseif exist([FilePath OSslash 'acqu3s'], 'file') == 2
        DataDim=3;
    elseif exist([FilePath OSslash 'acqu2s'], 'file') == 2
        DataDim=2;
    elseif exist([FilePath OSslash 'acqus'], 'file') == 2
        DataDim=1;
    else
        errordlg( {'Can not find the acqus file','It is normally found in the numbered data directory'},'File error')
        uiwait;
        brukerdata.acqus_error=1;
        return
    end
    %% Read in acquisition parameters
    if DataDim>=1
        brukerdata.acqus = GetParameters([FilePath  'acqus']);
        brukerdata.acqus.PULPROG(1)=[];
        brukerdata.acqus.PULPROG(end)=[];
        brukerdata.pulprog=brukerdata.acqus.PULPROG;
        logtext=['pulse programme is: ' brukerdata.pulprog];
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        
        % % % % % Added by Rituraj % % % % %
%         Most of these parameters are required for processing spen2GNAT 
%         function and also used by GNAT function for processing the fid  
        if contains(brukerdata.pulprog,'dosyuf')... 
                || contains(brukerdata.pulprog,'ufdosy')
            brukerdata.shape=readacqus([FilePath  'spnam1']);
            brukerdata.sf=str2double(brukerdata.acqus.BF1);
            brukerdata.sfo1=str2double(brukerdata.acqus.SFO1);
            brukerdata.N2 = brukerdata.acqus.L(4); % number of loops in the acquisition
            brukerdata.swtemp = str2double(brukerdata.acqus.SW_h);% spectral width
            brukerdata.NP = str2double(brukerdata.acqus.TD)/2; % total number of points in one fid
            brukerdata.dw = 1/(brukerdata.swtemp); % dwell time
            brukerdata.offsetconv=brukerdata.sfo1/brukerdata.sf-1; % offset
            brukerdata.D6=brukerdata.acqus.D(7); % acqusition gradient duration
            brukerdata.td=brukerdata.acqus.TD; % time domaine size
            brukerdata.d11=brukerdata.acqus.D(12);
            brukerdata.Ta= ((brukerdata.NP*brukerdata.dw)/(2*brukerdata.N2)); %Acquisiton Time
            brukerdata.sw=1/(2*brukerdata.Ta)/brukerdata.sf;
            brukerdata.expO1=str2double(brukerdata.acqus.O1);
            brukerdata.Filepath=FilePath;
        else
            brukerdata.sw=str2double(brukerdata.acqus.SW);
        end
        % % % % % end % % % % %
    end
    
    if DataDim>=2
        brukerdata.acqu2s = GetParameters([FilePath  'acqu2s']);
    end
    
    if DataDim>=3
        brukerdata.acqu3s = GetParameters([FilePath  'acqu3s']);
    end
    %% Read Processing parameters
    %get processing parameters, if there are any
    if DataDim>=1
        Procpath=[ProcPathRoot 'procs'];
        if exist(Procpath, 'file') == 2
            % % % % % Added by Rituraj % % % % %
            if contains(brukerdata.pulprog,'dosyuf')... 
                || contains(brukerdata.pulprog,'ufdosy')
                brukerdata.procs = GetParameters(Procpath);
                brukerdata.sp=(brukerdata.sfo1/brukerdata.sf-1)*1e6+0.5*brukerdata.sw*(brukerdata.sfo1/brukerdata.sf)-brukerdata.sw;
                brukerdata.lp=0;
                brukerdata.rp=0;
                brukerdata.isuf=1;
                % % % % % end % % % % %
            else
                brukerdata.procs = GetParameters(Procpath);
                brukerdata.sf=str2double(brukerdata.acqus.BF1);
                brukerdata.sfo1=str2double(brukerdata.acqus.SFO1);
                brukerdata.sp=(brukerdata.sfo1/brukerdata.sf-1)*1e6+0.5*brukerdata.sw*(brukerdata.sfo1/brukerdata.sf)-brukerdata.sw;
                brukerdata.lp=str2double(brukerdata.procs.PHC1);
                brukerdata.rp=str2double(brukerdata.procs.PHC0)*-1 - brukerdata.lp;
            end
            
        else
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='estimating processing parameters from acquisition parameters';
            brukerdata.logfile{end +1}=logtext;
            brukerdata.sf=str2double(brukerdata.acqus.BF1);
            brukerdata.sfo1=str2double(brukerdata.acqus.SFO1);
            brukerdata.sp=(brukerdata.sfo1/brukerdata.sf-1)*1e6+0.5*brukerdata.sw*(brukerdata.sfo1/brukerdata.sf)-brukerdata.sw;
            brukerdata.lp=0;
            brukerdata.rp=0;
            brukerdata.procs=[];
        end
    end
    
    if DataDim>=2
        Procpath=[ProcPathRoot 'proc2s'];
        if exist(Procpath, 'file') == 2
            brukerdata.proc2s = GetParameters(Procpath);
            
        else
            
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='No processing paramaters for the second dimension';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.proc2s=[];
        end
    end
    
    if DataDim>=3
        Procpath=[ProcPathRoot 'proc3s'];
        if exist(Procpath, 'file') == 2
            brukerdata.proc3s = GetParameters(Procpath);
            
        else
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='No processing paramaters for the second dimension';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.proc3s=[];
        end
    end
    %% Sort out parameters in DOSY Toolbox format
    
    if DataDim==1
        logtext='1D data';
        brukerdata.arraydim=1;
        brukerdata.ngrad=1;
        fidpath=[FilePath 'fid'];
        brukerdata.acqu2s=[];
        brukerdata.acqu3s=[];
    elseif DataDim==2
        logtext='2D data';
        brukerdata.acqu3s=[];
        % % % % % % Added bt Rituraj % % % % % %
        if contains(brukerdata.pulprog,'dosyuf')... 
                || contains(brukerdata.pulprog,'ufdosy')
            brukerdata.numscan=1;
            brukerdata.np=str2double(brukerdata.acqus.TD)/2;
            brukerdata.arraydim=str2num(brukerdata.acqu2s.TD);
        else
            brukerdata.arraydim=str2num(brukerdata.acqu2s.TD); %#ok<*ST2NM>
            brukerdata.np=str2double(brukerdata.acqus.TD)/2; % its place has been changed by Rituraj
        end
        fidpath=[FilePath 'ser'];
        brukerdata.sw1=str2double(brukerdata.acqu2s.SW);
    elseif DataDim==3
        logtext='3D data';
        brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
        fidpath=[FilePath 'ser'];
        brukerdata.sw1=str2double(brukerdata.acqu3s.SW);
    else
        error('Unknown data dimension')
    end
    disp(logtext)
    brukerdata.logfile{end +1}=logtext;
    %Get the parameters I need for now (in units I like)
    brukerdata.filename=FilePath(1:end-1);
    %     brukerdata.np=str2double(brukerdata.acqus.TD)/2;
    brukerdata.sfrq=str2double(brukerdata.acqus.SFO1);
    brukerdata.at=brukerdata.np/(brukerdata.sw*brukerdata.sfrq);
    
    if DataType==1
        if isempty(brukerdata.procs)
            error('cannot import processed data if the file "procs" is missing')
        end
        brukerdata.SI=round(str2double(brukerdata.procs.SI));
        brukerdata.at=0.5*brukerdata.at/(brukerdata.np/brukerdata.SI);
    end
    
    
    switch brukerdata.acqus.NUC1
        case '<1H>'
            brukerdata.gamma=267524618.573;
        case '<2H>'
            brukerdata.gamma=41066000;
        case '<10B>'
            brukerdata.gamma=28747000;
        case '<11B>'
            brukerdata.gamma=85847000;
        case '<13C>'
            brukerdata.gamma=67283000;
        case '<14N>'
            brukerdata.gamma=19338000;
        case '<15N>'
            brukerdata.gamma=-27126000;
        case '<17O>'
            brukerdata.gamma=-36281000;
        case '<19F>'
            brukerdata.gamma=251815000;
        case '<23Na>'
            brukerdata.gamma=70808000;
        case '<27Al>'
            brukerdata.gamma=69764000;
        case '<29Si>'
            brukerdata.gamma=-53190000;
        case '<31P>'
            brukerdata.gamma=108394000;
            
        otherwise
            logtext='unknown nucleus - defaulting to proton';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.gamma=267524618.573;
    end
    
    
    %% Import lists
    
    % Checking for a NATimport list
    if exist([FilePath 'NATimport'],'file') == 2
        % will use NATimport file to determine how to import the data
        logtext='Using NATimport file to determine import parameters';
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        fileid=fopen([FilePath  'NATimport'],'rt'); %#ok<NASGU>
    else
        %Using default import options
        logtext=['cannot open file: ' [FilePath  'NATimport']];
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        logtext='Using default import (normal for 1D, DOSY and Relaxation data';
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        if DataDim==1
            logtext='1D experiment';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.dosyconstant='non existing';
            brukerdata.ngrad=1;
            brukerdata.arraydim=1;
        elseif DataDim == 2 %likely a difusion or relaxation experiment
            logtext='2D experiment: trying diffusion and relaxation type';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.dosyconstant='non existing';
            brukerdata.ngrad=1;
            TD= str2num(brukerdata.acqu2s.TD);
            if exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment
                logtext='Diffusion experiment (difflist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
                brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                brukerdata.ngrad=length(brukerdata.Gzlvl);
                logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                %check for pulseprogram
                if strncmpi('Doneshot',brukerdata.pulprog,8)
                    logtext='Doneshot sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                    brukerdata.tau=brukerdata.acqus.CNST(18);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                elseif strfind(brukerdata.pulprog,'psyche_idosy')
                    logtext='monopolar sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                elseif strfind(brukerdata.pulprog,'oneshot')
                    logtext='oneshot sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                    brukerdata.tau=brukerdata.acqus.CNST(18);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                elseif strncmpi('stebp',brukerdata.pulprog,5) || strncmpi('ledbp',brukerdata.pulprog,5)
                    %bipolar sequence
                    logtext='bipolar sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21)-brukerdata.acqus.D(17)/2;
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(31)*1e-6;
                    brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*(brukerdata.DELTA-brukerdata.delta/3-brukerdata.tau/2);
                elseif strncmpi('ste',brukerdata.pulprog,3) || strncmpi('led',brukerdata.pulprog,3)
                    logtext='monopolar sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                else
                    logtext='unknown sequence; setting parameters to bipolar';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6; %probably assuming bipolar pulse
                    brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*(brukerdata.DELTA-brukerdata.delta/3);
                end
            elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                logtext='T2 experiment (vclist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.vclist = ImportList([FilePath 'vclist'], TD);
                if exist([FilePath 'vdlist'],'file') == 2
                    brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD);
                    logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                else
                    logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                    uiwait(gcf)
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                end
            elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                logtext='T1 experiment (vdlist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD);
                logtext='using vdlist directly for delay values';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                % % % % % % % % % % % % % % % % Added by Rituraj % % % % % % % % % % % %
            elseif contains(brukerdata.pulprog,'dosyuf')... 
                || contains(brukerdata.pulprog,'ufdosy')
                logtext='ultrafast dosy';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.DELTA=brukerdata.acqus.D(31);
                brukerdata.delta=(brukerdata.shape.SHAPE_LENGTH(1)*1e-6+brukerdata.acqus.D(12))/2;
                brukerdata.dosyconstant=brukerdata.gamma^2*brukerdata.delta^2*(brukerdata.DELTA-brukerdata.delta);
            else
                logtext='No recognizable files with increment values exist: importing empty arrays';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        else %higher order experiments, GDP 17 Jan 2018
            logtext='higher order experiment. try to sort experimental parameters'; disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
            brukerdata.ngrad=str2num(brukerdata.acqu2s.TD);
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.dosyconstant='non existing';
            if exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment
                logtext='3D Diffusion experiment (difflist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                if  strfind(brukerdata.pulprog,'psyche_idosy')
                    logtext='psyche_idosy sequence; setting paramaters to monopolar';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD= str2num(brukerdata.acqu2s.TD);
                    brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                    brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
                elseif strfind(brukerdata.pulprog,'oneshot')
                    logtext='oneshot sequence';
                    disp(logtext)
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                    brukerdata.tau=brukerdata.acqus.CNST(18);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                else
                    logtext='unknown sequence; setting parameters to bipolar';
                    disp(logtext)
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl2(isnan(Gzlvl2))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6; %probably assuming bipolar pulse
                    brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*(brukerdata.DELTA-brukerdata.delta/3);
                end
                if exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                    logtext='T2 experiment (vclist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vclist=ImportList([FilePath 'vclist'], TD2);
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(3) + brukerdata.acqus.P(3)*1e-6);
                    if exist([FilePath 'vdlist'],'file') == 2
                        brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                        logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                    else
                        logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                        warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                        uiwait(gcf)
                    end
                elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                    logtext='T1 experiment (vdlist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                else
                    logtext='only diffusion list found, I will sort like this';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                end
            elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                logtext='T2 experiment (vclist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                TD3= str2num(brukerdata.acqu3s.TD);
                TD2= str2num(brukerdata.acqu2s.TD);
                vclist3=ImportList([FilePath 'vclist'],TD3);
                vclist3(isnan(vclist3))=0;
                vclist2=ImportList([FilePath 'vclist'],TD2);
                vclist2(isnan(vclist2))=0;
                if nnz(vclist3)>nnz(vclist2)
                    brukerdata.vclist=vclist3;
                elseif nnz(vclist3)<nnz(vclist2)
                    brukerdata.vclist=vclist2;
                elseif nnz(vclist3)==nnz(vclist2)
                    if (length(vclist3)-nnz(vclist3))==0
                        brukerdata.vclist=vclist3;
                    else
                        brukerdata.vclist=vclist2;
                    end
                end
                brukerdata.ngrad=length(brukerdata.vclist);
                if exist([FilePath 'vdlist'],'file') == 2
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD2);
                    logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                else
                    %TD= str2num(brukerdata.acqu2s.TD);
                    %brukerdata.vclist = ImportList([FilePath 'vclist'], TD);
                    logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                    uiwait(gcf)
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                end
            elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                logtext='T1 experiment (vdlist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                TD3= str2num(brukerdata.acqu3s.TD);
                TD2= str2num(brukerdata.acqu2s.TD);
                vdlist3=ImportList([FilePath 'vdlist'],TD3);
                vdlist3(isnan(vdlist3))=0;
                vdlist2=ImportList([FilePath 'vdlist'],TD2);
                vdlist2(isnan(vdlist2))=0;
                if nnz(vdlist3)>nnz(vdlist2)
                    brukerdata.vdlist=vdlist3;
                elseif nnz(vdlist3)<nnz(vdlist2)
                    brukerdata.vdlist=vdlist2;
                elseif nnz(vdlist3)==nnz(vdlist2)
                    if (length(vdlist3)-nnz(vdlist3))==0
                        brukerdata.vdlist=vdlist3;
                    else
                        brukerdata.vdlist=vdlist2;
                    end
                end
                brukerdata.ngrad=length(brukerdata.vdlist);
                logtext='using vdlist directly for delay values';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            else
                logtext='Any list was found';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        end
    end
    
    %% Read in the binary data
    %check the format of the data
    %The byte order will be big from SGI (MIPS processor) and little from
    %Linux/Windows (Intel), this is stored as the status parameters BYTORDA
    %(acquisition data) and BYTORDP (processed data).  Processed data is always
    %stored as scaled 32 bit integers, but acquisition data can in principle be
    %present as double precision (64 bit) data - this is indicated by the
    %status parameter DTYPA (0=integer, 2=double).  The double state is
    %automatically triggered if the raw data would overflow (but even in the
    %modern age of large effective ADC resolution this is probably a rare event!)
    %a=brukerdata.acqus.BYTORDA
    % a=brukerdata.acqus.BYTORDP
    % a=brukerdata.acqus.DTYPA
    %a=brukerdata.acqus.DTYPP
    % brukerdata.acqus.BYTORDA='0';
    switch str2double(brukerdata.acqus.BYTORDA)
        case 0 %little endian
            byte_format='l';
        case 1  %big endian
            byte_format='b';
        otherwise
            error('unknown data format (BYTORDA)')
    end
    switch str2double(brukerdata.acqus.DTYPA)
        case 0 %32 bit integer
            byte_size='int32';
        case 1  %double
            byte_size='double';
        case 2  %double
            byte_size='double'; %MN 13March2017 Not sure if 1 or 2 does any difference to Matlab
        otherwise
            error('unknown data format (DTYPA)')
    end
    
    %Check that np is a multiple of 256 - otherwise correct
    % In a ser file each 1D fid will start at a 1024 byte block boundary even
    % if its size is not a multiple of 1k (256 points)
    corrpoints=rem(brukerdata.np*2,256);
    if corrpoints>0
        corrpoints=256-corrpoints;
        brukerdata.np=brukerdata.np+corrpoints/2;
    end
    
    FID=zeros(brukerdata.np,brukerdata.arraydim);
    
    
    %if importing the processed spectrum
    if DataType==1 %Processed data
        byte_format_proc='l';
        %Seems like it should always be l
        if DataDim==1
            ProcFIDpath=[ProcPathRoot '1r'];
            if exist(ProcFIDpath, 'file') == 2
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                realfid=fread(fileid_ProcFIDreal,brukerdata.SI*brukerdata.arraydim,byte_size);
                realfid=realfid*2^str2num(brukerdata.procs.NC_proc);
                FID=ifft(ifftshift(realfid));
                FID=FID(1:round(brukerdata.SI/2));
                fclose(fileid_ProcFIDreal);
            else
                error(['Processed data 1D data 1r does not exist in:' ProcFIDpath])
            end
            
        elseif DataDim==2
            ProcFIDpath=[ProcPathRoot '2rr'];
            
            if exist(ProcFIDpath, 'file') == 2
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                %  ProcFIDpath
                SI1=str2num(brukerdata.procs.SI);
                SI2=str2num(brukerdata.proc2s.SI);
                brukerdata.np=SI1;
                XDIM1=str2num(brukerdata.procs.XDIM);
                XDIM2=str2num(brukerdata.proc2s.XDIM);
                
                realfid=fread(fileid_ProcFIDreal,brukerdata.np*2*SI2,byte_size);
                realfid=realfid*2^str2num(brukerdata.proc2s.NC_proc);
                
                NoSM = SI1*SI2/(XDIM1*XDIM2);
                NoSM2 = SI2/XDIM2;
                
                %recipe from Nils Nyberg - works most of the time. Should
                %make it work 100% of the time when I feel like it.
                realfid = reshape(...
                    permute(...
                    reshape(...
                    permute(...
                    reshape(realfid,XDIM1,XDIM2,NoSM),...
                    [2 1 3]),...
                    XDIM2,SI1,NoSM2),...
                    [2 1 3]),...
                    SI1,SI2)';
                
                realfid=realfid';
                
                %This saves the processed 2D data for use with e.g. DOSY3D
                brukerdata.X2DSPEC=realfid;
                
                
                FID=realfid(:,1:brukerdata.arraydim);
                
                
                for p=1:brukerdata.arraydim
                    FID(:,p)=ifft(ifftshift(FID(:,p)));
                end
            else
                error(['Processed data 2D data 2rr does not exist in:' ProcFIDpath])
            end
        else
            error('Import of processed data for more that two dimensions is not implemented')
        end
        
        brukerdata.FID=FID(1:brukerdata.np/2,:);
        brukerdata.lp=0;
        brukerdata.rp=0;
        % % % % % % % Added by Rituraj % % % % % % %
    elseif contains(brukerdata.pulprog,'dosyuf')... 
                || contains(brukerdata.pulprog,'ufdosy')
        % to obtain input parameters and brukerdata structure specified in 
        % inputpspendosy_GNAT function and also retain the parameters,
        % specified in brukerdata structure (such parameters are expO1), 
        % from brukerimport_SPEN 
        [inputp,brukerdata]=inputpspendosy_GNAT(brukerdata);
        % to convert spen data into standard NMR data
        [~,brukerdata,~]=spen2GNAT...
            (FilePath,inputp.NBC_spen,inputp.refprof_y_n,brukerdata...
            ,inputp.acqref,inputp.chirpref,inputp.add_sm_proc);
    else %Raw FID data
        
        fileid_FID=fopen(fidpath,'r',byte_format);
        
        impfid=fread(fileid_FID,brukerdata.np*2*brukerdata.arraydim,byte_size);
        compfid=complex(impfid(1:2:end),-impfid(2:2:end));
        [m,~]=size(compfid);
        
        if m<brukerdata.np
            logtext='Actual points in FID and number of points resported (TD) are different. This can hapen after the use of the "convdta" macro. I will use the actual number and continue import';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.np=m;
            FID=zeros(brukerdata.np,brukerdata.arraydim);
        end
        
        for k=1:brukerdata.arraydim
            FID(:,k)=compfid((k-1)*brukerdata.np+1:k*brukerdata.np);
        end
        
        brukerdata.np=brukerdata.np-corrpoints/2;
        brukerdata.FID=FID((1:brukerdata.np),:);
        fclose(fileid_FID);
    end
    
    %% Sort out the digital filtering
    % This is needed for converion of digitally filtered data - se longer
    % explanation at the end of the file
    BrukDigital=[
        2       44.750       46.000       46.311;
        3       33.500       36.500       36.530;
        4       66.625       48.000       47.870;
        6       59.083       50.167       50.229;
        8       68.563       53.250       53.289;
        12       60.375       69.500       69.551;
        16       69.531       72.250       71.600;
        24       61.021       70.167       70.184;
        32       70.016       72.750       72.138;
        48       61.344       70.500       70.528;
        64       70.258       73.000       72.348;
        96       61.505       70.667       70.700;
        128       70.379       72.500       72.524;
        192       61.586       71.333            0;
        256       70.439       72.250            0;
        384       61.626       71.667            0;
        512       70.470       72.125            0;
        768       61.647       71.833            0;
        1024       70.485       72.063            0;
        1536       61.657       71.917            0;
        2048       70.492       72.031            0];
    
    % first check if GRPDLY exists and use that if so
    decim=1;
    dspfvs=1;
    if isfield(brukerdata.acqus,'GRPDLY') && str2double(brukerdata.acqus.GRPDLY)~=-1
        % % % % % % % % % % Added by Rituraj % % % % % % % % % %
        if contains(brukerdata.pulprog,'dosyuf')... 
                || contains(brukerdata.pulprog,'ufdosy')
            brukerdata.digshift=0;
        else
            % % % % % % % % % % end % % % % % % % % % %
            brukerdata.digshift=str2double(brukerdata.acqus.GRPDLY);
        end
        %disp('Hi')
    else
        %brukerdata.acqus.DECIM=0;
        
        if isfield(brukerdata.acqus,'DECIM')
            decim=str2double(brukerdata.acqus.DECIM);
            switch decim
                case 2
                    decimrow=1;
                case 3
                    decimrow=2;
                case 4
                    decimrow=3;
                case 6
                    decimrow=4;
                case 8
                    decimrow=5;
                case 12
                    decimrow=6;
                case 16
                    decimrow=7;
                case 24
                    decimrow=8;
                case 32
                    decimrow=9;
                case 48
                    decimrow=10;
                case 64
                    decimrow=11;
                case 96
                    decimrow=12;
                case 128
                    decimrow=13;
                case 192
                    decimrow=14;
                case 256
                    decimrow=15;
                case 384
                    decimrow=16;
                case 512
                    decimrow=17;
                case 768
                    decimrow=18;
                case 1024
                    decimrow=19;
                case 1536
                    decimrow=20;
                case 2048
                    decimrow=21;
                otherwise
                    logtext='unknown value of DECIM parameter in acqus - cannot set compensation for digital filtering';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    decim=0;
                    decimrow=Inf;
            end
        else
            logtext='no DECIM parameter in acqus - cannot set compensation for digital filtering';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            decim=0;
            decimrow=Inf;
        end
        if isfield(brukerdata.acqus,'DSPFVS')
            dspfvs=str2double(brukerdata.acqus.DSPFVS);
            switch dspfvs
                case 10
                    dspfvsrow=2;
                case 11
                    dspfvsrow=3;
                case 12
                    dspfvsrow=4;
                otherwise
                    logtext='unknown value of DSPVFS parameter in acqus - cannot set compensation for digital filtering';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    dspfvs=0;
                    dspfvsrow=0;
            end
        else
            logtext='no DSPFVS parameter in acqus - cannot set compensation for digital filtering';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            dspfvs=0;
        end
        if (decimrow>14) && (dspfvsrow==4)
            logtext='unknown combination of DSPVFS and DECIM parameters in acqus - cannot set compensation for digital filtering';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            dspfvs=0;
            decim=0;
        end
        
    end
    if (decim==0) || (dspfvs==0)
        %No digital filtering
        logtext='Parameters for digital filtering unknown - assumed to be data without digital filtering';
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        brukerdata.digshift=0;
    elseif (decim==1) && (dspfvs==1)
        %digital filtering set by GRPDLY
        %do nothing
    else
        brukerdata.digshift=BrukDigital(decimrow,dspfvsrow);
    end
    %brukerdata.digshift=0
    
    if DataType==1
        brukerdata.digshift=0;
    end
    
    brukerdata.digshift=round(brukerdata.digshift);
    
    
    brukerdata.np=length(brukerdata.FID(:,1));
    %% Check the acquisition order for 3D Data
    if DataDim==3
        brukerdata.procs;
        brukerdata.acqus;
        brukerdata.acqus.CNST(1);
        brukerdata.acqus.CNST(2);
        brukerdata.acqus.CNST(3);
        brukerdata.acqus.CNST(4);
        brukerdata.acqus.CNST(5);
        if isfield(brukerdata.acqus,'AQSEQ')
            if str2double(brukerdata.acqus.AQSEQ)==0
                %all is fine - this is the structure we expect
                logtext='Normal aquisition order';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            elseif str2double(brukerdata.acqus.AQSEQ)==1
                % reversed order - lets fix it
                logtext='Reversed aquisition order';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                orgfid=brukerdata.FID;
                
                % brukerdata.FID=brukerdata.FID.*0;
                
                array2=brukerdata.arraydim/brukerdata.ngrad;
                
                for k=1:array2
                    for p=1:brukerdata.ngrad
                        % brukerdata.FID(:,k+array2*(p-1))=orgfid(:,(k-1)*brukerdata.ngrad +p);
                        brukerdata.FID(:,(k-1)*brukerdata.ngrad +p)=orgfid(:,k+array2*(p-1));
                    end
                end
            else
                %we're not sure
                logtext='Unknown aquisition order';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        else
            logtext='no AQSEQ parameter. Cannot determine aquisition mode. Sticking with default';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
        end
    else
        %not a 3D set
    end
    
else
    brukerdata=[];
end
end
%% Utility Functions
function OutParams=GetParameters(FilePath)
%acquspath=[FilePath '/acqus'];
%Reading in all the acquisition or processing parameters from the acquXs or
% procXs files
fileid=fopen(FilePath,'rt');
if fileid==-1
    error(['cannot open file: ' acquspath])
end
kk=1;
parmflag=0;
while 1
    if parmflag==1
        parmflag=0;
    else
        parmline=fgetl(fileid);
    end
    if parmline==-1;  break;  end
    
    if length(parmline)<2
        logtext='probably a cut of line - I will continue anyway';
        disp(logtext)
    elseif strcmp(parmline(1:2),'##')
        if strcmp(parmline(3),'$')
            parmend=find(parmline=='=');
            if strcmp(parmline(parmend+2),'(')
                parmname=parmline(4:(parmend(1))-1);
                parmlist=[];
                while 1
                    parmline=fgetl(fileid);
                    if strcmp(parmline(1),'#') || strcmp(parmline(1),'<') || strcmp(parmline(1),'>') || strcmp(parmline(1),'$')
                        parmflag=1;
                        break
                    else
                        parmcell=textscan(parmline,'%q',1024);
                        tadam=parmcell{1};
                        tmplist=zeros(length(tadam),1);
                        for m=1:length(tadam)
                            tmplist(m)=str2num(cell2mat(tadam(m)));
                        end
                        parmlist=[parmlist; tmplist];                   %#ok<AGROW>
                    end
                end
                procpar.(parmname)=parmlist;
            else
                procpar.(parmline(4:(parmend(1))-1))=parmline((parmend+2):end);
            end
        else
            parmend=find(parmline=='=');
            procpar.(parmline(3:(parmend(1)-1)))=parmline((parmend+2):end);
        end
    elseif strcmp(parmline(1:2),'$$')
        str=['acqinfo' num2str(kk)];
        procpar.(str)=parmline(4:end);
        kk=kk+1;
    else
        %do nothing
    end
end
OutParams=procpar;
fclose(fileid);

end
function Vals = ImportList(filepath, TD)
% import gradient values from difflist
fileid_gzlvl=fopen(filepath,'rt');
if (fileid_gzlvl~=-1)
    Vals=zeros(1,TD);
    for k=1:TD
        parmline=fgetl(fileid_gzlvl);
        Vals(k)=str2double(parmline);
        if parmline==-1
            logtext='Warning! -number of  values in list does not correspond to the number ponits in TD';
            disp(logtext)
            break;
        end
    end
    
    fclose(fileid_gzlvl);
else
    logtext='no list file; cannot read invalues';
    disp(logtext)
    %brukerdata.logfile{end +1}=logtext;
end
end
%% INFORMATION ABOUT BRUKER DIGITAL FILTER
% For older versions of the Bruker hardware:
%
% A nice piece was found on the internet on how to calculate the number of points
% semi-automatically. Note that currently matNMR doesn't allow for the necessary
% negative-time apodization.
%
%
%    W. M. Westler and F.  Abildgaard
%    July 16, 1996
%
%    The introduction of digital signal processing by Bruker in their DMX
%    consoles also introduced an unusual feature associated with the data. The
%    stored FID no longer starts at its maximum followed by a decay, but is
%    prepended with an increasing signal that starts from zero at the
%    first data point and rises to a maximum after several tens of data points.
%    On transferring this data to a non-Bruker processing program such as FELIX,
%    which is used at NMRFAM, the Fourier transform leads to an unusable spectrum
%    filled with wiggles. Processing the same data with Bruker's Uxnmr
%    program yields a correct spectrum. Bruker has been rather reluctant
%    to describe what tricks are implemented during their processing protocol.
%
%    They suggest the data should be first transformed in Uxnmr and then inverse
%    transformed, along with a GENFID command, before shipping the data to another
%    program. Bruker now supplies a piece of software to convert the digitally
%    filtered data to the equivalent analog form.
%    We find it unfortunate that the vendor has decided to complicate
%    the simple task of Fourier transformation. We find that the procedure
%    suggested by Bruker is cumbersome, and more so, we are irritated since
%    we are forced to use data that has been treated with an unknown procedure.
%    Since we do not know any details of Bruker's digital filtration procedure
%    or the "magic" conversion routine that is used in Uxnmr, we have been forced
%    into observation and speculation. We have found a very simple, empirical
%    procedure that leads to spectra processed in FELIX that are identical,
%    within the noise limits, to spectra processed with Uxnmr. We deposit
%    this information here in the hope that it can be of some
%    use to the general community.
%    The application of a nonrecursive (or recursive) digital filter to time
%    domain data is accomplished by performing a weighted running average of
%    nearby data points. A problem is encountered at the beginning of
%    the data where, due to causality, there are no prior values. The
%    weighted average of the first few points, therefore, must include data
%    from "negative" time. One naive procedure, probably appropriate to NMR
%    data, is to supply values for negative time points is to pad the data with
%    zeros. Adding zeros (or any other data values) to the beginning of
%    the FID, however, shifts the beginning of the time domain data (FID) to
%    a later positive time. It is well known that a shift in the time
%    domain data is equivalent to the application of a frequency-dependent,
%    linear phase shift in the frequency domain. The 1st order phase shift
%    corresponding to a time shift of a single complex dwell point is 360 degrees
%    across the spectral width. The typical number of prepended points
%    found in DMX digitally filtered data is about 60 data points (see below),
%
%    the corresponding 1st order phase correction is ~21,000 degrees.
%    This large linear phase correction can be applied to the transformed data
%    to obtain a normal spectrum. Another, equivalent approach is to time
%    shift the data back to its original position. This results in the need
%    of only a small linear phase shift on the transformed data.
%    There is a question as what to do with the data preceding the actual
%    FID. The prepended data can be simply eliminated with the addition
%    of an equal number of zeros at the end of the FID (left shift). This
%    procedure, however, introduces "frowns" (some have a preference to refer
%    to these as "smiles") at the edge of the spectrum. If the sweep
%    width is fairly wide this does not generally cause a problem. The
%    (proper) alternative is to retain this data by applying a circular left
%    shift of the data, moving the first 60 or so points (see recommendations
%    below) to the end of the FID. This is identical to a Fourier transformation
%    followed by the large linear phase correction mentioned above. The
%    resulting FID is periodic with the last of the data rising to meet the
%    first data point (in the next period). Fourier transform of this
%    data results in an approximately phased spectrum. Further linear
%    phase corrections of up to 180 degrees are necessary. A zero fill applied
%    after a circular shift of the data will cause a discontinuity and thus
%    introduce sinc wiggles on the peaks. The usual correction for DC
%    offset and apodization of the data, if not done correctly, also results
%    in the frowns at the edges of the spectrum.
%
%    In our previous document on Bruker digital filters, we presented deduced
%    rules for calculating the appropriate number of points to be circular left
%    shifted. However, since then, newer versions of hardware (DQD) and software
%    has introduced a new set of values. Depending on the firmware versions
%    (DSPFVS) and the decimation rate (DECIM), the following lookup table will
%    give the circular shift values needed to correct the DMX data. The values
%    of DECIM and DSPFVS can be found in the acqus file in the directory containing
%    the data.
%
%     DECIM           DSPFVS 10       DSPFVS 11      DSPFVS 12
%
%       2              44.7500         46.0000        46.311
%       3              33.5000         36.5000        36.530
%       4              66.6250         48.0000        47.870
%       6              59.0833         50.1667        50.229
%       8              68.5625         53.2500        53.289
%      12              60.3750         69.5000        69.551
%      16              69.5313         72.2500        71.600
%      24              61.0208         70.1667        70.184
%      32              70.0156         72.7500        72.138
%      48              61.3438         70.5000        70.528
%      64              70.2578         73.0000        72.348
%      96              61.5052         70.6667        70.700
%     128              70.3789         72.5000        72.524
%     192              61.5859         71.3333
%     256              70.4395         72.2500
%     384              61.6263         71.6667
%     512              70.4697         72.1250
%     768              61.6465         71.8333
%    1024              70.4849         72.0625
%    1536              61.6566         71.9167
%    2048              70.4924         72.0313
%
%
%    The number of points obtained from the table are usually not integers.  The appropriate procedure is to circular shift (see protocol for details) by the integer obtained from truncation of the obtained value and then the residual 1st order phase shift that needs to be applied can be obtained by multiplying the decimal portion of the calculated number of points by 360.
%
%    For example,
%
%    If DECIM = 32, and DSPFVS = 10,
%    then #points 70.0156
%
%    The circular shift performed on the data should be 70 complex points and the linear
%    phase correction after Fourier transformation is approximately 0.0156*360= 5.62 degrees.
%
%    Protocol:
%
%       1. Circular shift (rotate) the appropriate number of points in the data indicated by
%       the  DECIM parameter. (see above formulae).
%
%       2. After the circular shift, resize the data to the original size minus
%       the number of shifted points. This will leave only the part of the
%       data that looks like an FID. Baseline correction (BC) and/or apodization
%       (EM etc.) should be applied only on this data, otherwise "In come the frowns."
%
%       Since the first part of the data (the points that are shifted) represents
%       negative time, a correct apodization would also multiply the shifted points
%       by a negative time apodization. The data size is now returned to
%       its original size to reincorporate the shifted points. There may
%       still be a discontinuity between the FID portion and the shifted points
%       if thelast point of the FID portion is not at zero. This will cause
%       sinc wiggles in the peaks.
%
%       3. Applying a zero fill to this data will lead to a discontinuity in the data
%       between the rising portion of the shifted points and the zero padding.
%       To circumvent this problem, the shifted points are returned (by circular
%       shift) to the front of the data, the data is zero filled, and then the
%       first points are shifted again to the end of the zero filled data.
%
%       4) The data can now be Fourier transformed and the residual calculated
%       1st order phase correction can be applied.
%
%
%
% For newer versions of the Bruker hardware:
%
%     For firmware versions of 20 <= DSPFVS <= 23 the GRPDLY parameter directly shows the number of
%     points that need to be shifted.
%
%     Thanks for Bruker for supplying this information!
%
%
%
%26-8-'00