metab_synth.m

The goal of SPID is to provide the user with tools capable to simulate, preprocess, process and clas

function [metabolite] = metab_synth(strFilename, fs, transmitterfreq, ref, tau, deltaf)

%function [metabolite] = metab_synth(strFilename, fs, transmitterfreq, ref, tau, deltaf,assign_damping) 
% strFilename='metabolite_specs_6.7/choline'; fs=7200; transmitterfreq=6E8; ref=4.7; tau=0; deltaf=[0]



%sparcify matrices!

fprintf('\nMaking %s...', strFilename);
time1=cputime;

warning off all
metab=xlsread(strFilename);
warning on all
%SET DAMPING CONSTANT FOR METABOLITE
%t2=0.02;
%t2=0.05;
%t2 = 0;

%SET SPIN ECHO DELAY
%tau=0.000;
%tau=0.030;
%tau=0.136;

%************************************************************************
%transmitterfreq=field_strength*10^7*42.58;
%transmitterfreq=5.998E9;
%transmitterfreq=5E9;

%t2=t2*2*pi;
tau=tau*2*pi;
%ref=4.81;
%ref=4.585;
bandwidth=fs;
sampint=1/bandwidth;
%zfill=8;

spin_num=length(metab)-2;

%construct chem shift vector PPM>Freq
for n=1:spin_num
    chem_shift(n)=(-metab(n,2)+ref)*(transmitterfreq/(10^6));
end
chem_shift=chem_shift.';




% %construct spin code
% for n=1:spin_num
%     spin_code(n)=metab(n,1);
% end
% spin_code=spin_code';

%construct equiv array
 for n=1:spin_num
     equiv_array(n)=metab(n,1);
 end
 equiv_array=equiv_array';

 % If the function has been passed an offset to add to each group...
 if nargin == 6      
     if length(deltaf)~=max(equiv_array)
%         fprintf('ERROR Incorrect number of deltafs');
         pause;
     end

     %update any chemical shift optimisation
     for r=1:max(equiv_array)
         for n=1:spin_num
             if equiv_array(n)==r
                 chem_shift(n)=chem_shift(n)+deltaf(r);
             end
         end
     end
 end
 
 %assignin('base', 'equiv_array', equiv_array);

 %construct half J coup matrix
for n=1:spin_num-1
    for m=n+1:spin_num
        j_coup(m,n)=metab(m,n+2)/(2*pi);
    end
end

%construct full J coup matrix
j_coup_full=zeros(spin_num+1);
for n=1:spin_num-1
    for m=n+1:spin_num
        j_coup_full(m+1,n+1)=metab(m,n+2)/(2*pi);
    end
end
j_coup_full=j_coup_full+j_coup_full.';
j_coup_full(1,2:spin_num+1)=chem_shift.';
j_coup_full(2:spin_num+1,1)=chem_shift.';

% %change to look for equivalences
% for n=2:spin_num+1
%     j_coup_cut1(n-1,1:spin_num+1)=j_coup_full(n,1:spin_num+1);
% end
% for n=1:spin_num
%     q=0;
%     for m=1:spin_num+1
%         %if j_coup_cut1(n,m)~=0;
%         if m~=n+1
%             q=q+1;
%             j_coup_cut2(n,q)=j_coup_cut1(n,m);
%         end   
%     end
% end
% j_coup_cut2=j_coup_cut2';
% %find equivalent spins
% p=1;
% equiv_array(1)=p;
% for m=1:spin_num-1
%     if j_coup_cut2(1:spin_num,m)==j_coup_cut2(1:spin_num,m+1)  
%         equiv_array(m+1)=p;   
%     else
%         p=p+1;
%         equiv_array(m+1)=p;
%     end
% end

%break
%construct chemical shift part of hamiltonian
chem_shift_ham=0;
for n=1:spin_num
    chem_shift_ham=chem_shift_ham+chem_shift(n)*prodop(spin_num,'z',n);
end

