WO2014125819A1 - 解析装置 - Google Patents

解析装置 Download PDF

Info

Publication number
WO2014125819A1
WO2014125819A1 PCT/JP2014/000730 JP2014000730W WO2014125819A1 WO 2014125819 A1 WO2014125819 A1 WO 2014125819A1 JP 2014000730 W JP2014000730 W JP 2014000730W WO 2014125819 A1 WO2014125819 A1 WO 2014125819A1
Authority
WO
WIPO (PCT)
Prior art keywords
discrete
spectrum
coupling coefficient
spline
negative
Prior art date
Application number
PCT/JP2014/000730
Other languages
English (en)
French (fr)
Inventor
賢 鎌田
Original Assignee
アトナープ株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by アトナープ株式会社 filed Critical アトナープ株式会社
Publication of WO2014125819A1 publication Critical patent/WO2014125819A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N27/00—Investigating or analysing materials by the use of electric, electrochemical, or magnetic means
    • G01N27/62—Investigating or analysing materials by the use of electric, electrochemical, or magnetic means by investigating the ionisation of gases, e.g. aerosols; by investigating electric discharges, e.g. emission of cathode
    • G01N27/622—Ion mobility spectrometry
    • G01N27/624—Differential mobility spectrometry [DMS]; Field asymmetric-waveform ion mobility spectrometry [FAIMS]

