20 July 2011

7. Processing 1D Bruker nmr data

Bruker 1D binary NMR files can be processed using a combination of cat, grep, sed, gawk and od, together with python and octave (w/ octave-optim) for some fancy line-fitting.

 brukdig2asc:
 #!/bin/bash
#usage: brukdig2asc
SW=`cat acqus | grep 'SW_h' | sed 's/\=/\t/g' | gawk '{print $2}'| tr -d '\n'`
TD=`cat acqus | grep 'TD=' | sed 's/\=/\t/g' | gawk '{print $2}'| tr -d '\n'`
O=`cat acqus | grep '$O1=' | sed 's/\=/\t/g' | gawk '{print $2}'`
SFO=`cat acqus | grep 'SFO1=' | sed 's/\=/\t/g' | gawk '{print $2}'`
#TIME=16384
#SWEEP=23809.5238095238
#AQ=`echo "1/(23809.5238095238/(16384/2))" | bc -lq`
cp fid fid.bin
ls fid.bin | cpio -o | cpio -i --swap -u
od -An -t dI -v -w8 fid.bin| gawk '{print NR,$1,$2}'| sed '1,64d' >fid.asc1
pynmr $SW $TD $O $SFO
makespec

pynmr:
#!/usr/bin/python2.6
import sys
#print str(sys.argv)
sweepwidth=float(sys.argv[1])
nopts=int(sys.argv[2])
centrefreq=float(sys.argv[3])
basefreq=float(sys.argv[4])

aq=1/(sweepwidth/(nopts/2))
#print str(sweepwidth),str(nopts)
f=open('fid.asc1','r')
g=open('fid.asc','w')
for line in f:
    line=line.rstrip('\n')
    line=line.split(' ')
#    print line
    freq=float(line[0])/(nopts/2)*sweepwidth+(centrefreq-sweepwidth/2)
    line[0]=(float(line[0])/(nopts/2))*aq
    g.write(str(line[0])+'\t'+str(line[1])+'\t'+str(line[2])+'\t'+str(freq)+'\n')
f.close
g.close

makespec:

#!/bin/bash
octave --silent --eval "fid=load('fid.asc');
#make xaxis
[nopts b]=size(fid);
aq=max(fid(:,1));
sw=nopts/aq;
freqx=linspace(0,sw,nopts)';

#apodizing
lb=5/10000;
fid(:,2)=fid(:,2).*exp(-lb.*freqx);
fid(:,3)=fid(:,3).*exp(-lb.*freqx);

#phasing
spec=[fid(:,1) real(fftshift(fft(fid(:,2)+i*fid(:,3)))) imag(fftshift(fft(fid(:,2)+i*fid(:,3))))];
[a b]=size(spec); spec(a/2,2:3)=[0 0];
phc=linspace(0,2*pi,180);
maxsig=0;k=1;
for n=1:180;
        localmax=max( real( (spec(:,2)+i*spec(:,3)).*exp(i*phc(n)) ));
        if (localmax>maxsig)
                maxsig=localmax;
                k=n;
        endif
endfor;
#simple baseline
absd=inline('m+t*0','t','m');
guess=0;
[f m kvg iter corp covp covr stdresid z r2]=leasqr(fid(:,4),real((spec(:,2)+i*spec(:,3)).*exp(i*phc(k))),guess,absd);
#disp(m)
#disp(sqrt(diag(covp)))

#make spectrum
spectrum=[fid(:,4) real((spec(:,2)+i*spec(:,3)).*exp(i*phc(k)))-m imag((spec(:,2)+i*spec(:,3)).*exp(i*(phc(k)+pi/2)))-m];

#fitting
pkg load optim
[a b]=max(spectrum(:,2));
centre=fid(b,4);
guess=[10 max(spec(:,2))]; #centre width height
#disp(guess)
lorentzian=inline('p(2)*(1/pi)*(p(1)/2)./((t-centre).^2+(0.5*p(1))^2)','t','p');
[f p r2]=leasqr(fid((b-150):(b+150),4),spectrum((b-150):(b+150),3),guess,lorentzian);

#filter out artefacts from fitting set
filtered=[0 0];
res=floor((max(fid(:,4))-min(fid(:,4)))/nopts);
for l=(b-ceil(5*p(1)/res)):(b+5*ceil(p(1))/res)
delta=lorentzian(fid(l,4),p)-spectrum(l,2);

if (delta>(lorentzian(fid(l,4),p))/1.2)
# do nothing
else
filtered=[filtered; fid(l,4) spectrum(l,2)];
endif
endfor

filtered=[ filtered(2:size(filtered(:,2)),1) filtered(2:(size(filtered(:,2))),2)  ];
[f p r2]=leasqr(filtered(:,1),filtered(:,2),p,lorentzian);

