546. NMR on Debian Wheezy: ACD/LABS NMR under Wine

I'm currently exploring options for processing NMR data on linux (see e.g. here for nmrnotebook). There are four categories of programs:
* proprietary/for-cost running natively on linux

* proprietary/non-open source but gratis running natively on linux
* gratis and open source running natively on linux

* windows programs that happen to run ok in wine

ACD/LABS belongs to the last category, and as such is not an ideal solution. However, at this point I'm willing to try just about anything.

So yes, this relies on you having wine installed, but sometimes we need to be pragmatic like Torvalds rather than dogmatic like Stallman (the world needs both of them).

---

Either way:

First register at http://www.acdlabs.com/account/register.php?redirect=/account/logout.php

Then download

 

---

Installation:

Install by running

wine ~/Downloads/nmrproc_academia12.exe

Go through the installation steps:

 

---

Running the software

For most of the steps below there's little need to provide a commentary:

I could only get it to fit by hitting 'Auto'

However, 'Auto' did respect my choice of function e.g. lorentzian vs mixed.

You can export the data as simple ascii files

 

In conclusion I must say that I'm fairly happy with this piece of software. I haven't used it for production purposes, so I can't vouch for it spitting out reliable data though. It's definitely worth exploring though.

294. Bruker 1D processing using octave/matlab

I wanted a set of scripts that behaved a little bit like the commands in bruker xwin-nmr/topspin, so that I could do some quick processing for visual inspection without having to do too much coding.

I also wanted some simple modules that I can plug into automated processing routines for large numbers of spectra (e.g. when doing kinetics).

So here are a few simple octave routines which should work in matlab as well. They won't change the world, but should be good enough for some basic 1D processing.

Because of the groupdelay (GRPDLY) used in by Bruker (see e.g. here (16th of June post) and here), you need to use the bruk2ana converter. There's little science behind the values which are applied since they are hardware specific.

Example workflow:

./bin2ascii experiment_1/1 fid
getpar experiment_1/1 my.par

octave:1> [fid,pars]=loadfid('fid.ascii','my.par');
octave:2> [zfid,pars]=zf(fid,pars);
octave:3> [test,phc1]=bruk2ana(zfid,pars);
octave:4> test=em(test,0.5);
octave:5> plot(test(:,1),test(:,3));
octave:6> spectrum=ft(test,pars);
octave:7> phased=apk(spectrum,phc1);
octave:8> final=absd(spectrum);
octave:9> pltspec(final)

---

Linux shell-scripts

---

bin2ascii

#!/bin/bash
#bin2ascii dir fid
cp $1/$2 $1/$2.bak
ls $1/$2.bak | cpio -o | cpio -i --swap -u
od -An -t dI -v -w8 $1/$2.bak| gawk '{print NR,$1,$2}' > $2.ascii

getpars

 #!/bin/bash
 #getpars $1 $2
 # $1 is the location (directory or .) and $2 is the root of the output file name
 SW=`cat $1/acqus | grep 'SW_h' | sed 's/\=/\t/g' | gawk '{print $2}'| tr -d '\n'`
 TD=`cat $1/acqus | grep 'TD=' | sed 's/\=/\t/g' | gawk '{print $2}'| tr -d '\n'`
 O=`cat $1/acqus | grep '$O1=' | sed 's/\=/\t/g' | gawk '{print $2}'`
 SFO=`cat $1/acqus | grep 'SFO1=' | sed 's/\=/\t/g' | gawk '{print $2}'`
 DECIM=`cat $1/acqus | grep 'DECIM=' | sed 's/\=/\t/g' | gawk '{print $2}'`
 DSPFVS=`cat $1/acqus | grep 'DSPFVS=' | sed 's/\=/\t/g' | gawk '{print $2}'`
 echo $SW > $2
 echo $TD >> $2
 echo $O >> $2
 echo $SFO >> $2
 echo $DECIM >> $2
 echo $DSPFVS >> $2

---

Octave scripts

---

loadfid.m

