Showing posts with label NMR. Show all posts
Showing posts with label NMR. Show all posts

31 January 2014

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.

545. NMR on Debian Wheezy: NMRNotebook by NMRTec

While it's possible -- even easy -- to write your own scripts for processing NMR data (e.g. here and here) there's still a value in having a GUI handy.


Either way, this entry is about NMRNotebook, which ultimately derives from another gratis program called Gifa which I experimented with in the early 00s. Apparently after Gifa v4 (article), NMRTec developed Gifa v5. After that there seems to have been a split in efforts (not a fork) where NMRTec developed NMRNotebook, and other people involved in Gifa developed NPK. NMRNotebook appears to use NPK as the underlying engine, but provides a nice GUI which is written in tk/tcl, java and python.

Note that Gifa 4 and Gifa 5 are still available for download, but require license files to run. I don't know if it's still possible to get licenses for Gifa. Based on my memory of running Gifa (4?) the interface was quite slick, and it had a terminal at the bottom, similar to xwinnmr.

License
Either way, NMRNotebook can be downloaded for free from NMRTec for academic users. You will need to register and then get a license, which can be had by following the instructions on their website and sending off an email with your details.

So I did that.

NOTE that the 1D line fitter is not included in the free version  -- it's priced at 1,000 Euro, which sounds a bit insane to me, in particular if you compare with the price of the software (Euro 100 for a student license, 750 for an industry one). Either way, you can write your own fitter in octave in fifteen minutes..


Installation

There are two files to download: NMRnotebook.sh and NNBMACROS.zip

Run NMRnotebook.sh

sh NMRnotebook.sh
NMRnotebook installer - starting installation... please wait NMRnotebook will be installed in your home directory. Unpacking ... Running post-install script ... done !! To run NMRnotebook, type '~/NMRnotebook/NMRnotebook' to create a shortcut named nnb in your home directory, type 'ln -s ~/NMRnotebook/NMRnotebook ~/nnb' To uninstall simply erase the directory ~/NMRnotebook and ~/.nmrnotebook Some examples can be found in the '~/NMRnotebook/examples' directory Thank you for using NMRnotebook. NMRtec software team
Run nmrnotebook as indicated:
~/NMRnotebook/NMRnotebook









Import the license by going to File/Open and select the nnb file you got via email.

Close the program and start it again.

Unzip the macros by
mkdir ~/work/nnbmacros -p
cp ~/Downloads/NNBMACROS.zip ~/work/nnbmacros
cd ~/work/nnbmacros
unzip NNBMACROS.zip

To run a macro, select it and run it:
To be able to launch nmrnotebook from GNOME create ~/.local/share/applications/nmrnotebook.desktop
[Desktop Entry] Name=NMRnotebook GenericName=nmrnotebook Comment=Software for processing of NMR data Exec=/home/me/NMRnotebook/NMRnotebook Terminal=false Type=Application Categories=Science Version=0.27

Overview
NMRNotebook can open varian, bruker etc. files. Bruker shown here


Window functions, showing EM and interative LB


Several different types of fourier transforms are available

My first spectrum

Zoom works similar to mestrec

Right-clicking brings up a menu that allows you to integrate, label, draw boxes, lines etc.

Normalised integral

NMRnotebook can accept SR values for calibration

Spectrum overview

NMRnotebook only has a few features available, but is sufficient for basic NMR processing. Line-fitting is the most serious omission.
 Overall I find NMRnotebook perfectly adequate for what I would see myself using it for -- simple processing of NMR spectra. Anything more serious and I'd use my own scripts -- but I'd do that anyway in order to be able to trust the data.

23 August 2013

499. Briefly: Drawing NMR sequences using metapost and Mark White's pulse.mp

Since it's Friday afternoon and I'm not likely to get anything useful done in the hour that remains before going home, I might as well put up another post.

Since posting http://verahill.blogspot.com.au/2013/08/498-briefly-drawing-nmr-pulse-sequences.html I've had a look at this: http://www.celos.net/comp/pulses/

And that actually is (almost -- let's not get carried away here) exactly what I have been looking for. The main issue was that it's meant for metapost -- or rather, the main issue was my unfamiliarity with metapost. Anyway, my life is know complete.

So here's how to get started...

First 'install' the metapost script:
sudo apt-get install texlive-metapost
mkdir ~/texmf/metapost
cd ~/texmf/metapost -p
wget http://www.celos.net/comp/pulses/pulses.mp
sudo texhash

Next, time to test-drive it
mkdir ~/tmp/pulse_test -p
cd ~/tmp/pulse_test
vim test.mp
input pulses.mp beginfig(0); initf; startline(rf,"RF"); xline(1); xpulse(0.5,0.5,"90"); ospan(-0.25,-1.0,1.5,"d1"); xline(0.75); xpulse(1.0,0.5,"180"); ospan(-0.5,-1.0,1.5,"d1"); xline(1.0); ospan(0,-1.0,1.5,"vd"); xline(1.5); xacq(2); endfig;
mpost test.mp mptopdf test.0 pdftops -eps test-0.pdf test-0.eps





Now, I somehow suspect you can embed metapost scripts directly in .tex documents, but from my brief testing I haven't quite managed to make it work.

19 December 2012

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

07 August 2011

13. phpSheduleIt - basic installation

phpScheduleIt is a nice web-driven instrument scheduling system, which I first come into contact with as an NMR user at the University of California at Davis, where it is used to manage the scheduling of the instruments in the NMR facility there.

You will need to have a LAMP server set up (see e.g. http://www.howtoforge.com/ubuntu_debian_lamp_server). You will need apache2, php and mysql to play nicely together. Basically follow steps one and two, then sudo /etc/init.d/apache2 restart.

Download phpScheduleIt from http://www.php.brickhost.com/

A basic way to get it up and running goes as follows:

1. create a directory under /var/www. In our example it will be /var/www/nmr
sudo mkdir /var/www/nmr. Unzip the file you downloaded from www.php.brickhost.com and put the files in the /var/www/nmr directory

2. in your terminal, go to /var/www/nmr and chmod o+x all files (including subdirectories)

3. cp /var/www/nmr/config/config.php.new /var/www/nmr/config/config.php

4. Edit it.
change $conf['app']['weburi'] = 'http://localhost/phpScheduleIt';
to

change $conf['app']['weburi'] = 'http://localhost/nmr'; 
edit other relevant items, such as adminEmail, defaultLanguage, timeFormat, emailType, defaultPassword etc,

The admin email is important, since this will become the administrator account.

5. go to http://localhost/nmr/install
If nothing happens, or you get an error message, it's time to start trouble-shooting. Going to the terminal and
cd /var/www/nmr/install
php index.php
might give some useful information. Most likely you have either a) made an error in config/config.php, or b) you haven't chmod:d the subdirectories.

6.  If all goes well, you are asked for your root user name and root mysql password (which you set during the installation). Thus, login as root with the correct password. Once you're in, hit the button saying 'create tables'. If all went well you can now delete the /var/www/nmr/install directory.

7. Go to http://localhost/nmr
You should be greeted with a log in window. Important: select create new user and enter the email address you gave as the administrator email in the config.php file. This will then automatically become an administrator account. It has nothing to do with the order you create accounts in - it's all about what email address you put in the config.php file. 

8.  Once you've set up the account you are able to create instruments, user groups, booking restrictions etc. 
Anyone who wants to use the facilities can create their own new registered users, BUT only the administrator(s) can assign resources to them. 

9. Be careful with your GMT setting...

This page has some information as well: http://www.m-osaka.com/jp/plaza/readme.html

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