MSBC多光谱校正算法处理拉曼

% % the simulation in the manuscript clc; clear; close all; %% ============== generate data =============================== N = 3351; t = (1: N)'; x = 4000 : -1 : 650; % wavenumber % eight Gaussian peaks with different widths sig1 = 50; sig2 = 10; sig3 = 25; d = 0.3* exp( -(t-600 ).^2/(sig1^2*2)) + 0.1* exp( -(t-1500).^2/(sig2^2*2))+... 0.2* exp( -(t-1540).^2/(sig2^2*2)) + 0.15*exp( -(t-1570).^2/(sig2^2*2))+... 0.15*exp( -(t-2100).^2/(sig3^2*2)) + 0.2* exp( -(t-2450).^2/(sig2^2*2))+... 0.3* exp( -(t-2500).^2/(sig2^2*2)) + 0.2* exp( -(t-2800).^2/(sig2^2*2)); % piecewise functions u = t / N; t1 = u(1:1000); t2 = u(1001:2000); t3 = u(2001:3000); t4 = u(3001:3351); y1 = 0.1*t1; % linear function y2 = -0.5*(t2-1500/N).^2; y2 = y1(end)-y2(1)+y2; % quadratic function y3 = 1*(t3-2201/N).^2 ; y3 = y2(end)-y3(1)+y3; % quadratic function y4 = 0.2*exp(t4); y4 = y3(end)-y4(1)+y4; % exponential function y12 = 0.13*t1; y22 = -0.3*(t2-1500/N).^2; y22 = y12(end)-y22(1)+y22; y32 = 1*(t3-2201/N).^2 ; y32 = y22(end)-y32(1)+y32; y42 = 0.25*exp(t4); y42 = y32(end)-y42(1)+y42; % construct two synthetic baselines base1 = [y1; y2; y3; y4] ; base2 = [y12; y22; y32; y42] ; % two baseline drift spectra data1 = d + base1; data2 = d + base2; data = [data1 data2]'; base = [base1 base2]'; clear y*; clear s*; clear t*; clear m; % figure show original spectra and drifted spectra % % Figure 1 figure('Position',[300,150,650,300]); subplot(121); plot(x,d,'k'); text(4700, 0.34, '(a)'); set(gca,'XDir','reverse'), set(gca,'XLim',[650 4000],'Ylim',[0,0.36]); xlabel('Wave number/cm^-^1'); ylabel('Absorbance'); subplot(122); plot(x,data1,'k--'); hold on; plot(x,data2,'k'); hold off; text(4700, 0.34, '(b)') set(gca,'XDir','reverse'),set(gca,'XLim',[650 4000],'Ylim',[0,0.36]) xlabel('Wave number/cm^-^1'); ylabel('Absorbance') %% =============================================================== %% ======= Multiple spectra baseline correction ================= lambda = 150; mu1 = 0.8*1e8; mu2 = 1e8; mu0 = [mu1 mu2]; p = 0; [z a w] = MSBC(data, lambda, mu0, p); RMSE1 = norm(base1-z(1,:)') RMSE2 = norm(base2-z(2,:)') RMSE = norm(data1-z(1,:)'-data2+z(2,:)') % % Figure 2 figure('Position',[300,150,650,300]); subplot(121); plot(x,base1,'k'); hold on; plot(x,z(1,:),'k:'); plot(x,base2,'k'); plot(x,z(2,:),'k:'); hold off; % legend('true baseline','estimated baseline') text(4350, 0.18, '(a)'); set(gca,'XDir','reverse'),set(gca,'XLim',[650 4000],'Ylim',[0,0.18]); xlabel('Wave number/cm^-^1'); ylabel('Absorbance'); ttt = 0:1:10; subplot(122); plot(ttt,a(1,:)','k');hold on; plot(ttt,a(2,:)','k:'); xlim([0 10]); ylim([0.9 1.1]); xlabel('Iterations'); ylabel('Relaxation factor'); text(-1.8, 1.1, '(b)') % % part of Fig.2(a), the local region 2500:-1:2000 xx = 2500:-1:2000; tt = 1500:2000; figure('Position',[300,150,620,300]); plot(xx,base1(tt),'k'); hold on; plot(xx,z(1,tt),'k:'); plot(xx,base2(tt),'k'); plot(xx,z(2,tt),'k:');hold off; set(gca,'XDir','reverse'),set(gca,'XLim',[2000 2500],'Ylim',[0.029,0.047]) set(gca,'XTick',[2000 2500]); set(gca,'YTick',[0.03 0.04 0.046]); %% ============================================================= %% ================ Asymmetric Least Squares =================== % grid search for the optimal parameters: lambda and mu for k = 1 : 2 data0 = data(k,:)'; base0 = base(k,:)'; % coarse tuning Err = []; pval = [0.1 0.01 0.001 0.0001]; muval = 10.^[1:10]; for p = pval for mu = muval bb = Aslsbase(data0, mu, p); err = norm( base0 - bb ); Err = [Err err]; end end pp = length(pval); mm = length(muval); ErrS = zeros(pp,mm); for i = 1: pp index = 1 : mm : pp*mm; ErrS(i,:) = Err(index(i): index(i)+mm-1); end mmin = min(min(ErrS)); [RI CI] = find(ErrS==mmin); po = pval(RI); muo = muval(CI); % for data1, p=0.0001,mu=1e6 % for data2, p=0.0001,mu=1e6 % fine tuning of mu, p is usually insensitive. Err2 = []; muval2 = [1e5: 1e5: 1e7]; for mu = muval2 bb1 = Aslsbase(data0, mu, po); err2 = norm(base0 - bb1); Err2 = [Err2 err2]; end mmin = min(min(Err2)); II = find(Err2==mmin); muoo = muval2(II); % for data1, mu=4e5 % for data1, mu=9e5 popt(k) = po; muopt(k) = muoo; end bb1 = Aslsbase(data1, muopt(1), popt(1)); bb2 = Aslsbase(data2, muopt(2), popt(2)); RMSE11 = norm(base1-bb1) % 0.0992 RMSE22 = norm(base2-bb2) % 0.0538 RMSE_2 = norm(data1-bb1-data2+bb2) % 0.0666 figure; plot(data1); hold on; plot(bb1,'r'); hold off; figure; plot(data2); hold on; plot(bb2,'r'); hold off;