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 fidgetpar 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 fidcp $1/$2 $1/$2.bakls $1/$2.bak | cpio -o | cpio -i --swap -uod -An -t dI -v -w8 $1/$2.bak| gawk '{print NR,$1,$2}' > $2.asciigetpars
#!/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 >> $2Octave 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 parametersD=[[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
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
Hiç yorum yok:
Yorum Gönder