function [fid,pars]=loadfid(infile,parfile)
%%Usage: [fid,pars]=loadfile(infile,parfile)
%%reads a pts, re, im ascii array 
%%generated by bin2ascii
%%and a parfile generated with genpar
 fid=load(infile);
 pars=load(parfile);
 t=linspace(0,(1/(pars(1)/(pars(2)/2))),pars(2)/2);
 fid=[t' fid(:,2) fid(:,3)];
end

zf.m

function [newfid,pars]=zf(fid,pars)
%% Usage:[newfid,updatedpars]=zf(fid,pars)
%% Doubles the number of points by zerofilling
 dims=size(fid);
 newfid=[fid' zeros(3,dims(1))]';
 pars(2)=pars(2)*2;
end

bruk2ana

function [fid,phc1]=bruk2ana(fid,pars)
%% Usage: [fid,phc1]=bruk2ana(fid,pars)
%% where phc1 is the first-order phase correction
%% Using https,//nmrglue.googlecode.com/svn-history/r44/trunk/nmrglue/fileio/bruker.py
%% and https,//ucdb.googlecode.com/hg/application/ProSpectND/html/dmx_digital_filters.html
%% The short version is: bruker fid data needs pre-processing and it's hardware dependent

%%D contains the digital filter parameters
D=[[10,2, 44.75];
[10,3, 33.5];
[10,4, 66.625];
[10,6, 59.083333333333333];
[10,8, 68.5625];
[10,12, 60.375];
[10,16, 69.53125];
[10,24, 61.020833333333333];
[10,32, 70.015625];
[10,48, 61.34375];
[10,64, 70.2578125];
[10,96, 61.505208333333333];
[10,128, 70.37890625];
[10,192, 61.5859375];
[10,256, 70.439453125];
[10,384, 61.626302083333333];
[10,512, 70.4697265625];
[10,768, 61.646484375];
[10,1024, 70.48486328125];
[10,1536, 61.656575520833333];
[10,2048,70.492431640625];
[11,2, 46.];
[11,3, 36.5];
[11,4, 48.];
[11,6, 50.166666666666667];
[11,8, 53.25];
[11,12, 69.5];
[11,16, 72.25];
[11,24, 70.166666666666667];
[11,32, 72.75];
[11,48, 70.5];
[11,64, 73.];
[11,96, 70.666666666666667];
[11,128, 72.5];
[11,192, 71.333333333333333];
[11,256, 72.25];
[11,384, 71.666666666666667];
[11,512, 72.125];
[11,768, 71.833333333333333];
[11,1024, 72.0625];
[11,1536, 71.916666666666667];
[11,2048, 72.03125];
[12,2, 46. ];
[12,3, 36.5];
[12,4, 48.];
[12,6, 50.166666666666667];
[12,8, 53.25];
[12,12, 69.5];
[12,16, 71.625];
[12,24, 70.166666666666667];
[12,32, 72.125];
[12,48, 70.5];
[12,64, 72.375];
[12,96, 70.666666666666667];
[12,128, 72.5];
[12,192, 71.333333333333333];
[12,256, 72.25];
[12,384, 71.666666666666667];
[12,512, 72.125];
[12,768, 71.833333333333333];
[12,1024, 72.0625];
[12,1536, 71.916666666666667];
[12,2048, 72.03125];
[13,2, 2.75]; 
[13,3, 2.8333333333333333];
[13,4, 2.875];
[13,6, 2.9166666666666667];
[13,8, 2.9375];
[13,12, 2.9583333333333333];
[13,16, 2.96875];
[13,24, 2.9791666666666667];
[13,32, 2.984375];
[13,48, 2.9895833333333333];
[13,64, 2.9921875];
[13,96, 2.9947916666666667];];

 h=find(D(:,2)==pars(5));
 j=find(D(h,1)==pars(6));
 magickey=D(h(j),3);
 chop=floor(magickey);

 phc1=(magickey-chop)*2*pi; %the first-order phase correction gets mangled by bruker

 tmp=size(fid); %matlab workaround. rows/columns would be more elegant
 
 newfid=[fid(chop:tmp(1),2:3)' fid(1:chop-1,2:3)']';
 fid=[fid(:,1) newfid(:,1) newfid(:,2)];
 
end

em.m

function fid=em(fid,lb)
%%Usage: fid=em(fid,lb)
%%Exponential multiplication window function
%%Increases Signal-to-noise at the expense
%%of resolution
 fid(:,2)=fid(:,2).*exp(-lb.*fid(:,1));
 fid(:,3)=fid(:,3).*exp(-lb.*fid(:,1));
end

gm.m

function fid=gm(fid,lb)
%%Usage: fid=gm(fid,lb)
%%Gaussian multiplication window function
%%Increases Signal-to-noise at the expense
%%of resolution
 fid(:,2)=fid(:,2).*exp(-(lb.*fid(:,1)).^2);
 fid(:,3)=fid(:,3).*exp(-(lb.*fid(:,1)).^2);
end

de.m

function fid=de(fid,lb,gm)
%%Usage fid=de(fid,lb,gm)
%%Double-exponential window function
%%Increases resolution at the expense of
%%signal-to-noise
 at=max(fid(:,1));
 defun= @(lb,gm,t) (exp(-(t.*lb-gm*at))).^2;
 fid(:,2)=fid(:,2).*defun(lb,gm,fid(:,1));
 fid(:,3)=fid(:,3).*defun(lb,gm,fid(:,1));
end

traf.m

function fid=traf(fid,lb)
%%Usage: fid=traf(fid,lb)
%%TRAF window function
%%Increases resolution at the expense of
%%the signal-to-noise
 at=max(fid(:,1));
 traffun= @(lb,t) (exp(-t.*lb)).^2./((exp(-t.*lb)).^3+(exp(-at*lb)).^3);
 fid(:,2)=fid(:,2).*traffun(lb,fid(:,1));
 fid(:,3)=fid(:,3).*traffun(lb,fid(:,1));
end

ft.m

function spectrum=ft(fid,pars)
%%Usage: spectrum=ft(fid,pars)
%% Spectrum is a complex array with the frequency in
%%the first column and the real and imaginary parts
%%in the second column
%%pars(3)=centrefreq, pars(1)=SW
 spectrum=fftshift(fft(fid(:,3)+i*fid(:,2)));
 tmp=size(spectrum);%matlab workaround
 freq=linspace(pars(3)+pars(1)/2,pars(3)-pars(1)/2,tmp(1));
 spectrum=[freq' spectrum];
endfunction

apk.m

function spectrum=apk(spectrum,phc1)
%%Usage spectrum=apk(spectrum,phc1)
%%Spectrum is a complex matrix with
%%the frequency in the first column
%%and the complex spectrum in the 
%%second column. phc1 is the first order
%%phase correection

 tmp=size(spectrum);
 m=720;
 ph=linspace(-2*pi,2*pi,m);
 maxsig=0;k=1;
 minsig=-inf;
 for n=1:m;
  spex=real( (spectrum(:,2)).*exp(i*(ph(n)+phc1*i/tmp(1))) );
  localmin=min(spex);
        localmax=max(spex);
        if (localmin>minsig) 
                minsig=localmin;
                k=n;
        end
 end
 ph0=ph(k);
 spectrum(:,2)=spectrum(:,2).*exp(i*(ph0+phc1*i/tmp(1)));
end

altapk.m

function [spectrum,ph]=altapk(spectrum,phc0,phc1)
%%Usage -spectrum,ph]=altapk(spectrum,phc0,phc1)
%%Spectrum is a complex matrix with the frequency in the first column
%%and the complex spectrum in the second column. phc0 and phc1 are the first order
%%phase correction parameters, respectively, and are used as initial guesses.
%%This is an implementation of Chen, Weng, Goh and Garland, J. Mag. Res., 2002, 158, 164-168 and depends on entropy.m.

    ph=[phc0;phc1];
        ph=minimize("entropy",{ph,spectrum});

        %compute spectrum with optimal phase params
        pts=linspace(1,size(spectrum(:,2),1),size(spectrum(:,2),1));
        phi=(ph(1)+ph(2).*pts./max(pts))';
        spectrum(:,2)=spectrum(:,2).*exp(i*phi);
end

entropy.m

function E=entropy(ph,spectrum)
%%Used by altapk.m
    pts=linspace(1,size(spectrum(:,2),1),size(spectrum(:,2),1));
        penalty=5.53;

    phi=(ph(1)+ph(2).*pts./max(pts))';
        size(phi);
        spectrum(:,2)=spectrum(:,2).*exp(i*phi);
        R=real(spectrum(:,2));
        size(R);
    Rm=firstderiv(R);
        size(Rm);
    h=abs(Rm)/sum(abs(Rm));
        size(h);

        negs= imag((R).^(1/2));
        negs(find(negs>1))=1;
    P= @(R) penalty.*sum((negs).*R.^2);

    E=-sum(h.*log(h))+P(R);
end

absd.m

function spectrum=absd(spectrum)
%%Usage spectrum=absd(spectrum)
%%Simple (linear) baseline correction
        bsline=@(m) sum(abs(real(spectrum(:,2))-m));
        guess=0;
        p=0;
        newm=minimize(bsline,p);
        spectrum(:,2)=spectrum(:,2)-newm;
end

pltspec.m

function pltspec(spectrum)
%%Usage: pltspec(spectrum)
%%Where spectrum is a complex matrix
%%with the frequency in the first column
%%and the complex spectrum (a+i*b) in the
%%second column
 plot(spectrum(:,1),real(spectrum(:,2)))
end

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;"