import numpy as np import matplotlib.pyplot as plt from numpy.linalg import inv from pytictoc import TicToc t = TicToc() t.tic() # Define a Gaussian peak function def gaussian(x, center, amplitude, width): # This generates a single peak with a gaussian lineshape, # with a specific center, amplitude, and width. # width is the FWHM (full width at half maximum) return amplitude*np.exp(-((x-center)/(0.60056120439323*width))**2) # Create a set of spectra G representing the components of the mixture NumPoints=1000 X=np.linspace(200,800,NumPoints) center1=450 G1=gaussian(X,center1,1,100) center2=500 G2=gaussian(X,center2,1,120) center3=550 G3=gaussian(X,center3,1,130) E=np.c_[G1,G2,G3] # compose the separate G vectors into a matrix E noiseamp=.01 # Specify the amplitudes of the components in the mixture amp1=1 # Amplitude (height) of peak 1 amp2=.5 # Amplitude (height) of peak 2 amp3=.3 # Amplitude (height) of peak 2 # Compute the simulated noisy mixture spectrum A=amp1*G1+amp2*G2+amp3*G3+noiseamp*np.random.randn(len(X)) t.tic() # Start clock for execution timimg # Apply normal equation to compute concentration vector C C = inv(E.T.dot(E)).dot(E.T).dot(A) plt.plot(X,A) plt.plot(X,C[0]*G1) plt.plot(X,C[1]*G2) plt.plot(X,C[2]*G3) PercentDifference1=100*(amp1-C[0])/amp1 PercentDifference2=100*(amp2-C[1])/amp2 PercentDifference3=100*(amp3-C[2])/amp3 plt.text(center1,C[0],PercentDifference1, fontsize=22) plt.text(center2,C[1],PercentDifference2, fontsize=22) plt.text(center3,C[2],PercentDifference3, fontsize=22) plt.title('Python: Percent error of peak heights measured by Classical Least Squares') t.toc() # stop clock and print elapsed time print(C) # Print measured amplitudes