Definitions

  • the present invention relates to an apparatus and method for analyzing a spectrum.
  • a unit for analyzing a sample disclosed in International Publication No. WO2012 / 056709 is obtained by supplying a sample to an ion mobility sensor that measures the ionic strength of a chemical substance passing through an electric field controlled by at least two parameters.
  • a curve fitting method using Gaussian characteristics is known as an analysis method of measurement data that may have a plurality of peaks. There is a further need for a method that can analyze the presence of multiple peaks at high speed.
  • One aspect of the present invention is a spectral sampling value obtained from a sensor and defined along the first coordinate axis of the spectrum, having a positive value at a finite number of sampling points, and 0 (zero,
  • a first arithmetic unit that calculates an inner product with a local non-negative discrete function that is zero) and represents a spectrum by linear combination approximation of a plurality of local non-negative discrete functions, and a plurality of spectra based on the coupling coefficient
  • a separation unit that decomposes into a distribution corresponding to a local non-negative discrete function.
  • the first arithmetic unit fixes the coupling coefficient to 0 (zero) when the obtained coupling coefficient is negative.
  • Spectra obtained from sensors that capture quantities related to physical phenomena are often a combination of several non-negative distributions.
  • spectra obtained from sensors that output information about the presence of chemicals such as ion mobility sensors, mass analyzers, gas chromatography, Raman, and other spectroscopic analyzers are non-negative in relation to the presence of some chemicals. It can be interpreted as a combination of distributions. Therefore, by preparing an appropriate local non-negative discrete function and calculating the inner product and calculating the coupling coefficient, the spectrum can be approximated at high speed without curve fitting to infinity, and the distribution contained in the spectrum can be analyzed .
  • the coupling coefficient of a plurality of local non-negative discrete functions does not become negative. For this reason, when the coupling coefficient is negative, by fixing the coupling coefficient to 0 (zero), the spectrum can be approximated by coupling with a distribution that more closely matches the physical phenomenon to be measured.
  • the local non-negative discrete function preferably includes a discrete B-spline that is an integer m-fold convolution accumulation of a rectangular discrete function having a value of 1 at consecutive integer n sampling points and a value of 0 at other sampling points. .
  • a discrete B-spline that is an integer m-fold convolution accumulation of a rectangular discrete function having a value of 1 at consecutive integer n sampling points and a value of 0 at other sampling points.
  • the first arithmetic unit includes an inner product arithmetic unit that calculates an inner product of the spectrum and the discrete B-spline, and the separation unit decomposes the spectrum into a plurality of B-spline distributions based on the coupling coefficient. It is desirable to include.
  • the first arithmetic unit includes a first digital filter unit including a transfer function of the following expression (A) and a second digital filter unit including a transfer function of the following expression (B).
  • a discrete B-spline can calculate the inner product at high speed using these transfer functions.
  • the analysis apparatus obtains a coupling coefficient candidate in which at least one of n and m is different by the first arithmetic unit, evaluates an error between the linear combination approximation of the discrete B-spline by the coupling coefficient candidate and the spectrum, It is desirable to further have a second arithmetic unit for obtaining the coupling coefficient.
  • An integer n and / or m suitable for approximating the spectrum can be selected, and peak separation performance can be improved.
  • the first arithmetic unit may include a unit that selects a discrete B-spline having a large positive inner product with respect to the spectrum and sequentially obtains a coupling coefficient when the discrete B-spline is expressed by linear combination approximation.
  • SPmax is the number of samplings calculated in the inner product calculation unit. There is almost no meaning to use a discrete B-spline having convolution accumulation exceeding the maximum value of the sampling points, and the operation time can be shortened by using integers m and n satisfying the equation (C). n ⁇ m ⁇ SPmax (C)
  • An example of the spectrum to be analyzed is the output of the ion mobility sensor, and an example of the first coordinate axis is the compensation voltage.
  • Another aspect of the present invention is a method for analyzing a spectrum obtained from a sensor.
  • the analysis method includes the following steps. 1. Compute the inner product of a spectrum sampling value and a local non-negative discrete function to obtain a coupling coefficient for expressing the spectrum by linear combination approximation of multiple local non-negative discrete functions. 2. Decomposing the spectrum into distributions corresponding to multiple local non-negative discrete functions based on the coupling coefficient. Obtaining the coupling coefficient includes fixing the coupling coefficient to 0 if the obtained coupling coefficient is negative.
  • the local non-negative discrete function includes a discrete B-spline that is an integer m-fold convolution accumulation of a rectangular discrete function having a value of 1 at n consecutive integer sampling points and a value of 0 (zero) at other sampling points. It is desirable.
  • determining the coupling coefficient includes computing an inner product of the spectrum and the discrete B-spline, and separating includes decomposing the spectrum into a plurality of B-spline distributions based on the coupling coefficient. Is desirable.
  • calculating the inner product preferably includes calculating the transfer functions of the above formulas (A) and (B).
  • obtaining the coupling coefficient includes evaluating an error between the linear combination approximation of the discrete B-spline and the spectrum by a coupling coefficient candidate in which at least one of the integer n and the integer m is different, and obtaining the coupling coefficient. Is desirable.
  • a program is a method in which a processor calculates an inner product of a spectrum sampling value and a local non-negative discrete function to obtain a coupling coefficient when the spectrum is expressed by a linear combination approximation of a plurality of local non-negative discrete functions; Decomposing the spectrum into a distribution corresponding to a plurality of local non-negative discrete functions based on the coupling coefficient, and determining the coupling coefficient is to fix the coupling coefficient to 0 when the obtained coupling coefficient is negative including.
  • the block diagram which shows schematic structure of an analyzer.
  • the figure which shows a mode that an inner product is calculated using a digital filter.
  • the figure which shows a mode that a turning coupling approximation is calculated from a coupling coefficient using a digital filter.
  • the flowchart which shows the outline of an analysis method.
  • FIG. 1 shows an example of an analysis device for analyzing the output (measurement data, IMS data) of the ion mobility sensor 1.
  • This analysis device 10 may analyze IMS data 19 input directly from the ion mobility sensor 1 or analyze IMS data 19 stored in an appropriate storage (recording medium) 2. There may be.
  • the analysis apparatus 10 includes an arithmetic unit (first arithmetic unit) 11 that obtains a coupling coefficient for approximating the IMS data 19 by discrete B-splines, and an evaluation unit (first arithmetic unit that evaluates an approximation situation based on the obtained coupling coefficient. Selected as a linear combination of discrete B-splines approximating IMS data 19 by the evaluation unit 12 and a library 13 in which setting values used by these units 11 and 12 are stored. Peak separation unit 14 for separating and recognizing peaks included in the B-spline approximation, and classifying the peaks based on attributes such as position, height, shape, and trend of the separated peaks (distribution). An estimation unit 15 for estimating the measurement object and a database including data on the separated peaks And a 16.
  • the analysis apparatus 10 is typically configured using computer resources including a processor (CPU) and a memory, and is recorded as software (program or program product) on an appropriate recording medium or a computer network such as the Internet. Provided through the wear.
  • the analysis device 10 may be provided by an ASIC, an LSI, a reconfigurable semiconductor circuit device (processor), an optical circuit device, or the like.
  • an asymmetric field ion mobility spectrometer FIMS
  • DMS differential electric mobility spectrometer
  • An example of the FAIMS sensor 1 is an ion mobility analyzer described in International Publication WO2006 / 013396 or International Publication WO2007 / 014303. These are monolithic or micromachine type compact mass spectrometers, which are sufficiently portable.
  • the FAIMS sensor 1 includes a drift chamber that moves the ionized measurement object while being influenced by an electric field, and a detector that detects the ionized measurement object (charge of the measurement object) that has passed through the drift chamber.
  • the software-controlled electric field generated by the voltage Vd and the voltage Vc fluctuates positively and negatively at a specific period, and the chemical substance to be detected is filtered by the filtering effect of the electric field, For example, it collides with the detector at the msec level and is measured as ion intensity (ion current) Ic.
  • Vd voltage Dispersion Voltage, dispersion voltage, electric field voltage (Vd)
  • Vc voltage Compensation Voltage, compensation voltage
  • Data (measurement data, IMS data) 19 obtained from the ion mobility sensor is at least three-dimensional data including the value of the ion current.
  • the value of the ionic current may form a plurality of peaks when the Vc voltage is changed. Therefore, there is a possibility that the chemical composition (chemical substance contained in the sample) of the sample detected by the ion mobility sensor 1 can be specified based on the peak contained in the IMS data 19.
  • FIG. 2 shows the m-floor B-spline.
  • FIG. 2A shows a B-spline obtained by convolving and integrating a rectangular function once, and is a first-order B-spline.
  • 2B shows a B-spline on the second floor
  • FIG. 2C shows a B-spline on the third floor
  • FIG. 2D shows a B-spline on the fourth floor.
  • FIG. 3 shows a state in which a discrete B-spline is defined by m-order convolution accumulation of n sample sequences (sampling points) obtained from a rectangular function.
  • An m-th order discrete B-spline is obtained as an output of a digital filter defined by a transfer function by performing z-transform as shown in the following equation (0).
  • Ion mobility analysis is a method for discriminating analysis of chemical molecules. Because it can analyze trace amounts of chemical substances, it can be used to analyze odors and scents and detect poisonous substances, drugs, and explosives. The analysis procedure is roughly divided into physical processing and computational processing.
  • FIG. 4 (a) and 4 (b) schematically show physical processing in the ion mobility sensor 1.
  • FIG. The ionized chemical molecules move inside the ion mobility sensor 1 and are caught by the positive electrode 1a or the negative electrode 1b.
  • a typical ionization method is electron beam irradiation from a radioactive element, and may be a method such as UV irradiation or corona discharge. While the molecules ionized inside the ion mobility sensor 1 are moving, as described above, the molecules move zigzag by the dispersion voltage (Vd voltage) and the compensation voltage (Vc voltage), and the conditions of the Vd voltage and the Vc voltage are satisfied.
  • Vd voltage dispersion voltage
  • Vc voltage compensation voltage
  • IMS data 19 of the ion mobility sensor 1.
  • An example of the IMS data 19 is a spectrum (a spectrum distribution, a spectrum, a waveform distribution, and an output distribution) in which the ion current Ic varies due to a change in the Vc voltage.
  • the ionized molecules are influenced by the Vd voltage and the Vc voltage in proportion to their charge / mass, and collide with air molecules on the way and proceed while diffusing. For this reason, on average, it proceeds at a speed peculiar to the types of ions 3a to 3c, and varies with time to reach an electrode, for example, the positive electrode 1a.
  • the variation is assumed to follow a normal distribution. This is because if the position change caused by each collision varies according to the same distribution, the distribution obtained through infinite collisions becomes a Gaussian distribution.
  • FIG. 4C shows an example of the IMS data 19.
  • This waveform is a spectrum (spectral distribution) of the ion current obtained when the Vc voltage is changed while the Vd voltage is constant.
  • This current waveform is called a profile, and the first coordinate axis is the Vc voltage (compensation voltage). Since one Gaussian distribution corresponds to one kind of ion, the profile is modeled as an unloaded multiple sum of a plurality of Gaussian distributions 19a to 19c as shown in FIG. 4 (d).
  • the conventional calculation process is aimed at identifying each distribution constituting it from a given profile.
  • the type of ion molecule is determined by the average value of each distribution, and the amount is determined by height. Since the number, average, and variance of the Gaussian distribution are unknown, it is a common technique to search for these parameters by the steepest descent method. The search is a sequential process and requires many repeated operations until convergence. For this reason, the current calculation processing takes a lot of time even if a high-performance CPU is used.
  • processing time is reduced, physical phenomena rarely extend from negative infinity to positive infinity, and local non-negative discrete functions are more suitable for approximating distributions based on physical phenomena. There are many cases.
  • an integer m-order B-spline representing a variation in velocity obtained through a finite number of (m) collisions according to a uniform distribution can be used.
  • B-splines Although it is not the original use of B-splines as a curve representation, if the integer m is increased, the B-splines approach a Gaussian distribution. For example, according to the principle of ion mobility analysis, the shape of the distribution cannot also be a perfect normal distribution with an infinite base. It is not unnatural to use B-splines instead of the normal distribution.
  • a discrete B-spline obtained by replacing a uniform continuous distribution with a uniform discrete distribution sampled at n points can be generated by only addition and subtraction (T. Saramaki, Y. Neuvo and S. K. Mitra, Design of computationally efficient interpolated FIR filters, IEEE Trans. Circuits & Syst., 35 (1), 70-88, 1988).
  • the discrete B-spline converges to the original B-spline by increasing the number of sample points n, and the waveform can be approximated to a discrete B-spline by a high-speed digital filter (K. Ichige and M.
  • the high-speed digital filter is used for n fixed at various values, assuming that the spread of ions reaching the same period is the same. Approximate brute force execution in parallel. The case where a high-precision approximation with a small number of discrete B-splines is obtained is taken as the analysis result.
  • One method that can be used to determine approximate accuracy is sparse approximation. (For sparse approximation, M. Vetterli, P. Marziliano and T. Blu, Sampling signals with finite rate of innovation, IEEE Trans. Signal Process., SP-50, 1417-1428, 2002).
  • the m-th order B-spline is defined as an integer m-fold convolution integral of the rectangular function shown in FIG.
  • the central limit theorem shows that this converges to a normal distribution at the limit where the integer m is infinite.
  • a discrete B-spline that is an m-fold convolution accumulation of n sample value sequences obtained from a rectangular function converges to a B-spline in the limit that makes the integer n infinite.
  • this discrete B-spline is used as a distribution.
  • a discrete B-spline is defined as an integer m-fold convolution integral of a rectangular discrete function having a value of 1 at an integer n sampling points and a value of 0 at other sampling points.
  • the integer m and the integer n are integers of 1 or more and have the following relationship with the number of samples (sampling number) SPmax used for analysis. n ⁇ m ⁇ SPmax (C)
  • the discrete B-spline on the first floor is defined by the following equation (1).
  • the z conversion is expressed by the following equation (2).
  • the number of samples (number of samples) is not infinite but finite, that is, the maximum value SPmax, and the range of the variable k is (0, 1, 2,... M ⁇ n ⁇ 1). .
  • the inner product other than the sampling value is 0, an infinite sum expression may be used in this specification.
  • this inner product can be calculated as the following equation (8).
  • This operation corresponds to the following expression (11) when z conversion is performed.
  • the inner product calculation unit 11d of the calculation unit 11 of the analysis apparatus 10 includes a first digital filter unit 11a including the transfer function of Expression (A) and a second digital filter unit 11b including the transfer function of Expression (B). Including. Although overflow occurs during the calculation in these digital filter units 11a and 11b, these filter units 11a and 11b have a function of performing a “2” complement operation having the number of bits that can accommodate the final result of the calculation. And the final result is correct.
  • z conversion of ( ⁇ bm [k]) is Bm (z ⁇ 1 ). Accordingly, the z-transformed expression corresponds to the following expression (16).
  • This inner product is obtained by applying the output obtained by inputting the profile p [r] to a digital filter having a transfer function “((1-z ⁇ n ) / (1 ⁇ z ⁇ 1 )) m ” ((n ⁇ 1) Obtained by sampling at m + r). Similar to the above, this processing can be executed by the first digital filter unit 11a and the second digital filter unit 11b which are composed of two parts.
  • the inner product of the profile p [k] and the discrete B-spline bm [k-r] is calculated by a sum-of-products operation as shown in the above equation (14) in a normal way of thinking.
  • the equations (15) and (16) it can be understood that the calculation can be performed by the digital filter units 11a and 11b configured only by addition and subtraction. Therefore, in the inner product calculation unit 11d, the calculation configuration can be simplified, and the inner product calculation can be processed at high speed.
  • the profile p [r] is input to the digital filter units 11a and 11b of the inner product calculation unit 11d, and the inner product ⁇ p [•], bm [• -r]> at ((n ⁇ 1) m + r) is calculated. It shows how to do.
  • the m-th order accumulation calculated by the first digital filter unit 11a which is the first half does not depend on the integer n, and the integer n is included only in the m-th order difference calculated by the second digital filter unit 11b which is the second half. Yes. Therefore, the m-th order accumulation may be calculated only once for each profile, and the second digital filter unit 11b only calculates the m-th order difference while changing the integer n. It is possible to extract the inner product of the profile and the profile. The computation required for calculating one inner product is almost m subtractions, and the processing time can be shortened.
  • the first arithmetic unit 11 further includes a least square approximation unit (least square approximation unit) 11c that performs non-negative coefficient least square approximation.
  • a profile p [k] is expressed by a function q [k] obtained by linearly combining K discrete B-splines arranged at an interval “1 (one)” with a coupling coefficient c [l].
  • the linear combination approximation q [k] is expressed by the following equation (17), and the square error is calculated by the following equation (18).
  • the coefficient (coupling coefficient) c [l] that minimizes the square error (square error) E is determined by solving the following linear equation (19) according to a normal least square approximation method.
  • the spectrum 19 to be analyzed is the output of the ion mobility sensor 1. Therefore, it can be assumed that the coefficient is proportional to the number of ions. For this reason, the constraint that the coefficient c [l] is non-negative is imposed. Even if the output is other than the output of the ion mobility sensor 1, the coefficient c [l] is non-negative if it is a physical quantity such as a spectrum (distribution) indicating the presence of a chemical substance.
  • the least square approximation unit 11c fixes negative coefficients among the coefficients obtained by the normal least square approximation to determine the coefficients under this condition, and sets the corresponding discrete B-spline to 0 (zero).
  • FIG. 6 shows the linear combination approximation q [k] obtained for approximating the profile from the coupling coefficient c [k] obtained above using the first digital filter unit 11a and the second filter unit 11b. It shows how to calculate. Since the coupling coefficient is limited to non-negative, the number of coupling coefficients used in the approximation linear combination is limited, and the peaks included in the profile p [k] can be substantially separated.
  • the analysis unit 10 further includes an evaluation unit (second arithmetic unit) 12 that evaluates the linear combination approximation q [k] obtained from the coupling coefficient obtained by the first arithmetic unit 10.
  • the evaluation unit 12 obtains in advance coupling coefficients as candidate coupling coefficients for various integers n and / or integers m, and evaluates the resulting linear coupling approximation q [k] against the profile p [k].
  • the optimum coupling coefficient from among the coupling coefficient candidates and thereby the linear coupling approximation is output. That is, the evaluation unit 12 evaluates the brute force of the linear combination approximation of the discrete B-spline to the profile for various assumed m and n. Since the inner product calculation of the discrete B-spline can be processed in a short time as described above, it is practically possible to evaluate the brute force.
  • the evaluation unit 12 includes a first evaluation unit 12a that evaluates a square error and a second evaluation unit 12b that performs a sparsity evaluation.
  • the first evaluation unit 12a calculates a square error. Calculation of the linear combination approximation q [k] from the coupling coefficient c [k] is performed in a short time by the first digital filter unit 11a and the second digital filter unit 11b as shown in FIG. Can be executed. The square error between the linear combination approximation q [k] obtained in this way and the given profile p [k] is obtained by accumulating (p [k] ⁇ q [k]) 2 .
  • the first evaluation unit 12a calculates the mean square error E1 by the following equation (20).
  • the second evaluation unit 12b sets the case where the square error E1 is sufficiently reduced by the first evaluation unit 12a from among the coupling coefficient candidates obtained for various integers n as candidates for analysis results.
  • the second evaluation unit 12b searches for a combination of coupling coefficients that provides a sparse good approximation result from among the coupling coefficient candidates. This is based on the premise that it is an essential condition for a correct analysis that the profile p can be approximated by a small number of discrete B-splines among the obtained candidates.
  • the second evaluation unit 12b selects and outputs a case where the evaluation amount E2 is maximum as an analysis result.
  • FIGS. 7A to 7J show an example of the profile p [k] in the range of a local window of sample values of length 128 points for various integers n, where the integer m is fixed to 4.
  • the approximate result is shown.
  • the profile is the output (ion current, IMS data, spectrum) 19 of the ion mobility sensor 1, and the coordinate axis (first coordinate axis) is the compensation voltage Vc.
  • a solid line indicates a profile
  • a broken line indicates an approximate curve (linear coupling approximation) q [k]
  • a one-dot chain line indicates a discrete B-spline constituting the approximate curve q [k].
  • the approximate error E1 obtained by the first evaluation unit 12a is 10% or less when the integer n is 6 to 12, and larger than 15% in other cases. Therefore, the second evaluation unit 12b calculates the evaluation amount E2 for the approximate curve (linear combination approximation) q [k] in the former case, that is, when the integer n is 6 to 12. When the integer n is 6 to 12, the evaluation amount E2 of the sparsity is maximized when the integer n is 12. This case is adopted as a sparse good approximation, and the second evaluation unit 12b outputs. In this example, a result that two peaks are included in the profile p is output.
  • FIG. 8 is an example of an approximation result obtained by sampling one profile 19 in its entirety.
  • the optimal integer n may change over the course of the first coordinate axis.
  • the optimum integer n can be estimated for each window of the local compensation voltage Vc value.
  • the optimum values of the integer n for the left peak and the right peak are 12 and 14, respectively.
  • the peak separation unit 14 of the analysis apparatus 10 updates the peak database 16 with the peaks separated (detected peaks) based on the above results.
  • the peak classification and estimation unit 15 determines the peak type using the evaluation rules provided by the peak classification rule table and the peak type (type) rule table of the peak database 16.
  • the evaluation rules include the theoretical relationship between the peak characteristics of the reactant (RIP) and the peak characteristics of the monomer, dimer and trimer derived from the chemical substance.
  • the classification and estimation unit 15 includes a function of estimating a chemical substance contained in the sample gas measured by the ion mobility sensor 1.
  • the separation and estimation unit 15 includes a selection function for selecting a chemical candidate included in the sample gas based on the identified peak indicating at least one of a monomer, a dimer, and a trimer, and the identified peak and the third It includes a selection function that obtains a correlation with the parameter concentration and selects a chemical substance candidate contained in the sample gas.
  • FIG. 9 is a flowchart showing an outline of the analysis method in the analysis unit 10.
  • the analysis method (analysis process) is recorded on a suitable recording medium as a program (program product) for controlling a computer or a processor constituting the analysis unit 10, or is provided via a computer network and mounted on the analysis unit 10. .
  • the analysis process may be executed locally, may be executed on a cloud or on a server via a computer network, and may be processed in a distributed manner between the local and the cloud.
  • a local non-negative discrete function is set in step 51.
  • a discrete B-spline (formula (3)) will be used for explanation, but the local non-negative discrete function includes a raised cosine (Raised Cosine) and a spheroidal wavefunction (Prolate Spheroidal Wave Function). Any function that has a positive value at a finite number of sampling points and 0 (zero, zero) at other sampling points may be used.
  • step 52 the spectrum data (IMS data) 19 which is the output of the sensor 1 is acquired, and the profile p [•] is set.
  • step 53 a digital filter corresponding to the expressions (A) and (B) is prepared to calculate the inner product of the profile p [•] and the discrete B-spline bm [• -r], and the inner product is calculated at high speed.
  • the discrete B-spline bm [•] is an integer m-fold convolution accumulation of a rectangular discrete function having a value of 1 at integer n sampling points and a value of 0 at other sampling points.
  • a non-negative coupling coefficient is obtained by least square approximation.
  • the profile p [k] is a function obtained by linearly combining K discrete B-splines arranged at an interval “1” with a coupling coefficient c [l]. Calculate the square error when approximated by q [k].
  • the coefficient c [l] of the linear combination of the non-negative local discrete function is also made non-negative. Therefore, the negative coefficient first obtained by the least square approximation is fixed to 0 (zero), the corresponding discrete B-spline is excluded, and the least square approximation is performed so that all coefficients are non-negative. Repeat until.
  • step 55 the processing of steps 53 and 54 is performed for various integers n and m set in advance.
  • the integers n and m are selected within the range of the above formula (C).
  • the integers n and m may be limited to those suitable for approximation by linear combination depending on the measurement object of the sensor, the type and characteristics of the sensor, and the processing speed can be further improved.
  • the evaluation unit 12 evaluates an error from the profile p [k] from the linear combination approximation q [k] obtained as described above, and obtains a coupling coefficient.
  • the mean square error E1 (formula (20)) is small and the sparse approximation evaluation amount E2 (formula (23)) is large can be used.
  • the evaluation amount E0 is a combination of the mean square error E1 and the evaluation amount E2 of sparse approximation, and the smaller the evaluation amount E0, the better the analysis result.
  • the peak separation unit 14 corresponds the profile (spectrum) p [k] to the discrete B-spline of each coupling coefficient based on each coupling coefficient of the linear coupling approximation q [k] with good error evaluation result. Separate (decompose) into B-spline distribution.
  • the estimation unit 15 classifies the separated peaks and estimates the measurement object (chemical substance, combination thereof, binding, etc.) of the IMS data 19.
  • FIG. 10 shows a different example of step 54 for obtaining the non-negative coefficient by the method of least squares in the least square approximation unit 11c.
  • this example analysis method
  • not all of the discrete B-splines bm [ ⁇ -r] are used for approximation from the beginning, but selected from those having a large positive inner product with respect to the profile p [k] to be approximated.
  • Including the method of adding sequentially As a result, the processing time can be shortened and the risk of numerical instability can be avoided.
  • step 61 initial setting for repeated processing is performed.
  • step 62 an integer r that gives the maximum ⁇ p [•], bm [• ⁇ r]> is selected.
  • the integer r selected first is r1.
  • step 63 a coupling coefficient c [r1] of the least square approximation of p1 [ ⁇ ] by the bm [ ⁇ ⁇ r1] is obtained.
  • step 64 if this coefficient c [r1] is non-negative, a linear combination approximation q [•] is obtained. On the other hand, if the coefficient c [r1] obtained first is negative, it is determined that the only coefficient c [r1] is 0. Therefore, in step 67, the linear combination approximation is not obtained. For this reason, the calculation is terminated (step 68).
  • the mean square error E1 is calculated in step 69, and if the mean square error E1 is equal to or less than the threshold value, a linear combination approximation q is obtained, and the calculation ends (step 75). . If it is determined that a sufficiently accurate linear combination approximation q has been obtained using the given integers n and m, then it may proceed to error estimation (step 56) or peak separation (step 57), and other integers n and When obtaining the linear combination approximation q for m, the process is repeated while changing those values.
  • step 63 the least square approximation coupling coefficients c [r1] and c [r2] of the profile p [•] with the discrete B-splines bm [• -r1] ”and bm [• -r2] are obtained. If any of these coefficients are negative in step 64, they are fixed at 0 in step 65, and discrete B-splines corresponding to those coefficients are excluded in step 66. Returning to step 63, the least square approximation is repeated until the remaining coefficients, in this case, the two coefficients are non-negative, to obtain the linear combination approximation q [•]. If the two coefficients are both 0, the calculation ends in step 68, and if the square error (square error) E1 becomes sufficiently small, the calculation ends in step 75.
  • step 62 r which gives the next maximum inner product is selected.
  • Coupling coefficients c [r1] to c [r (i + 1)] of least square approximation of profile residuals pi [ ⁇ ] by B-splines bm [ ⁇ ⁇ r1] to bm [ ⁇ ⁇ r (i + 1)] Is obtained (step 63).
  • step 64 and 65 If any of these coefficients are negative, they are fixed to 0 (steps 64 and 65), the corresponding discrete B-splines are excluded (step 66), and the least squares approximation is expressed as i + 1 Repeat until all the coefficients (remaining coefficients) are non-negative to obtain the linear combination approximation q [•]. If all the coefficients are 0 or the square error E1 is sufficiently small, the calculation is terminated (steps 67, 68, 69, 75).
  • the linear combination approximation q [•] is set as the (i + 1) -th linear combination approximation q (i + 1) [•], and the linear combination approximation q (i + 1) [•] of the profile p [•]
  • the present invention has been described with an example in which the calculation mechanism related to the calculation method of the discrete B-spline is first applied to the ion mobility analysis. It is as explained above.
  • the arrangement interval of the discrete B-splines does not need to be fixed to an integer n, and the finest interval “1” is appropriate in the above-described application scene.
  • numerical instability may be caused by performing the least square approximation using an almost linearly dependent function system.
  • the numerical examples tried so far do not significantly affect the final processing result.
  • a discrete B-spline which is a good alternative to the Gaussian distribution
  • a system using discrete B-splines can be provided as a system that can be easily realized by hardware processing such as a processor.
  • One factor that is suitable for hardware implementation is that the inner product of a given waveform and a discrete B-spline and the load sum of the discrete B-spline can be calculated with a very simple digital filter.
  • the profile waveform is summed with the sum of the loads of discrete B-splines arranged at fine intervals, and the least square approximation is performed under the constraint that the load is non-negative. Detected when load increases. Since discrete B-splines arranged at fine intervals are almost linearly dependent, their least-square approximation is numerically unstable, but they are discretely given a negative load in the numerical calculation process. Since the B-spline is removed, the approximate result can be finally obtained stably. Since one-time approximation can be performed at high speed, discrete B-splines of various widths can be tried, and a case where a good sparse approximation using a large and small load is obtained can be adopted as an analysis result.
  • the evaluation unit 12 may set the single evaluation amount E0 for evaluating the goodness of analysis by combining the mean square error E1 and the evaluation amount E2 of the sparsity, and is applied to various current waveforms. Numerical stability, analytical performance, and processing speed may be evaluated. Further, the spectrum (profile) to be analyzed is not limited to the output of the ion mobility sensor 1, and the present invention is the output of other types of sensors such as a mass analyzer, Raman spectroscopic analysis, and infrared spectroscopic analysis. It can also be applied to.
  • the analysis target of the analysis apparatus 10 is a measurement result of a phenomenon estimated to follow a normal distribution such as a measurement result of a phenomenon involving a finite number of collisions according to a uniform distribution, and a measurement result of a phenomenon such as spectroscopy or absorption. Also good. Further, the first coordinate axis of the spectrum is not limited to the compensation voltage, and may be another physical quantity, time, or the like.

Landscapes

  • Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Chemical & Material Sciences (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Electrochemistry (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

第１の座標軸に沿って得られるスペクトルのサンプリング値と、連続する整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０となる矩形離散関数の整数ｍ回畳み込み累積である離散Ｂ－スプラインとの内積を演算し、スペクトルを離散Ｂ－スプラインの線形結合近似で表現する際の結合係数を求める第１の演算ユニット（１１）と、結合係数に基づきスペクトルを複数のＢ－スプライン分布に分解するピーク分離ユニット（１４）とを含む解析装置（１０）を提供する。

Description

解析装置
　本発明は、スペクトルを解析する装置および方法に関するものである。
　国際公開ＷＯ２０１２／０５６７０９号公報に開示されたサンプルを分析するユニットは、少なくとも２つのパラメータにより制御される電界を通過する化学物質のイオン強度を測定するイオン移動度センサーにサンプルを供給して得られた測定データに含まれるデータの２次元表現であって、第１のパラメータを変化させ、他のパラメータを固定したときのイオン強度を示す２次元表現の中に存在するピークを検出する機能ユニットと、検出されたピークと他の２次元表現の中に存在するピークとの連続性および生滅に基づいて検出されたピークを分類する機能ユニットと、分類されたピークに基づいて前記サンプルに含まれる化学物質を推定する機能ユニットとを有する。
　複数のピークが存在する可能性がある測定データの解析方法として、ガウシアン特性を用いたカーブフィッティング方法が知られている。さらに高速で複数のピークの存在を解析できる方法が求められている。
　本発明の一態様は、センサーから得られるスペクトルのサンプリング値と、スペクトルの第１の座標軸に沿って定義され、有限個のサンプリング点で正の値を持ち、他のサンプリング点では０（ゼロ、零）である局所非負離散関数との内積を演算し、スペクトルを複数の局所非負離散関数の線形結合近似で表現する際の結合係数を求める第１の演算ユニットと、結合係数に基づきスペクトルを複数の局所非負離散関数に対応する分布に分解する分離ユニットとを有する解析装置である。第１の演算ユニットは、得られた結合係数が負の場合は結合係数を０（ゼロ）に固定する。
　物理的な現象に関わる量（測定値、物理量）を捉えるセンサーから得られるスペクトルは、非負のいくつかの分布の結合であることが多い。特に、化学物質の存在に関する情報を出力するセンサー、たとえば、イオン移動度センサー、質量分析器、ガスクロマトグラフィ、ラマンなどの分光分析装置から得られるスペクトルは、いくつかの化学物質の存在に関連した非負の分布の結合であると解釈できる。したがって、適当な局所非負離散関数を用意して内積を計算して結合係数を求めることにより、無限遠までのカーブフィッティングを行わずに、高速でスペクトルを近似でき、スペクトルに含まれる分布を解析できる。
　さらに、物理的な現象を要因とするスペクトルを近似する際に、複数の局所非負離散関数の結合係数がマイナスになることはないと仮定してよい。このため、結合係数が負の場合は結合係数を０（ゼロ）に固定することにより、測定対象の物理的な現象に、よりマッチした分布の結合でスペクトルを近似できる。
　局所非負離散関数の例は、ハン窓関数などのレイズド・コサイン（Raised Cosine）、回転楕円体波動関数（Prolate Spheroidal Wave Function）、離散Ｂ－スプラインなどである。局所非負離散関数は、連続する整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０となる矩形離散関数の整数ｍ回畳み込み累積である離散Ｂ－スプラインを含むことが望ましい。曲線表現というＢ－スプラインの元来の用途ではないが、整数ｍを増やせばＢ－スプラインはガウス分布に近づき、Ｂ－スプラインを分布関数として利用できる。
　この場合、第１の演算ユニットは、スペクトルと離散Ｂ－スプラインとの内積を演算する内積演算ユニットを含み、分離ユニットは、結合係数に基づきスペクトルを複数のＢ－スプライン分布に分解するピーク分離ユニットとを含むことが望ましい。
　第１の演算ユニットは、以下の式（Ａ）の伝達関数を含む第１のディジタルフィルタユニットと、以下の式（Ｂ）の伝達関数を含む第２のディジタルフィルタユニットとを含むことが望ましい。離散Ｂ－スプラインはこれらの伝達関数を用いて内積を高速に演算できる。
（１／（１－ｚ-1））m・・・（Ａ）
（１－ｚ-n）m　　　　・・・（Ｂ）
　さらに、解析装置は、第１の演算ユニットにより、ｎおよびｍの少なくともいずれかが異なる結合係数候補を求め、結合係数候補による離散Ｂ－スプラインの線形結合近似とスペクトルとの誤差を評価して、結合係数を求める第２の演算ユニットをさらに有することが望ましい。スペクトルを近似するのに適した整数ｎおよび／またはｍを選択でき、ピーク分離性能を向上できる。
　第１の演算ユニットは、スペクトルに対して大きな正の内積を持つ離散Ｂ－スプラインを選んで逐次、線形結合近似で表現する際の結合係数を求めるユニットを含んでいてもよい。
　整数ｎおよび整数ｍの組み合わせの一例は以下の式（Ｃ）を満足するものである。なお、ＳＰｍａｘは、内積演算ユニットにおいて計算されるサンプリング数である。サンプリング点の最大値を超えた畳み込み累積をもつ離散Ｂ－スプラインを用いる意味はほとんどなく、式（Ｃ）を満足する整数ｍおよびｎを用いることにより演算時間を短縮できる。
ｎ×ｍ≦ＳＰｍａｘ　　　　・・・（Ｃ）
　解析対象のスペクトルの一例はイオン移動度センサーの出力であり、第１の座標軸の一例は補償電圧である。
　本発明の他の態様の１つは、センサーから得られるスペクトルの解析方法である。解析方法は以下のステップを含む。
１．スペクトルのサンプリング値と局所非負離散関数との内積を演算し、スペクトルを複数の局所非負離散関数の線形結合近似で表現する際の結合係数を求めること。
２．結合係数に基づきスペクトルを複数の局所非負離散関数に対応する分布に分解すること。
結合係数を求めることは、得られた結合係数が負の場合は結合係数を０に固定することを含む。
　局所非負離散関数は、連続する整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０（ゼロ）となる矩形離散関数の整数ｍ回畳み込み累積である離散Ｂ－スプラインを含むことが望ましい。この場合、結合係数を求めることは、スペクトルと離散Ｂ－スプラインとの内積を演算することを含み、分離することは、結合係数に基づきスペクトルを複数のＢ－スプライン分布に分解することを含むことが望ましい。
　また、内積を演算することは、上記式（Ａ）および（Ｂ）の伝達関数を演算することを含むことが望ましい。また、結合係数を求めることは、整数ｎおよび整数ｍの少なくともいずれかが異なる結合係数候補による離散Ｂ－スプラインの線形結合近似とスペクトルとの誤差を評価して、結合係数を求めることを含むことが望ましい。
　本発明のさらに異なる他の態様の１つは、センサーから得られるスペクトルを解析するプログラムである。プログラム（プログラム製品）は、プロセッサが、スペクトルのサンプリング値と局所非負離散関数との内積を演算し、スペクトルを複数の局所非負離散関数の線形結合近似で表現する際の結合係数を求めることと、結合係数に基づきスペクトルを複数の局所非負離散関数に対応する分布に分解することとを有し、結合係数を求めることは、得られた結合係数が負の場合は結合係数を０に固定することを含む。
解析装置の概略構成を示すブロック図。 Ｂ－スプラインを説明する図。 離散Ｂ－スプラインを説明する図。 イオン移動度センサー内のイオン（分子）の移動と、スペクトルとの関係を示す図。 内積を、ディジタルフィルタを用いて計算する様子を示す図。 結合係数から旋回結合近似をディジタルフィルタを用いて計算する様子を示す図。 Ｂ－スプライン近似のいくつかの例を示す図。 Ｂ－スプライン近似の異なる例を示す図。 解析方法の概略を示すフローチャート。 異なる解析方法の概略を示すフローチャート。
発明の実施の形態
　図１にイオン移動度センサー１の出力（測定データ、ＩＭＳデータ）を解析する解析装置の一例を示している。この解析装置１０は、イオン移動度センサー１から直に入力されたＩＭＳデータ１９を解析するものであってもよく、適当なストレージ（記録媒体）２に格納されたＩＭＳデータ１９を解析するものであってもよい。
　解析装置１０は、ＩＭＳデータ１９を離散Ｂ－スプラインにより近似するための結合係数を求める演算ユニット（第１の演算ユニット）１１と、求められた結合係数による近似の状況を評価する評価ユニット（第２の演算ユニット）１２と、これらのユニット１１および１２により使用される設定値などが格納されたライブラリ１３と、評価ユニット１２によりＩＭＳデータ１９を近似する離散Ｂ－スプラインの線形結合として選択されたＢ－スプライン近似に含まれるピークを分離および認識するピーク分離ユニット１４と、分離されたピーク（分布）の位置・高さ・形状・動向などの属性に基づいてピークを分類し、ＩＭＳデータ１９の測定対象を推定する推定ユニット１５と、分離されたピークに関するデータが含まれるデータベース１６とを含む。
　解析装置１０は、典型的には、プロセッサ（ＣＰＵ）とメモリとを含むコンピュータ資源を用いて構成され、ソフトウェア（プログラムまたはプログラム製品）として、適当な記録媒体に記録され、またはインターネットなどのコンピュータネットウェアを介して提供される。解析装置１０は、ＡＳＩＣ、ＬＳＩあるいは再構成可能な半導体回路装置（プロセッサ）、光回路装置などにより提供されてもよい。
　イオン移動度センサー１の典型的なものとしては、非対称電界イオン移動度スペクトロメータ（ＦＡＩＭＳ）または微分型電気移動度スペクトロメータ（ＤＭＳ）が知られている。ＦＡＩＭＳセンサー１の一例は、国際公開ＷＯ２００６／０１３３９６または国際公開ＷＯ２００７／０１４３０３に記載されたイオン移動度分析装置である。これらはモノリシックな、またはマイクロマシン型のコンパクトな質量分析装置であり、十分に携帯できるものである。
　ＦＡＩＭＳセンサー１は、イオン化された測定対象に電場の影響を与えながら移動させるドリフトチャンバと、ドリフトチャンバを通過したイオン化された測定対象（測定対象の電荷）を検出する検出器とを含む。ドリフトチャンバにおいては、電圧Ｖｄおよび電圧Ｖｃにより生成されるソフトウェア制御された電界が特定の周期でプラス・マイナスに変動し、その電界のフィルタリング効果により、検出目標の化学物質がフィルタリングされ、短期間、たとえば、ｍｓｅｃレベルで検出器に衝突し、イオン強度（イオン電流）Ｉｃとして測定される。
　Ｖｄ電圧（Ｄｉｓｐｅｒｓｉｏｎ Ｖｏｌｔａｇｅ、分散電圧、電界電圧（Ｖｄ））は交流成分であり、Ｖｃ電圧（Ｃｏｍｐｅｎｓａｔｉｏｎ Ｖｏｌｔａｇｅ、補償電圧）は直流成分である。ＦＡＩＭＳセンサー１では、これら２つの変量に応じて変化するイオン強度をイオン電流として検出することが可能である。検出されたイオン電流の波形は化学物質により異なることが多いので、化学物質の特定を行うことが可能である。イオン移動度センサーから得られるデータ（測定データ、ＩＭＳデータ）１９は、イオン電流の値を含めると少なくとも３次元のデータとなる。また、化学物質によっては、たとえばイオン移動度特性により、Ｖｄ電圧を固定した場合に、Ｖｃ電圧を変えるとイオン電流の値は複数のピークを形成することがある。したがって、ＩＭＳデータ１９に含まれるピークに基づき、イオン移動度センサー１で検出しているサンプルの化学組成（サンプルに含まれる化学物質）を特定できる可能性がある。
　図２に、ｍ階のＢ－スプラインを示している。図２（ａ）は、矩形関数を１回畳み込み積分したＢ－スプラインであり、１階のＢ－スプラインである。図２（ｂ）は、２階のＢ－スプラインを示し、図２（ｃ）は、３階のＢ－スプラインを示し、図２（ｄ）は、４階のＢ－スプラインを示す。ｍを無限大にした極限では正規分布に収束することは中心極限定理で明らかである。
　図３に、矩形関数から得たｎ個の標本数列（サンプリング点）のｍ階畳み込み累積により離散Ｂ－スプラインが定義される様子を示している。ｍ階の離散Ｂ－スプラインは、以下の式（０）に示すようにｚ変換することにより、伝達関数で定義されるディジタルフィルタの出力として得られる。
Figure JPOXMLDOC01-appb-M000001
　イオン移動度分析は、化学物質分子の判別分析の手法である。微量の化学物質を分析できるので、匂い・香りの分析、毒物・薬物・爆薬の探知に利用できる。分析手順は、物理的処理と計算的処理に大別される。
　図４（ａ）および（ｂ）にイオン移動度センサー１における物理的処理を模式的に示している。イオン化された化学分子がイオン移動度センサー１の内部で移動し、正極１ａまたは負極１ｂにキャッチされる。典型的なイオン化法は、放射性元素からの電子線照射であり、ＵＶ照射またはコロナ放電などの方法であってもよい。イオン移動度センサー１の内部でイオン化された分子が移動中に、上述したように、分散電圧（Ｖｄ電圧）と、補償電圧（Ｖｃ電圧）とによりジグザグに動き、Ｖｄ電圧とＶｃ電圧との条件により電極１ａまたは１ｂに到達した分子のみがイオン移動度センサー１の測定データ（ＩＭＳデータ）１９として測定される。ＩＭＳデータ１９の一例は、Ｖｃ電圧の変化によりイオン電流Ｉｃが変動するスペクトル（スペクトル分布、スペクトラム、波形分布、出力分布）である。
　イオン化された分子は、その電荷／質量に比例してＶｄ電圧とＶｃ電圧との影響を受け、途中で空気分子に衝突して拡散しながら進む。このため、平均的にはイオンの種類３ａ～３ｃに特有な速度で進み、時間的にばらついて電極、たとえば正極１ａに到達する。ばらつきは正規分布に従うと仮定されている。各衝突で生じる位置変化が同一分布に従ってばらつくならば、無限回の衝突を経て得られる分布がガウス分布となるからである。
　図４（ｃ）にＩＭＳデータ１９の一例を示している。この波形はＶｄ電圧を一定としてＶｃ電圧を変えたときに得られるイオン電流のスペクトル（スペクトル分布）であり、この電流波形をプロファイルと呼び、第１の座標軸はＶｃ電圧（補償電圧）になる。１つの種類のイオンに１つのガウス分布が対応するので、プロファイルは図４（ｄ）のように複数１９ａ～１９ｃのガウス分布の非負荷重和としてモデル化される。
　従来型の計算処理は、与えられたプロファイルからそれを構成する各々の分布を同定することを目標としている。イオン分子の種類が各分布の平均値によって判明し、量は高さで判明する。ガウス分布の個数・平均・分散が未知であるため、これらのパラメータを最急降下法で探索することが一般的な手法である。探索は逐次処理であり、また、収束までに多くの繰り返し演算を要する。このため、現状の計算処理は高性能ＣＰＵを用いても多くの時間を要している。
　物理的処理の小型・高速化の一例は、２００６年にオウルストーン社からリリースされたＦＡＩＭＳチップによって成されたものである。センサチップの大きさは数十ミリメートル角、物理的処理の速度は数ミリ秒である。したがって、計算処理を効率化することによりＦＡＩＭＳなどの化学物質分析センサー１および解析装置１０を含めたシステム全体を携帯機器に組込むことも可能になる。
　そのためには、計算処理アルゴリズムを、伝統的なガウス混合分布モデルに基づく逐次探索からチップ化に適するような新たな計算処理に変更することが望ましい。そのアプローチの特徴は、以下の４点にまとめられる。
　（１）負の無限大から正の無限大までの広がりを有するガウス分布の代わりに、有限個のサンプリング点で正の値を持ち、他のサンプリング点では０（ゼロ）である局所非負離散関数を用いる。処理時間が短縮されるとともに、物理現象が負の無限大から正の無限大までの広がりを有することは稀であり、局所非負離散関数の方が物理現象に基づく分布を近似するために適している場合も多い。典型的な局所非負離散関数の１つとして、一様分布に従う有限回(ｍ回)の衝突を経て得られる速度のばらつきを表す整数ｍ階のＢ－スプラインを用いることができる。曲線表現というＢ－スプラインの元来の用途ではないが、整数ｍを増やせばＢ－スプラインはガウス分布に近づく。たとえば、イオン移動度分析の原理によれば、分布の形も、また、無限の裾野をもつ完璧な正規分布であるはずはない。正規分布の代わりに、Ｂ－スプラインを用いることも不自然ではない。
　（２）Ｂ－スプラインの高速計算法を利用する。一様連続分布をｎ点で標本化した一様離散分布で置き換えて得られる離散Ｂ－スプラインが加算と減算だけで生成できる（T. Saramaki, Y. Neuvo and S. K. Mitra, Design of computationally efficient interpolated FIR filters, IEEE Trans. Circuits & Syst., 35(1), 70-88, 1988  参照）。標本点数ｎを増やすことによって離散Ｂ－スプラインが元々のＢ－スプラインに収束し、高速ディジタルフィルタにより波形を離散Ｂ－スプライン近似できる（K. Ichige and M. Kamada, An approximation for discrete B-splines in time domain, IEEE Signal Process. Lett., 4, 82-84, 1997、および、K. Ichige, M. Kamada and R. Ishii, A simple scheme of decomposing and reconstructing continuous-time signals by B-splines, IEICE Trans. Fundamentals, E81-A, 2391-2399, 1998を参照）。
　（３）近似に用いる離散Ｂ－スプラインを極めて細かい間隔で配置する。曲線表現というＢ－スプラインの元来の用途では、表現に必要十分な間隔でＢ－スプラインを配置することが常識であった。イオン移動度分析では分布の中心の位置が未知であるが、位置を探索すれば、遅い逐次計算に後戻りしてしまう。そこで、予め離散Ｂ－スプラインを細かい間隔で配置しておき、プロファイルの近似にとって主要な離散Ｂ－スプラインの荷重だけが大きくなって浮き彫りになることを狙う。
　（４）未知パラメータである分布の分散を逐次的に探索するのではなく、同時期に到達するイオンの広がりは同程度であると仮定して、様々な値に固定したｎについて高速ディジタルフィルタによる近似を総当り的に並列実行する。少数の離散Ｂ－スプラインによる高精度な近似が得られる場合を分析結果とする。近似精度の判定に利用できる方法の１つとしてはスパース近似を挙げることができる（スパース近似については、M. Vetterli, P. Marziliano and T. Blu, Sampling signals with finite rate of innovation, IEEE Trans. Signal Process., SP-50, 1417-1428, 2002を参照）。
　ｍ階のＢ－スプラインは、図２（ａ）に示される矩形関数の整数ｍ回畳込み積分として定義されている。これが、整数ｍを無限大にする極限で正規分布に収束することは、中心極限定理が示すとおりである。さらに、矩形関数から得たｎ個の標本値列のｍ回畳込み累積である離散Ｂ－スプラインが、整数ｎを無限大にする極限でＢ－スプラインに収束することも知られている。解析装置１０の演算ユニット１１においては、この離散Ｂ－スプラインを分布として用いる。離散Ｂ－スプラインは、整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０となる矩形離散関数の整数ｍ回畳込み積分として定義される。
　なお、整数ｍおよび整数ｎは１以上の整数であり、解析に供される標本数（サンプリング数）ＳＰｍａｘに対し以下の関係となる。
ｎ×ｍ≦ＳＰｍａｘ　　　　・・・（Ｃ）
　スペクトルの長さ（サンプリングポイント数）ＳＰｍａｘが５１２である場合、その長さよりも幅の広いＢ－スプラインを用いてもよいが、内積計算には寄与しないので処理速度的には無駄となる。
　１階の離散Ｂ－スプラインは以下の（１）式で定義される。また、そのｚ変換は、以下の（２）式で表現される。
Figure JPOXMLDOC01-appb-M000002
　したがって、再帰的に定義されるｍ階の離散Ｂ－スプラインおよびそのｚ変換は、以下に示す式（３）および式（４）によりそれぞれ表現される。
Figure JPOXMLDOC01-appb-M000003
　離散Ｂ－スプライン（ｂｍ[k-l]）と（ｂｍ[k-r]）の内積は、以下の式（５）のように変形される。以下において、「＜関数ａ,関数ｂ＞」は、関数ａと関数ｂの内積をとって数値１つを求める演算を示す。たとえば、＜p[・-l]，ｂｍ[・-r]＞は、式（５）のように変数ｋについてマイナス無限大からプラス無限大まで数列の和（総和）を基本的には示し、内積の値が変数ｋに依存しないため、変数ｋは「・」にて代用する。また、実際の内積計算においてサンプリング数（標本数）は無限ではなく有限、すなわち、最大値ＳＰｍａｘであり、変数ｋの範囲は（０，１，２，・・・ｍ×ｎ－１）である。しかしながら、サンプリング値以外の内積は０になるので、本明細書においては無限和の表現を用いることがある。
Figure JPOXMLDOC01-appb-M000004
　ここで、「~bｍ[k]＝ｂｍ[-k]」とおけば、これは、さらに以下の式（６）のように変形される。
Figure JPOXMLDOC01-appb-M000005
　また、そのｚ変換は、以下の式（７）のように表現される。
Figure JPOXMLDOC01-appb-M000006
　このため、この内積は以下の式（８）として計算できる。
Figure JPOXMLDOC01-appb-M000007
　ｚ変換が以下の式（９）である係数（結合係数）ｃ[l]によるｍ階Ｂ－スプラインｂｍ[k-l]の線形結合ｑ[k]は、以下の式（１０）なる畳込みで表される。
Figure JPOXMLDOC01-appb-M000008
　この演算は、ｚ変換したときの以下の式（１１）に対応する。
Figure JPOXMLDOC01-appb-M000009
　したがって、「（（１－ｚ-n）／（１－ｚ-1））m」なる伝達関数をもつディジタルフィルタに結合係数ｃ[k]を入力して得られる出力を時刻ｋにおいてサンプリングすれば線形結合ｑ[k]が求められる。このディジタルフィルタは、式（Ａ）のｍ重累積と、式（Ｂ）のｍ階差分に分けて実行可能である。
（１／（１－ｚ-1））m・・・（Ａ）
（１－ｚ-n）m　　　　・・・（Ｂ）
　解析装置１０の演算ユニット１１の内積演算ユニット１１ｄは、式（Ａ）の伝達関数を含む第１のディジタルフィルタユニット１１ａと、式（Ｂ）の伝達関数を含む第２のディジタルフィルタユニット１１ｂとを含む。これらのディジタルフィルタユニット１１ａおよび１１ｂにおける計算途中にオーバーフローは生じるが、これらのフィルタユニット１１ａおよび１１ｂは、計算の最終結果を収容できるビット数を擁する、「２」の補数演算を行う機能を備えており、最終結果は正しい値となる。
　なお、最終結果である線形結合ｑ[k]の絶対最大値は、離散Ｂ－スプラインが非負であることから導かれる以下の不等式（１２）により、結合係数の絶対値｜ｃ[l]｜の最大値のｎm倍を越えない。
Figure JPOXMLDOC01-appb-M000010
　解析装置１０に供給されるＩＭＳデータ１９のスペクトルの離散的なサンプリング値、すなわちプロファイルをｐ[k]で表すと、そのｚ変換は以下の式（１３）に示したように表現される。
Figure JPOXMLDOC01-appb-M000011
　プロファイルｐ[k]と、離散Ｂ－スプラインｂｍ[k-r]の内積は以下の式（１４）で示され、式（１５）のように変形できる。
Figure JPOXMLDOC01-appb-M000012
　ここで、（~ｂｍ[k]）のｚ変換はＢｍ(z-1)となる。したがって、ｚ変換した表現では、以下に示す式（１６）に対応する。この内積は、「（（１－ｚ-n）／（１－ｚ-1））m」なる伝達関数をもつディジタルフィルタにプロファイルｐ[r]を入力して得られる出力を時刻（(ｎ-1)ｍ＋ｒ）においてサンプリングすれば求められる。この処理は、上記と同様に、２つの部分から成る第１のディジタルフィルタユニット１１ａおよび第２のディジタルフィルタユニット１１ｂで実行可能である。
Figure JPOXMLDOC01-appb-M000013
　プロファイルｐ[k]と離散Ｂ－スプラインｂｍ[k-r]との内積は、通常の考え方では、上記の式（１４）の通りに積和演算で計算する。しかしながら、式（１５）と（１６）によれば、加算と減算だけで構成されるディジタルフィルタユニット１１ａおよび１１ｂで計算できることが分かる。したがって、内積演算ユニット１１ｄにおいては演算用の構成を簡略化でき、内積演算を高速で処理できる。
　図５に、内積演算ユニット１１ｄのディジタルフィルタユニット１１ａおよび１１ｂにプロファイルｐ[r]を入力し、（(ｎ-1)ｍ＋ｒ）における内積＜p[・]，ｂｍ[・-r]＞を計算する様子を示している。前半部分である第１のディジタルフィルタユニット１１ａで計算するｍ階累積は整数ｎに依存せず、後半部分である第２のディジタルフィルタユニット１１ｂで計算するｍ階差分にだけ整数ｎが含まれている。そのため、ｍ階累積の計算は各プロファイルに対して１度だけでよく、第２のディジタルフィルタユニット１１ｂで整数ｎを変えながらｍ階差分を計算するだけで、異なる整数ｎについての離散Ｂ－スプラインとプロファイルとの内積を取り出すことが可能である。１つの内積の計算に要する演算は、ほとんどｍ回の減算だけであり、処理時間を短縮できる。
　第１の演算ユニット１１は、さらに、非負係数最小自乗近似を行う最小自乗近似ユニット（最小二乗近似ユニット）１１ｃを含む。この最小自乗近似ユニット１１ｃでは、プロファイルｐ[k]を、間隔「１（one）」で配置されたＫ個の離散Ｂ－スプラインを結合係数ｃ[l]で線形結合した関数ｑ[k]で近似したときの自乗誤差を計算する。線形結合近似ｑ[k]は、以下の式（１７）で表され、自乗誤差は以下の式（１８）で計算される。
Figure JPOXMLDOC01-appb-M000014
　自乗誤差（二乗誤差）Ｅが最小になるような係数（結合係数）ｃ[l]は、通常の最小自乗近似法に従って以下の線形方程式（１９）を解くことによって定まる。
Figure JPOXMLDOC01-appb-M000015
　解析の対象であるスペクトル１９は、イオン移動度センサー１の出力である。したがって、係数はイオンの個数に比例すると想定できる。このため、係数ｃ[l]には非負であるという制約条件が課せられる。イオン移動度センサー１の出力以外であっても、化学物質の存在を示すスペクトル（分布）などの物理量であれば、係数ｃ[l]は非負である。最小自乗近似ユニット１１ｃは、この条件下で係数を定めるために通常の最小自乗近似で得られた係数のうち負となったものを０（ゼロ）に固定し、それに対応する離散Ｂ－スプラインを除外して最小自乗近似を、全ての係数が非負になるまで繰り返す機能を含む。線形方程式（１９）を解くことは、浮動小数点演算の四則演算を要するが、解くべき線形方程式はバンド行列で表現できる。また、対象となる関数がだんだん減っていくので、それほど処理時間が増加する方向にはならない。離散Ｂ－スプライン同士の内積＜ｂｍ[・-k]、ｂｍ[・-l]＞は、予め計算してライブラリ１３に保存しておくことも可能である。
　細かい間隔で配置した離散Ｂ－スプラインの列は、理論的には線形独立であるものの、実際にはほとんど線形従属である。このため、荷重を決定する処理が数値的に不安定となる可能性が指摘され得る。しかしながら、テストデータに対する数値実験では、多数の離散Ｂ－スプラインを用いる初期の段階では数値的に不安定になるものの、負の係数を得た離散Ｂ－スプラインが除外されていく過程で近接する離散Ｂ－スプラインの個数が減っていく。このために最終的には計算が安定することが観測されている。すなわち、最小自乗近似ユニット１１ｃにおいて結合係数を非負に制限していることが数値的に不安定になることを抑制していると判断される。
　図６に、第１のディジタルフィルタユニット１１ａおよび第２のフィルタユニット１１ｂを用いて、上記で求められた結合係数ｃ[k]からプロファイルを近似するために求めた線形結合近似ｑ[k]を計算する様子を示している。結合係数を非負に制限しているため、近似用の線形結合で用いられる結合係数の数は限られ、プロファイルｐ[k]に含まれるピークを実質的に分離できる。
　解析ユニット１０は、さらに、第１の演算ユニット１０により求められた結合係数により得られる線形結合近似ｑ[k]を評価する評価ユニット（第２の演算ユニット）１２を含む。評価ユニット１２は、さまざまな整数ｎおよび／または整数ｍに対して結合係数を結合係数候補として事前に求め、それらにより得られる線形結合近似ｑ[k]をプロファイルｐ[k]に対して評価し、結合係数候補の中から最適な結合係数と、それにより線形結合近似を出力する。すなわち、この評価ユニット１２は、想定される様々なｍとｎとについて、プロファイルに対する離散Ｂ－スプラインの線形結合近似を総当たりで評価する。離散Ｂ－スプラインの内積演算は上述したように短時間で処理できるので、総当たりで評価することが現実的に可能となる。
　評価ユニット１２は自乗誤差の評価を行う第１の評価ユニット１２ａと、スパーシティ（Sparsity）評価を行う第２の評価ユニット１２ｂとを含む。第１の評価ユニット１２ａは自乗誤差を求める。結合係数ｃ[k]から線形結合近似ｑ[k]を計算することは、式（１０）に従って図６に示したように第１のディジタルフィルタユニット１１ａおよび第２のディジタルフィルタユニット１１ｂによって短時間に実行できる。これにより得られた線形結合近似ｑ[k]と、与えれたプロファイルｐ[k]の自乗誤差は、(ｐ[k]-ｑ[k])2を累積して計算することにより求められる。第１の評価ユニット１２ａは、平均自乗誤差Ｅ１を以下の式（２０）により計算する。
Figure JPOXMLDOC01-appb-M000016
　第２の評価ユニット１２ｂは、様々な整数ｎについて得られた結合係数候補の中から第１の評価ユニット１２ａにより自乗誤差Ｅ１が十分に小さくなった場合を分析結果の候補とする。第２の評価ユニット１２ｂは、結合係数候補の中から、さらに、スパースな良い近似結果が得られる結合係数の組み合わせを探る。これは、得られた候補の中で、少数の離散Ｂ－スプラインでプロファイルｐを近似できていることが正しい分析の本質的な条件であるとの前提によるものである。
　プロファイルｐ[k]も線形結合近似ｑ[k]も非負である。このため、自乗誤差が極めて小さければ、それらの面積Σｋ（ｐ[k]）とΣｋ（ｑ[k]）とは殆ど同じになるはずである。そこで、離散Ｂ－スプラインｂｍ[k-l]を、その面積「Σｋ（ｂｍ[k-l]）＝ｎm」で割って正規化した近似表現を以下の式（２１）のように考える。
Figure JPOXMLDOC01-appb-M000017
　修正された係数ｎmｃ[l]の和は、以下の式（２２）のようにほぼ一定になる。
Figure JPOXMLDOC01-appb-M000018
　この状況において、スパーシティ（sparsity）は、線形結合近似ｑ[k]が大きな少数の塊で構成されていることを意味する。その評価量の１つは、以下の式（２３）に示す評価量Ｅ２である。第２の評価ユニット１２ｂは、評価量Ｅ２が最大となる場合を分析結果として選び、出力する。
Figure JPOXMLDOC01-appb-M000019
　図７（ａ）～（ｊ）に、整数ｍを４に固定し、様々な整数ｎについて、長さ１２８点のサンプル値という局所的な窓の範囲で、プロファイルｐ[k]の１つの例を近似した結果を示している。プロファイルは、イオン移動度センサー１の出力（イオン電流、ＩＭＳデータ、スペクトル）１９であり、座標軸（第１の座標軸）は補償電圧Ｖｃである。図中で、実線はプロファイルを示し、破線は近似曲線（線形結合近似）ｑ[k]を示し、一点鎖線は近似曲線ｑ[k]を構成する離散Ｂ－スプラインを表す。
　図７（ａ）～（ｊ）において、第１の評価ユニット１２ａにより得られる近似誤差Ｅ１は整数ｎが６～１２の場合に１０％以下であり、その他の場合には１５％より大きい。したがって、第２の評価ユニット１２ｂは、前者の場合、すなわち、整数ｎが６～１２の場合の近似曲線（線形結合近似）ｑ[k]について評価量Ｅ２を計算する。整数ｎが６～１２のうちで、スパーシティの評価量Ｅ２が最大になるのは、整数ｎが１２の場合である。この場合をスパースな良い近似として採用し、第２の評価ユニット１２ｂは出力する。この例では、プロファイルｐには２つのピークが含まれているとの結果が出力される。
　図８は、１つのプロファイル１９全体についてサンプリングし、それらの値に対する近似結果の例である。最適な整数ｎは、第１の座標軸の経過にしたがって変わってもよい。たとえば、最適な整数ｎは局所的な補償電圧Ｖｃ値の窓ごとに推定することも可能である。図８に示した例では、左側のピークと右側のピークに対する整数ｎの最適値は、それぞれ１２と１４となった。
　解析装置１０のピーク分離ユニット１４は、上記の結果に基づいて分離されたピーク（検出されたピーク）によりピークデータベース１６を更新する。ピーク分類および推定ユニット１５は、ピークデータベース１６のピーク分類ルールテーブルおよびピーク種別（タイプ）ルールテーブルにより提供された評価ルールを用い、ピークの種別を判断する。
　評価ルールには、リアクタント（ＲＩＰ）のピーク特性と、化学物質から派生したモノマー、ディマーおよびトリマーのピーク特性の理論上の関係が含まれる。さらに、分類および推定ユニット１５は、イオン移動度センサー１が測定しているサンプルガスに含まれる化学物質を推定する機能を含む。分離および推定ユニット１５は、モノマー、ディマーおよびトリマーの少なくともいずれかのピークを示す特定されたピークに基づきサンプルガスに含まれる化学物質の候補を選択する選択機能と、特定されたピークと第３のパラメータである濃度との相関を求め、サンプルガスに含まれる化学物質の候補を選択する選択機能などを含む。
　図９に、解析ユニット１０における解析方法の概略をフローチャートにより示している。解析方法（解析処理）は、解析ユニット１０を構成するコンピュータまたはプロセッサを制御するプログラム（プログラム製品）として適当な記録媒体に記録され、またはコンピュータネットワークを介して提供され、解析ユニット１０に実装される。解析処理は、ローカルで実行されてもよく、コンピュータネットワークを介してクラウド上、あるいはサーバー上で実行されてもよく、ローカルとクラウドとで分散して処理されてもよい。
　演算ユニット１１において解析処理を開始する際に、ステップ５１において、局所非負離散関数を設定する。以下においては、離散Ｂ－スプライン（式（３））を用いて説明するが、局所非負離散関数は、レイズド・コサイン（Raised Cosine、ハン窓関数）、回転楕円体波動関数（Prolate Spheroidal Wave Function）などの、有限個のサンプリング点で正の値を持ち、他のサンプリング点では０（ゼロ、零）である関数であればよい。
　ステップ５２において、センサー１の出力であるスペクトルデータ（ＩＭＳデータ）１９を取得し、プロファイルｐ[・]を設定する。ステップ５３において、プロファイルｐ[・]と離散Ｂ－スプラインｂｍ[・-r]との内積を計算するため式（Ａ）および（Ｂ）に対応するディジタルフィルタを準備し、高速で内積を計算する環境を整える。離散Ｂ－スプラインｂｍ[・]は、整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０となる矩形離散関数の整数ｍ回畳み込み累積である。
　ステップ５４においてディジタルフィルタを用い、プロファイル（スペクトル）ｐ[k]を複数の離散Ｂ－スプラインｂｍ[k-r]の線形結合近似ｑ[k]（式（１７））で表現する際の結合係数であって、非負の結合係数を最小自乗近似により求める。具体的には、最小自乗近似ユニット１１ｃで説明したように、プロファイルｐ[k]を、間隔「１」で配置されたＫ個の離散Ｂ－スプラインを結合係数ｃ[l]で線形結合した関数ｑ[k]で近似したときの自乗誤差を計算する。プロファイルｐ[k]はセンサー１の出力であるので、非負の局所離散関数の線形結合の係数ｃ[l]も非負とする。したがって、最初に最小自乗近似で得られた係数のうち負となったものを０（ゼロ）に固定し、それに対応する離散Ｂ－スプラインを除外して最小自乗近似を、全ての係数が非負になるまで繰り返す。
　ステップ５５において、予め設定されたさまざまな整数ｎおよびｍについてステップ５３および５４の処理を行う。整数ｎおよびｍは、上記式（Ｃ）の範囲内で選択される。整数ｎおよびｍは、センサーの測定対象、センサーの種類および特性などにより、線形結合で近似するのに適したものに限定してもよく、さらに処理速度を向上できる。
　ステップ５６において、評価ユニット１２が、上記により得られた線形結合近似ｑ[k]の中から、プロファイルｐ[k]との誤差を評価して、結合係数を求める。上述したように、平均自乗誤差Ｅ１（式（２０））が小さいことと、スパース近似の評価量Ｅ２（式（２３））が大きいこととを使うことができる。また、分析の良さを評価する単一の評価量Ｅ０として以下の式（２４）を用いてもよい。評価量Ｅ０は、平均自乗誤差Ｅ１とスパース近似の評価量Ｅ２とを組み合わせたものであり、評価量Ｅ０が小さいほど分析結果が良いことを示している。
Figure JPOXMLDOC01-appb-M000020
　ステップ５７において、ピーク分離ユニット１４が、誤差評価結果の良い線形結合近似ｑ[k]の各結合係数に基づきプロファイル（スペクトル）ｐ[k]を、それぞれの結合係数の離散Ｂ－スプラインに対応するＢ－スプライン分布に分離（分解）する。ステップ５８において、推定ユニット１５が、分離されたピークを分類し、ＩＭＳデータ１９の測定対象（化学物質、それらの組合せ、結合など）を推定する。
　図１０に、最小自乗近似ユニット１１ｃにおいて、非負係数を最小自乗法により求めるステップ５４の異なる例を示している。この例（解析方法）は、離散Ｂ－スプラインｂｍ[・-r]の全てを始めから近似に用いるのではなく、近似対象のプロファイルｐ[k]に対して大きな正の内積をもつものから選んで、逐次追加する手法を含む。これにより、処理時間を短縮でき、さらに、数値的に不安定となるリスクを避けることもできる。
　ステップ６１において、繰り返し処理のための初期設定を行う。次に、ステップ６２において、最大の＜ｐ[・]，ｂｍ[・-r]＞を与える整数ｒを選択する。たとえば、第１番目に選択される整数ｒをｒ１とする。ステップ６３において、そのｂｍ[・-r1]によるｐ１[・]の最小二乗近似の結合係数ｃ[r1]を求める。
　ステップ６４において、この係数ｃ[r1]が非負ならば、線形結合近似ｑ[・]を求める。一方、最初に求められる係数ｃ[r1]が負ならば、唯一の係数ｃ[r1]が０と判断されることになるので、ステップ６７において線形結合近似が求められない。このため計算を終了する（ステップ６８）。
　係数ｃ[r1]が正であれば、ステップ６９において平均自乗誤差Ｅ１を計算し、平均自乗誤差Ｅ１が閾値以下であれば良い線形結合近似ｑが得られたとして計算を終了する（ステップ７５）。所定の整数ｎおよびｍを用いて十分に精度の高い線形結合近似ｑが得られたと判断される場合は誤差評価（ステップ５６）またはピーク分離（ステップ５７）に進んでもよく、他の整数ｎおよびｍについて線形結合近似ｑを求める場合は、それらの値を変えて処理を繰り返す。
　ステップ６９において平均自乗誤差Ｅ２が閾値以上であれば、第１回目の線形結合近似ｑ[・]を線形結合近似ｑ１[・]とし、ステップ７０において、プロファイルｐ[・]に対する線形結合近似ｑ１[・]の残差「ｐ１[・]＝ｐ[・]－ｑ１[・]」を計算し、ステップ７１において整数ｉを（ｉ＋１）に増加する。つぎに、ステップ６２にもどって、次の最大の＜ｐ１[・]，ｂｍ[・-r]＞を与える整数ｒを選択し、それをｒ２とする。
　ステップ６３において、離散Ｂ－スプラインｂｍ[・-r1]」とｂｍ[・-r2]によるプロファイルｐ[・]の最小自乗近似の結合係数ｃ[r1]とｃ[r2]とを求める。ステップ６４において、これらの係数のうち負になるものがあれば、ステップ６５において、それらを０に固定し、ステップ６６において、それらの係数に対応する離散Ｂ－スプラインを除外する。ステップ６３に戻り、最小自乗近似を、残された係数、この場合は２つの係数が非負になるまで繰り返して、線形結合近似ｑ[・]を求める。２つの係数がどちらも０になれば、ステップ６８で計算を終了し、自乗誤差（二乗誤差）Ｅ１が十分に小さくなれば、ステップ７５で計算を終了する。
　上記以外の場合は、ステップ７０からの処理を繰り返す。すなわち、この段階で得られた線形結合近似ｑ[・]を２回目の線形結合近似ｑ２[・]とし、プロファイルｐ[・]に対する線形結合近似ｑ２[・]の残差ｐ２[・]（＝ｐ[・]－ｑ２[・]）を計算する。次にステップ６２に戻って次の最大の内積を与えるｒを選択する。
　一般には、残差ｐｉ[・]（＝ｐ[・]－ｑｉ[・]（ｉ＝１，２，３, …））に対して、最大の＜ｐｉ[・]，ｂｍ[・-r]＞を与えるrを選択し、それをｒ（ｉ＋１）とする（ステップ６２）。Ｂ－スプラインｂｍ[・-r1]～ｂｍ[・-r(i+1)]によるプロファイル残差ｐｉ[・]の最小二乗近似の結合係数ｃ[r1]～ｃ[r(i+1)]を求める（ステップ６３）。これらの係数のうち負になるものがあれば、それを０に固定し（ステップ６４および６５）、それに対応する離散Ｂ－スプラインを除外して（ステップ６６）、最小自乗近似を、ｉ＋１個の全ての係数（残った係数）が非負になるまで繰り返して、線形結合近似ｑ[・]を求める。すべての係数が０になるか、自乗誤差Ｅ１が十分に小さくなれば、計算を終了する（ステップ６７、６８、６９、７５）。
　計算を継続できる場合は、この線形結合近似ｑ[・]を（ｉ＋１）番目の線形結合近似ｑ（ｉ＋１）[・]とし、プロファイルｐ[・]に対する線形結合近似ｑ（ｉ＋１）[・]の残差ｐ（ｉ＋１）[・]（＝ｐ[・]－ｑ（ｉ＋１）[・]）を計算する（ステップ７０）。上記を繰り返して、計算が終了したときには、すべての係数が０の場合には近似が不能であると判定し、そうでない場合には自乗誤差Ｅ１が十分に小さくなった良い近似が得られたことになる。
　以上に説明したように、上記においては、離散Ｂ－スプラインの計算法に関する計算機構をイオン移動度分析に初めて適用した例で本発明を説明しているが、他の局所非負離散関数であってもよいことは説明した通りである。離散Ｂ－スプラインの配置間隔は整数ｎに固定されている必要はなく、上記の本応用場面では最も細密な間隔「１」が適切であった。また、ほとんど線形従属な関数系による最小自乗近似を行うことには、数値的な不安定性を生じる恐れがあるといわれている。しかしながら、上記のような非負の条件を導入することにより、いままでに試した数値例では、最終的な処理結果には大きな影響が出ていない。
　上記の例では、ガウス分布の代わりにその良い代替である離散Ｂ－スプラインを例として、イオン移動度分析で得られるプロファイル波形の分析に役立つハードウェア実現向きの手法を構成できることを説明した。離散Ｂ－スプラインも用いたシステムは、特に、プロセッサなどのハードウェア処理で実現しやすいシステムとして提供できる。ハードウェア実現に適するのは、与えられた波形と離散Ｂ－スプラインとの内積および離散Ｂ－スプラインの荷重和が極めて簡単なディジタルフィルタで計算できることが１つの要因である。
　このシステムにおいては、細かい間隔で配置した離散Ｂ－スプラインの荷重和でプロファイル波形を、荷重が非負であるという制約条件の下で最小自乗近似することによって、プロファイル中に埋もれている主要な成分の荷重が大きくなって検出される。細かい間隔で配置された離散Ｂ－スプラインは、ほとんど線形従属であるため、それらによる最小二乗近似は数値的には不安定になるが、数値計算の過程で負になった荷重を与えられた離散Ｂ－スプラインを除いていくため、最終的には安定的に近似結果が得られる。１回の近似が高速に行えるので、様々な幅の離散Ｂ－スプラインを試して、大きくて少ない荷重を用いた良いスパース近似が得られた場合を分析結果として採用できる。
　上記において、評価ユニット１２は、平均自乗誤差Ｅ１とスパーシティの評価量Ｅ２とを組み合わせて、分析の良さを評価する単一の評価量Ｅ０を設定してもよく、様々な電流波形に適用して数値的安定性・分析性能・処理速度を評価してもよい。また、解析する対象であるスペクトル（プロファイル）は、イオン移動度センサー１の出力に限定されず、本発明は、質量分析器、ラマン分光分析、赤外分光分析などの他のタイプのセンサーの出力にたいしても適用できる。解析装置１０の解析対象は、一様分布に従う有限回の衝突を伴うような現象の測定結果などの正規分布に従うと推定される現象の測定結果、分光、吸光などの現象の測定結果であってもよい。また、スペクトルの第１の座標軸は補償電圧に限定されず、他の物理量、時間などであってもよい。

Claims (13)

  1. 　センサーから得られるスペクトルのサンプリング値と、前記スペクトルの第１の座標軸に沿って定義され、有限個のサンプリング点で正の値を持ち、他のサンプリング点では０である局所非負離散関数との内積を演算し、前記スペクトルを複数の局所非負離散関数の線形結合近似で表現する際の結合係数を求める第１の演算ユニットであって、得られた結合係数が負の場合は結合係数を０に固定する第１の演算ユニットと、
    　前記結合係数に基づき前記スペクトルを前記複数の局所非負離散関数に対応する分布に分解する分離ユニットとを有する解析装置。
  2. 　請求項１において、局所非負離散関数は、連続する整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０となる矩形離散関数の整数ｍ回畳み込み累積である離散Ｂ－スプラインを含み、
    　前記第１の演算ユニットは、前記スペクトルと前記離散Ｂ－スプラインとの内積を演算する内積演算ユニットを含み、
    　前記分離ユニットは、前記結合係数に基づき前記スペクトルを複数のＢ－スプライン分布に分解するピーク分離ユニットを含む、解析装置。
  3. 　請求項２において、前記第１の演算ユニットは、以下の式（Ａ）の伝達関数を含む第１のディジタルフィルタユニットと、
    　以下の式（Ｂ）の伝達関数を含む第２のディジタルフィルタユニットとを含む、解析装置。
    （１／（１－ｚ-1））m・・・（Ａ）
    （１－ｚ-n）m　　　　・・・（Ｂ）
  4. 　請求項２または３において、前記第１の演算ユニットにより、前記整数ｎおよび前記整数ｍの少なくともいずれかが異なる結合係数候補を求め、前記結合係数候補による離散Ｂ－スプラインの線形結合近似と前記スペクトルとの誤差を評価して、前記結合係数を求める第２の演算ユニットをさらに有する、解析装置。
  5. 　請求項２ないし４のいずれかにおいて、前記第１の演算ユニットは、前記スペクトルに対して大きな正の内積を持つ離散Ｂ－スプラインを選んで逐次、線形結合近似で表現する際の結合係数を求めるユニットを含む、解析装置。
  6. 　請求項２ないし５のいずれかにおいて、前記整数ｎおよび前記整数ｍと、前記内積演算ユニットにおいて計算されるサンプリング数ＳＰｍａｘとは以下の関係を満たす、解析装置。
    ｎ×ｍ≦ＳＰｍａｘ
  7. 　請求項１ないし６のいずれかにおいて、前記スペクトルはイオン移動度センサーの出力を含む、解析装置。
  8. 　センサーから得られるスペクトルの解析方法であって、
    　前記スペクトルのサンプリング値と、前記スペクトルの第１の座標軸に沿って定義され、有限個のサンプリング点で正の値を持ち、他のサンプリング点では０である局所非負離散関数との内積を演算し、前記スペクトルを複数の局所非負離散関数の線形結合近似で表現する際の結合係数を求めることと、
    　前記結合係数に基づき前記スペクトルを前記複数の局所非負離散関数に対応する分布に分解することとを有し、
    　前記結合係数を求めることは、得られた結合係数が負の場合は結合係数を０に固定することを含む、解析方法。
  9. 　請求項８において、局所非負離散関数は、連続する整数ｎ個のサンプリング点では値が１となり、他のサンプリング点では値が０となる矩形離散関数の整数ｍ回畳み込み累積である離散Ｂ－スプラインを含み、
    　前記結合係数を求めることは、前記スペクトルと前記離散Ｂ－スプラインとの内積を演算することを含み、
    　前記分離することは、前記結合係数に基づき前記スペクトルを複数のＢ－スプライン分布に分解することを含む、解析方法。
  10. 　請求項９において、前記内積を演算することは、以下の式（Ａ）および（Ｂ）の伝達関数を演算することを含む、解析方法。
    （１／（１－ｚ-1））m・・・（Ａ）
    （１－ｚ-n）m　　　　・・・（Ｂ）
  11. 　請求項９または１０において、前記結合係数を求めることは、前記整数ｎおよび前記整数ｍの少なくともいずれかが異なる結合係数候補による離散Ｂ－スプラインの線形結合近似と前記スペクトルとの誤差を評価して、前記結合係数を求めることを含む、解析方法。
  12. 　請求項９ないし１１のいずれかにおいて、前記結合係数を求めることは、前記スペクトルに対して大きな正の内積を持つ離散Ｂ－スプラインを選んで逐次、線形結合近似で表現する際の結合係数を求めることを含む、解析方法。
  13. 　センサーから得られるスペクトルを解析するプログラムであって、
    　プロセッサが、前記スペクトルのサンプリング値と、前記スペクトルの第１の座標軸に沿って定義され、有限個のサンプリング点で正の値を持ち、他のサンプリング点では０である局所非負離散関数との内積を演算し、前記スペクトルを複数の局所非負離散関数の線形結合近似で表現する際の結合係数を求めることと、
    　前記結合係数に基づき前記スペクトルを前記複数の局所非負離散関数に対応する分布に分解することとを有し、
    　前記結合係数を求めることは、得られた結合係数が負の場合は結合係数を０に固定することを含む、プログラム。
PCT/JP2014/000730 2013-02-13 2014-02-13 解析装置 WO2014125819A1 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2013025894A JP2016075474A (ja) 2013-02-13 2013-02-13 解析装置
JP2013-025894 2013-02-13

Publications (1)

Publication Number Publication Date
WO2014125819A1 true WO2014125819A1 (ja) 2014-08-21

Family

ID=51353831

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/000730 WO2014125819A1 (ja) 2013-02-13 2014-02-13 解析装置

Country Status (2)

Country Link
JP (1) JP2016075474A (ja)
WO (1) WO2014125819A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20170066464A (ko) * 2014-10-14 2017-06-14 스미스 디텍션-워트포드 리미티드 방법 및 장치
WO2018066587A1 (en) * 2016-10-04 2018-04-12 Atonarp Inc. System and method for accurately quantifying composition of a target sample
CN107966428A (zh) * 2016-10-19 2018-04-27 西派特（北京）科技有限公司 一种提高微型拉曼光谱仪分辨率的方法
WO2018097129A1 (en) * 2016-11-23 2018-05-31 Atonarp Inc. System and method for determining set of mass to charge ratios for set of gases

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004286515A (ja) * 2003-03-20 2004-10-14 Olympus Corp 分光画像から複数の固有スペクトル成分の存在割合を求める方法
JP2005509135A (ja) * 2001-04-03 2005-04-07 ザ、テクサス、エイ、アンド、エム、ユーニヴァーサティ、システィム 周波数領域光子伝播計測により懸濁系内の粒子を特徴付けする方法
JP2009526359A (ja) * 2006-02-10 2009-07-16 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 照明装置の管理
US20090248320A1 (en) * 2008-03-27 2009-10-01 Nellcor Puritan Benett Llc System And Method For Unmixing Spectroscopic Observations With Nonnegative Matrix Factorization
JP2011524578A (ja) * 2008-06-20 2011-09-01 ナショナル・アイシーティ・オーストラリア・リミテッド 反射率スペクトルのコンパクト表現

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005509135A (ja) * 2001-04-03 2005-04-07 ザ、テクサス、エイ、アンド、エム、ユーニヴァーサティ、システィム 周波数領域光子伝播計測により懸濁系内の粒子を特徴付けする方法
JP2004286515A (ja) * 2003-03-20 2004-10-14 Olympus Corp 分光画像から複数の固有スペクトル成分の存在割合を求める方法
JP2009526359A (ja) * 2006-02-10 2009-07-16 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 照明装置の管理
US20090248320A1 (en) * 2008-03-27 2009-10-01 Nellcor Puritan Benett Llc System And Method For Unmixing Spectroscopic Observations With Nonnegative Matrix Factorization
JP2011524578A (ja) * 2008-06-20 2011-09-01 ナショナル・アイシーティ・オーストラリア・リミテッド 反射率スペクトルのコンパクト表現

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
BRABLEC A ET AL.: "Deconvolution of spectral line profiles: solution of th inversion problem", JOURNAL OF PHYSICS D APPLIED PHYSICS, vol. 32, no. 15, pages 1870 - 1875 *
ICHIGE K ET AL.: "A SIMPLE SCHEME OF DECOMPOSITION AND RECONSTRUCTION OF CONTINUOUS- TIME SIGNALS BY B-SPLINES", PROCEEDINGS OF THE 1998 IEEE INTERNATIONAL SYMPOSIUM ON CIRCUITS AND SYSTEMS, vol. 5, 31 May 1998 (1998-05-31), pages 453 - 456 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20170066464A (ko) * 2014-10-14 2017-06-14 스미스 디텍션-워트포드 리미티드 방법 및 장치
CN107078018A (zh) * 2014-10-14 2017-08-18 史密斯探测-沃特福特有限公司 方法和装置
JP2017536536A (ja) * 2014-10-14 2017-12-07 スミスズ ディテクション−ワトフォード リミテッド 方法及び装置
KR102551534B1 (ko) * 2014-10-14 2023-07-04 스미스 디텍션-워트포드 리미티드 방법 및 장치
WO2018066587A1 (en) * 2016-10-04 2018-04-12 Atonarp Inc. System and method for accurately quantifying composition of a target sample
US11227753B2 (en) 2016-10-04 2022-01-18 Atonarp Inc. System and method for accurately quantifying composition of a target sample
CN107966428A (zh) * 2016-10-19 2018-04-27 西派特（北京）科技有限公司 一种提高微型拉曼光谱仪分辨率的方法
CN107966428B (zh) * 2016-10-19 2020-01-03 西派特（北京）科技有限公司 一种提高微型拉曼光谱仪分辨率的方法
WO2018097129A1 (en) * 2016-11-23 2018-05-31 Atonarp Inc. System and method for determining set of mass to charge ratios for set of gases
JP2019536037A (ja) * 2016-11-23 2019-12-12 アトナープ株式会社 一組のガスに対応する一組の質量電荷比を求めるためのシステムおよび方法
US11183376B2 (en) 2016-11-23 2021-11-23 Atonarp Inc. System and method for determining set of mass to charge ratios for set of gases

Also Published As

Publication number Publication date
JP2016075474A (ja) 2016-05-12

Similar Documents

Publication Publication Date Title
CN104303258B (zh) 用于得到增强的质谱数据的方法和装置
US8346487B2 (en) Methods of automated spectral peak detection and quantification without user input
Gutierrez-Osuna et al. Transient response analysis of an electronic nose using multi-exponential models
Trygg O2‐PLS for qualitative and quantitative analysis in multivariate calibration
Golightly et al. Bayesian sequential inference for nonlinear multivariate diffusions
US8428889B2 (en) Methods of automated spectral peak detection and quantification having learning mode
EP2322922B1 (en) Method of improving the resolution of compounds eluted from a chromatography device
US9251122B2 (en) Method and apparatus for estimating a molecular mass parameter in a sample
Liu et al. Selective iteratively reweighted quantile regression for baseline correction
WO2014125819A1 (ja) 解析装置
JPH02297060A (ja) 時間変化信号データ分析装置及び方法と、パターン認識手段生成方法と、測定データのピーク分類方法及び装置と、データの属性の特性決定方法及び装置と、ピーク特性付け装置
CN103201621B (zh) 基于由传感器测量出的数据找出候选的方法和装置
US20120089342A1 (en) Methods of Automated Spectral and Chromatographic Peak Detection and Quantification without User Input
CN106529008B (zh) 一种基于蒙特卡罗及lasso的双集成偏最小二乘建模方法
Eiceman et al. Pattern recognition analysis of differential mobility spectra with classification by chemical family
Galvao et al. An application of subagging for the improvement of prediction accuracy of multivariate calibration models
US9989505B2 (en) Mass spectrometry (MS) identification algorithm
EP4102509A1 (en) Method and apparatus for identifying molecular species in a mass spectrum
JP2013181910A (ja) 質量分析システム
US20040254741A1 (en) Method and apparatus for modeling mass spectrometer lineshapes
Urban et al. Noise and baseline filtration in mass spectrometry
CN114364971B (zh) 用于确定样品中存在的组分的光谱装置和方法
Zolfonoun Determination of 6 Li and 7 Li abundances by inductively coupled plasma-optical emission spectrometry combined with derivative spectroscopy, multivariate curve resolution and multivariate calibration methods
Zhang et al. Denoising of Signals, Signal Enhancement, and Baseline Correction in Chromatographic Science
Gualdron et al. Novel feature selection method based on Stochastic Methods Coupled to Support Vector Machines using H-NMR data (data of olive and hazelnut oils)

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 14751830

Country of ref document: EP

Kind code of ref document: A1

DPE1 Request for preliminary examination filed after expiration of 19th month from priority date (pct application filed from 20040101)
NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14751830

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP