

%import 1D varian spectroscopy data file and execute FFT

fid = fopen('fid','r','ieee-be');

np=524288; %total # points
sfrq=138.2708416; %actual center freq in MHz
sw=500000  ; %spectral width in Hz
lb=0; %line broadening in Hz
ph0 = 15*(pi/180); %0th order phase correction in radians
ph1 = 5*(pi/180); %1st order phase correction in radians


% Read datafileheader
nblocks   = fread(fid,1,'int32');
ntraces   = fread(fid,1,'int32');
np        = fread(fid,1,'int32');
ebytes    = fread(fid,1,'int32');
tbytes    = fread(fid,1,'int32');
bbytes    = fread(fid,1,'int32');
vers_id   = fread(fid,1,'int16');
status    = fread(fid,1,'int16');
nbheaders = fread(fid,1,'int32');


% s_data    = bitget(status,1);
% s_spec    = bitget(status,2);
s_32      = bitget(status,3);
s_float   = bitget(status,4);
% s_complex = bitget(status,5);
% s_hyper   = bitget(status,6);

% reset output structures

RE = zeros(np/2,ntraces*nblocks);
IM = zeros(np/2,ntraces*nblocks);

inx = 1;
for b = 1:nblocks
    sprintf('read block %d\n',b);
    % Read a block header
    scale     = fread(fid,1,'int16');
    bstatus   = fread(fid,1,'int16');
    index     = fread(fid,1,'int16');
    mode      = fread(fid,1,'int16');
    ctcount   = fread(fid,1,'int32');
    lpval     = fread(fid,1,'float32');
    rpval     = fread(fid,1,'float32');
    lvl       = fread(fid,1,'float32');
    tlt       = fread(fid,1,'float32');
    for t = 1:ntraces
        %We have to read data every time in order to increment file pointer
        if s_float == 1
            data = fread(fid,np,'float32');
        elseif s_32 == 1
            data = fread(fid,np,'int32');
        else
            data = fread(fid,np,'int16');
        end
        RE(:,inx) = data(1:2:np);
        IM(:,inx) = data(2:2:np);
        inx = inx + 1;     
    end %trace loop
end  % done reading one block


fidr=RE;%real part of FID
fidi=IM;%imaginary part of FID

%parameters from the 'procpar' file:

np=524288; %total # points

for i=1:np/2;    %perform line broadening pt by pt 
    fidr(i)=fidr(i)*exp(-t/5000);
    fidi(i)=fidi(i)*exp(-t/5000);
end

np=524288; %total # points

fidcom=fidr + 1i*fidi;       %recombine Re and Im parts    

for n=1:(np/2)
    fidcom(n)=fidcom(n)*exp(1i*(ph0+n*ph1));
end

spectrum=abs(fftshift(fft2(fidcom)));
phasedSpectrum = fftshift(fft2(fidcom));



PPMpp=((sw/sfrq))/(np/2);
axis=PPMpp:PPMpp:PPMpp*np/2;
axis=axis-PPMpp*np/4;


figure; 
plot(axis,flip(abs(phasedSpectrum)),'LineWidth',1.0);  %plot the chemical shift spectrum
title('^129Xe spectrum','FontSize',16,'FontName','Arial'); 
bx=gca;
bx.XAxis.FontName='Arial';
bx.XAxis.FontSize=16;
bx.YAxis.FontName='Arial';
bx.YAxis.FontSize=16;
xlabel('Chemical Shift (ppm)','FontSize',16,'FontName','Arial'); 
ylabel('a.u.','FontSize',16,'FontName','Arial');
set(bx, 'xdir', 'reverse');



