General Functions
Peak shapes
The following functions are useful when generating peaks with various shapes. See the examples for using them during peak fitting for instance.
#
Spectra.gaussiennes — Method.
gaussiennes(amplitude::Array{Float64},centre::Array{Float64},hwhm::Array{Float64},x::Array{Float64};style::String = "None")
gaussiennes, written in the plural french form there, is a function that allows to build gaussian peaks. The gaussian function used there is:
y = amplitude x exp(-ln(2) x [(x-centre)/hwhm]^2 ) You can enter the amplitude, centre and half-width at half-maximum (hwhm) values as arrays of float 64 (even containing one float value), without specifying style. hwhm is proportional to the standard deviation sigma: hwhm= sqrt(2xln(2)) x sigma
that is used in a normal distribution (see function normal_dist).
INPUTS:
amplitude: Array{Float64} containing the peaks amplitudes; centre: Array{Float64} containing the peaks centres; hwhm: Array{Float64} containing the peaks half-width at middle heights (hwhm); x: Array{Float64} containing the x axis values;
OPTIONS:
style: ASCIIString = "None", see examples below.
OUTPUTS:
y_calc: Array{Float64} containing the calculated y values; y_peaks: Array{Float64} containing the y values of the different peaks.
Examples
To have four gaussian peaks centered at 800, 900, 1000 and 1100 cm-1 with hwhm of 50 cm-1 on a Raman spectrum, you will enter:
y_calc, y_peaks = gaussiennes([1.0,1.0,1.0,1.0], [800.0,900.0,1000.0,1100.0], [50.0,50.0,50.0,50.0], x)
and y_peaks will contain in 4 columns the 4 different y values of the peaks, and y_calc their sum (the total model). Now, if you want to calculate more complex models, such as for instance contructing how the Raman peaks of water vary with pressure, you might like to parametrize the variations of the peak parameters rather than just fitting each spectrum. This will provide more robust fits of the spectra, as you will fit them together, and will also force you to find the correct underlying mathematical assumption.
The gaussiennes function allows you to do that. If you specify style = "poly", you can enter arrays for the amplitudes, centres and half-widths at half-maximum (hwhm) of the peaks, with in each column the coefficients for the polynomial variations of this parameters. The second column of x will need to contain the second variable for those polynomial functions.
Let's say for instance that we have one peak at 900 cm-1 in a pure material. It's frequency seems to linearly shift with increasing the amount of hydrogen in this material, but it's intensity is non-linearly increasing, following a quadratic variation. It's width seems constant.
How to write that with gaussiennes? Well, first you need to construct a relevant x axis: first column contains the frequency, and the second one contains the chemical variable value. In our case, we want to model the peak between 800 and 1000 cm-1, for 1 wt% H. So we have an x array build like:
frequency = collect(800:1:1000) x = ones(length(frequency),2) x[:,1] = frequency[:] x[:,2] = 1.0
Ok, now lets build our y peaks:
amplitudes = [1.0 0.1 0.1] frequencies = [900.0 2.0] hwhm = 20.0 y_calc, y_peaks = gaussiennes(amplitudes, frequencies, hwhm, x)
This should provide you how the shape of the peak is as a function of both the frequency and the chemical composition there. If you want to go further, you might just want to stick gaussiennes in a loop, and play with creating various peaks with changing the chemical parameter in the x[:,2] column!
#
Spectra.lorentziennes — Method.
lorentziennes(amplitude::Array{Float64},centre::Array{Float64},hwhm::Array{Float64},x::Array{Float64};style::String = "None")
INPUTS:
amplitude: Array{Float64} containing the peaks amplitudes; centre: Array{Float64} containing the peaks centres; hwhm: Array{Float64} containing the peaks half-width at middle heights (hwhm); x: Array{Float64} containing the x axis values;
OPTIONS:
style: ASCIIString = "None", see examples in the gaussiennes documentation.
OUTPUTS:
y_calc: Array{Float64} containing the calculated y values; y_peaks: Array{Float64} containing the y values of the different peaks.
#
Spectra.pearson7 — Method.
pearson7(a1::Array{Float64},a2::Array{Float64},a3::Array{Float64},a4::Array{Float64},x::Array{Float64};style::String = "None")
a Pearson 7 peak with formula a1 ./ (1 + ((x-a2)./a3).^2 .* (2.0.^(1./a4) - 1.0))
INPUTS:
a1: Array{Float64} ; a2: Array{Float64} ; a3: Array{Float64} ; a4: Array{Float64} ; x: Array{Float64} containing the x axis values;
OPTIONS:
style: ASCIIString = "None", see examples in the gaussiennes documentation.
OUTPUTS:
y_calc: Array{Float64} containing the calculated y values; y_peaks: Array{Float64} containing the y values of the different peaks.
#
Spectra.pseudovoigts — Method.
pseudovoigts(amplitude::Array{Float64},centre::Array{Float64},hwhm::Array{Float64},lorentzian_fraction::Array{Float64},x::Array{Float64};style::String = "None")
A mixture of gaussian and lorentzian peaks.
INPUTS:
amplitude: Array{Float64} containing the peaks amplitudes; centre: Array{Float64} containing the peaks centres; hwhm: Array{Float64} containing the peaks half-width at middle heights (hwhm); lorentzian_fraction: Array{Float64}, containing the lorentzian fraction of the pseudovoigt function. Should be comprised between 0 and 1; x: Array{Float64} containing the x axis values;
OPTIONS:
style: ASCIIString = "None", see examples in the gaussiennes documentation.
OUTPUTS:
y_calc: Array{Float64} containing the calculated y values; y_peaks: Array{Float64} containing the y values of the different peaks.
#
Spectra.normal_dist — Method.
normal_dist(nd_amplitudes::Array{Float64},nd_centres::Array{Float64},nd_sigmas::Array{Float64},x::Array{Float64})
The real normal distribution / gaussian function
INPUTS:
amplitude: Array{Float64} containing the peaks amplitudes; centre: Array{Float64} containing the peaks centres; hwhm: Array{Float64} containing the peaks half-width at middle heights (hwhm); lorentzian_fraction: Array{Float64}, containing the lorentzian fraction of the pseudovoigt function. Should be comprised between 0 and 1; x: Array{Float64} containing the x axis values;
OPTIONS:
style: ASCIIString = "None", see examples in the gaussiennes documentation.
OUTPUTS:
y_calc: Array{Float64} containing the calculated y values; y_peaks: Array{Float64} containing the y values of the different peaks.
Peak measurement
#
Spectra.peakmeas — Method.
peakmeas(x::Array{Float64},y::Array{Float64};smoothing = "yes", filter = :SavitzkyGolay, M=5,N=2,ese_y=1.,y_smo_out=false)
The peakmeas function allows performing measurements of the position, width, intensity and centroïd of a dominant peak in a provided x-y signal.
It smooths the signal with a Savitzky-Golay filter prior to measuring the peak position, width and intensity. It is advised to check that the M and N values of the Savitzky-Golay filter are adequate for your problem before trusting the results from peakmeas. For that, just use the y_smo_out option.
half-width at half-maximum are calculated as the width of the peak at half its maximum intensity. This calculation is not affected by any asumption of peak symmetry (no fitting is done).
Centroïd is calculated as sum(y./sum(y).*x).
INPUTS:
x: Array{Float64}, the x values; y: Array{Float64}, the y values.
OPTIONS:
smoothing, String, triggers the smoothing of the spectrum if set to yes (default value); filter, Symbol, the filter that will be used. See the smooth function documentation; M=5, the M parameter for smoothing y with a Savitzky-Golay filter. See smooth function documentation; N=2, the M parameter for smoothing y with a Savitzky-Golay filter. See smooth function documentation; y_smo_out=false, the smoothed signal. Signal will be smoothed if set to true, using the SavitzkyGolayFilter function with the M and N values. y_smo output will also be provided.
OUTPUTS:
intensity: Float64, the intensity of the peak; position: Float64, the position of the peak; hwhm: Float64, the half-width at half-maximum of the peak; centroïd: Float64, the centroïd of the peak; if y_smo_out is set to true, then another output is provided: y_smo: Array{Float64}, the smoothed y signal.
Integration
Spectra.jl provides functions that allow one to integrate the area under a region of a spectrum, or to calculate the area under Gaussian, Lorentzian or other bands.
#
Spectra.trapz — Method.
trapz{Tx<:Number, Ty<:Number}(x::Vector{Tx}, y::Vector{Ty})
Trapezoidal integration.
INPUTS:
x: Vector{Float64} containing the x values; y: Vector{Float64} containing the y values.
OUTPUTS:
area: Vector{Float64}, the trapezoidal integration value.
This function is particularly helpful to calculate the area under a portion of a spectrum, and can be used for various purposes (normalisation, area comparison, etc.).
#
Spectra.bandarea — Method.
bandarea(Amplitude::Array{Float64},HWHM::Array{Float64}; peak_shape = "Gaussian", error_switch = "no", eseAmplitude::Array{Float64} = [0.0], eseHWHM::Array{Float64} = [0.0])
This function replaces the function gaussianarea in the version <0.1.9 of Spectra.jl. It allows to calculate the area under a specific band, with different shapes. For now, only Gaussian bands are supported, but other band shapes will be added soon. (This explains why gaussianarea is deprecated in favor of a more generic function)
gaussianarea allows to calculate the area under a gaussian peak from its half-width at half maximum (hwhm) and its amplitude, with the possibility of calculating the error based on the inputs of the errors on hwhm and amplitude. Call it as:
area, esearea = band(Amplitude,HWHM; peak_shape = "Gaussian", error_switch = "no", eseAmplitude = [0.0], eseHWHM = [0.0])
INPUTS:
Amplitude: Array{Float64}, contains the amplitudes (intensity) of the band(s); HWHM: Array{Float64}, contains the half width at half maximum of the peaks;
OPTIONS
peak_shape: String, indicates the shape of the component. Only "Gaussian" is supported for now; error_switch: String, should be "yes" or "no". If "yes", the arrays containing the errors affecting the band amplitude and widhts should be provided in eseAmplitude and eseHWHM (see below); eseAmplitude: Array{Float64}, an array that contains the errors affecting Amplitude; eseHWHM: Array{Float64}, an array that contains the errors affecting HWHM;
OUTPUTS:
area: Array{Float64}, an array that contains the areas; if error_switch is set to "yes", then a second output is provided: esearea: Array{Float64}, an array that contains the propagated errors affecting the areas calculations.
Polynomials
#
Spectra.poly — Method.
poly(p::Vector{Float64},x::Array{Float64})
This function just allows to build a polynomial curve.
INPUTS:
p: Vector{Float64}, containing the polynomial parameters. For a linear curve, p = [1.0,1.0], for a second order polynomial, p = [1.0,1.0,1.0], etc.; x: Array{Float64}, containing the x values for calculation.
Output:
y: Array{Float64}, containing the result of calculation.
#
Spectra.polyfit — Method.
polyfit(x::Array{Float64}, y::Array{Float64}, n::Int64)
Fit a polynom of degree n to x-y data. Code from https://rosettacode.org/wiki/Polynomial_regression#Julia
INPUTS:
x: Array{Float64}, the x values; y: Array{Float64}, the y values; n: Int64, the polynomial degree: 0 for a constant up to as high as wanted.
OUTPUTS:
coefs: Array{Float64}, containing the polynomial coefficients.
Splines
Not all the splines packages provide the same performances for data smoothing and interpolation. By experience, the Dierckx spline package ("Dspline" option in the baseline() function) provides a good starting point, but is not as usefull as other spline packages.
The csaps function of Matlab uses the SMOOTH Fortran library, and provides better smoothing capabilities for noisy data. Similarly, the GCVSPL Fortran package from Woltring (1986) also provides a very robust way to smooth and interpolate noisy data.
Starting from Spectra v0.3.4, the gcvspline Python module (https://github.com/charlesll/gcvspline) is used in the smooth and baseline function. It behaves exactly as the previous wrapping of GCVSPL.f in Julia, such that this should be transparent to users.