%construct jcoupling part of hamiltonian
j_coup_ham=0;
for n=1:spin_num-1
    for m=n+1:spin_num
        j_coup_ham=j_coup_ham+2*pi*j_coup(m,n)*prodop(spin_num,'x',n)*prodop(spin_num,'x',m);
        j_coup_ham=j_coup_ham+2*pi*j_coup(m,n)*prodop(spin_num,'y',n)*prodop(spin_num,'y',m);
        j_coup_ham=j_coup_ham+2*pi*j_coup(m,n)*prodop(spin_num,'z',n)*prodop(spin_num,'z',m);
    end
end

total_ham=(chem_shift_ham+j_coup_ham);
%total_ham = diag(diag(total_ham));
%find eigenvalues and vectors
[eig_vec,eig_val]=eig(total_ham);

%fprintf('calculated hamiltonian\n');

for n=1:2^spin_num
    eig_val_array(n)=eig_val(n,n);
end

%fprintf('calculating amplitudes\n');

%does this twice for choline???
for r=1:max(equiv_array)
    q=0;
    %pulse seq
    pulse_seq=0;
    for n=1:spin_num
        if equiv_array(n)==r
            pulse_seq=pulse_seq+(prodop(spin_num,'z',n));
        end
    end
    %first90x
    pulse_seq=rotop(spin_num,'x',90)*pulse_seq*rotop(spin_num,'x',-90);
    %evolve
    pulse_seq=expm(-i*total_ham*tau/2)*pulse_seq*(inv(expm(-i*total_ham*tau/2)));
    %180x
    pulse_seq=rotop(spin_num,'x',180)*pulse_seq*rotop(spin_num,'x',-180);
    %evolve
    pulse_seq=expm(-i*total_ham*tau/2)*pulse_seq*(inv(expm(-i*total_ham*tau/2)));
    %Ix operator sum
    Ixsum=0;
    for p=1:spin_num
        Ixsum=Ixsum+prodop(spin_num,'x',p);
    end
    Ixsum=sparse(Ixsum);
    coherence=-eig_vec'*pulse_seq*eig_vec;
    coupled_coherence=eig_vec'*Ixsum*eig_vec;
    a=2*i*exp(-i*(pi))*coherence.*coupled_coherence;

    for n=1:2^spin_num
        for m=1:2^spin_num
            %ignore smallest peaks
            if abs(a(n,m))>0.01
                q=q+1;
                a_cut(q,r)=a(n,m);
                peak_freq(q,r)=(-eig_val_array(m)+eig_val_array(n))*(2*pi);
                %remove +1 peaks
                if peak_freq(q,r)*chem_shift(find(equiv_array == r, 1))<0
                    a_cut(q,r)=0;
                end
            end
        end
    end
    arraylength(r)=q;
end
    
%fprintf('calculated frequencies\n');
cpt_time=cputime-time1;
%fprintf('calculation took %6.2f seconds\n',cpt_time);

a_cut=a_cut/(2^spin_num);
%%%%%%% EDITTED

% biggest=zeros(size(a_cut,2),1);
% 
% for g = 1:1:size(a_cut,2)
%     group = a_cut(:,g);
%     
%     for n=1:length(group)
%         if (angle(group(n))>-3.14) || (angle(group(n))<-3.15)
%             if abs(group(n))>biggest(g)
%                biggest(g)=abs(group(n));
%             end
%         end
%     end
%     
%     for n=1:length(group)
%         if abs(group(n))<biggest(g)*0.3;%1.0001
%             a_cut(n,g)=0;
%         end
%     end      
% end
%%%%%%%%%%%%%%%%%%%%%%%%


% make the output matrix of the parameters
metabolite = [];
for r=1:max(equiv_array)
    for n=1:arraylength(r)
      if abs(a_cut(n,r)) > eps
	metabolite = [ metabolite; ... 
		       r a_cut(n,r) (peak_freq(n,r)/(2*pi)).' 0 ...
		     ];
      end
    end
end


% if nargin > 6
%    % metabolite(:,4) = 3 + 7*rand(1,1)*ones(size(metabolite(:,4)));
%     metabolite(:,4) = 6.5*ones(size(metabolite(:,4)));
% end
fprintf('done.');