#disp(p')
#disp(r2)
params=[centre centre/67.8 max(lorentzian(fid(:,4),p)) p(1) 1.000 p(2)];
disp(params)
#save
spex=[fid(:,4) real((spec(:,2)+i*spec(:,3)).*exp(i*phc(k)))-m imag((spec(:,2)+i*spec(:,3)).*exp(i*(phc(k)+pi/2)))-m lorentzian(fid(:,4),p)];
save spectrum.dat spex;"

9 comments:

  1. Hi. I apologize for digging up an old post, but I was wondering whether you could clarify the part where you split the data into real and imaginary components. I recall reading that the Bruker ADC alternates sampling real and imaginary points so the two don't actually line up in the time domain. Do you know if this is correct? If so, would this have any impact on the calculation? Thanks.

    ReplyDelete
    Replies
    1. No worries.

      The separation of the imaginary and real components is done here:
      od -An -t dI -v -w8 fid.bin| gawk '{print NR,$1,$2}'| sed '1,64d' >fid.asc1

      The data is definitely /stored/ using the same time coordinate, so there's no hint that Im and Re are misaligned. Having said that, I actually don't know.

      As for effect: try this octave script:
      ########################################
      nopts=1*1024; acqu=10^-0;
      t=linspace(0,acqu,nopts);
      nu0=0;
      freqaxis=linspace(nu0-1/(2*(acqu/nopts)),nu0+1/(2*(acqu/nopts)),nopts);
      fid=zeros(nopts,3);
      fid(:,1)=t';

      offset=10^(-4)

      nu=0; T1=0.06; T2=0.06;
      ft= @(t) [-sin(2*pi*(nu0-nu).*t).*exp(-(t./T2)) cos(2*pi*(nu0-nu).*(t+offset)).*exp(-((t+offset)./T2)) (1-exp(-t./T1))];

      sim=ft(t');

      fid(:,2)=sim(:,1); fid(:,3)=sim(:,2); fid(:,4)=sim(:,3);

      spectrum=zeros(nopts,1);
      spectrum(:,1)=fftshift((fft( fid(:,2)+i*fid(:,3) )));
      plot(freqaxis,imag(spectrum(:,1)))

      ########################
      You can change the offset which will shift the misalignment between the imaginary and real fids.

      The effect is a separation of magnetisation -- you get an dispersive signal at -freq in addition to an attenuated signal at freq. Try changing the offset from 10^(-6) to 10^(-3).

      Delete
    2. I'll be writing a new set of scripts for processing (I've learned a bit more in the past year), so check back in a week or so and I'll hopefully have a set of minimal scripts that mimimic the xwin-nmr workflow (i.e. one script for binary fid to ascii fid, one for fourier transform (ft), one for phase correction (pk) etc.)

      Anyway, what I had forgotten is that there's obviously no time information in the fid -- just data points. Regardless, using the octave script above you can get a functional understanding of what an offset would do.

      tail fid.ascii
      16375 792 522
      16376 8718 9299
      16377 10011 1589
      16378 12421 -5299
      16379 12131 16347
      16380 4781 3898
      16381 9357 1508
      16382 -5350 4289
      16383 2749 5987
      16384 5414 10767

      Delete
    3. Hey there lindqvist, did you write a new script to convert the fid file to a plotable spectre?
      I am trying to create a shell script to process it for 3k files, but I am having difficulties on generate the spectre.

      Delete
    4. Roberto, I have a more recent post at http://verahill.blogspot.com.au/2012/12/bruker-1d-processing-using-octavematlab.html

      You obviously don't have to use all the different scripts such as zero filling, phase correction etc.

      What part are you having issues with? What do you get? What works so far?

      Delete
    5. First of all thank you very much by answering so fast!! :)

      So, basically what I did was run the command "od -An -i -v -w4 ./$file " where $file is the variable which I am iterating on a find logic.

      I am saving the result of it on a new file.txt.

      After that I load the txt file at octave and run the command "abs(fft(FILE))" in order to have an array with the spectre data.

      The bad new is that the plotted result doesn't seems nothing like a real nmr spectre... :(

      Delete
    6. Two thoughts:
      * if you aren't already, make sure that you are reading the data as being complex (i.e. real + imag part). Otherwise you'll get two peaks on either side of the centre of the spectrum.
      * try using fftshift e.g. abs(fftshift(fft(fid))) ) -- this will work if you find that the edges of the spectum contain the signal data (spectrum looks like a smiley face), and will shift the data to the centre.

      Delete
    7. Hey there, here am I again hehe...

      Following your tips on my own script I am getting a mirrored graphic which doesn't seems a spectrum yet.

      Is there any way that I can send you the picture of the graphic that I am getting?

      Delete
  2. This comment has been removed by a blog administrator.

    ReplyDelete