Nuclear magnetic resonance spectrogram noise reduction method and system based on baseline recognition

Abstract

The invention discloses a baseline recognition-based nuclear magnetic resonance spectrogram noise reduction method, which comprises the following steps of: obtaining nuclear magnetic resonance spectrogram data y through a nuclear magnetic resonance spectrometer, and processing the nuclear magnetic resonance spectrogram data y by adopting a continuous wavelet transform method and a sliding window method in sequence to obtain a peak region and a baseline region of the nuclear magnetic resonance spectrogram data, the method comprises the following steps: identifying a peak region of nuclear magnetic resonance spectrogram data y, dividing the identified peak region into a positive peak region and a negative peak region, fitting a spectrogram baseline z according to the positive peak region and the negative peak region, acquiring a new baseline region and a new peak region from the nuclear magnetic resonance spectrogram data y according to the spectrogram baseline z, and performing noise reduction processing on the new baseline region by adopting a filtering method to obtain a noise reduction result. The baseline region and the peak region obtained after noise reduction jointly form nuclear magnetic resonance spectrogram data after noise reduction. The technical problems that an existing method based on matrix decomposition is high in calculation complexity, has the risk of misjudging a weak peak as noise and is poor in robustness can be solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of nuclear magnetic resonance data post-processing, and more specifically, relates to a nuclear magnetic resonance spectrum denoising method and system based on baseline identification. Background Art

[0002] Nuclear magnetic resonance (NMR) spectroscopy is a crucial component of analytical chemistry. Its high-resolution and non-invasive nature have made it an indispensable tool in fields such as structural biology, analytical chemistry, and metabolomics. However, due to its inherent population distribution, NMR spectroscopy exhibits inherently lower sensitivity compared to other characterization methods such as mass spectrometry and infrared spectroscopy. Consequently, various NMR spectroscopy noise reduction methods have emerged.

[0003] The existing NMR spectrum denoising methods mainly include the following: the first is a method based on deep learning (such as the DN-Unet deep neural network designed by Wu Ke et al.), which combines the encoder-decoder and convolutional neural network structure to design a deep neural network model, and directly suppresses the noise in the NMR spectrum by learning a large amount of data; the second is a method based on matrix decomposition (such as the Cadzow algorithm), which uses time domain data to form a Hankel or Toeplitz matrix, and uses singular value decomposition (SVD) to factorize the matrix, and achieves signal dimensionality reduction and denoising by discarding singular values ​​and singular vectors that are highly correlated with noise; the third is the traditional denoising method, which reduces noise by smoothing the spectral data.

[0004] However, the above-mentioned existing NMR spectroscopy noise reduction algorithms all have some non-negligible defects:

[0005] First, the training process of the above-mentioned deep learning-based methods is highly dependent on a large amount of high-quality measured NMR data. However, due to sample and instrument time limitations, it is very difficult to obtain sufficient and high-quality real data. Its actual application effect is limited by the quality of the training data, and it is difficult to achieve ideal noise reduction effects on complex spectra containing negative peaks, broad peaks, etc.

[0006] Second, the matrix decomposition-based method has high computational complexity and risks misclassifying weak peaks as noise, resulting in poor robustness.

[0007] Third, while suppressing noise, the above-mentioned traditional noise reduction methods will inevitably lose the details of information-rich peak features, thereby reducing the resolution of the NMR spectrum and making it difficult to simultaneously take into account signal integrity and noise suppression. Summary of the Invention

[0008] In response to the above defects or improvement needs of the prior art, the present invention provides a nuclear magnetic resonance spectrum denoising method and system based on baseline identification, which aims to accurately identify the peaks and baselines of the spectrum, retain the original rich peak information and effectively suppress the noise level, thereby solving the technical problems of high computational complexity, the risk of misjudging weak peaks as noise, and poor robustness of the existing matrix decomposition-based methods, and the training process of the existing deep learning-based methods is highly dependent on a large amount of high-quality measured NMR data. Due to sample and instrument time limitations, it is very difficult to obtain sufficient and high-quality real data, and its actual application effect is subject to the technical problem of the quality of the training data. In addition, the existing traditional denoising methods will inevitably lose the details of the information-rich peak features while suppressing noise, thereby reducing the resolution of the nuclear magnetic resonance spectrum and making it difficult to simultaneously take into account the signal integrity and noise suppression.

[0009] To achieve the above object, according to one aspect of the present invention, a method for reducing noise in nuclear magnetic resonance spectra based on baseline recognition is provided, comprising the following steps:

[0010] (1) obtaining nuclear magnetic resonance spectrum data y through a nuclear magnetic resonance instrument;

[0011] (2) processing the nuclear magnetic resonance spectrum data y obtained in step (1) by using a continuous wavelet transform method and a sliding window method to obtain a peak region and a baseline region of the nuclear magnetic resonance spectrum data, and dividing the identified peak region into a positive peak region and a negative peak region;

[0012] (3) fitting the spectrum baseline z according to the positive peak area and the negative peak area obtained in step (2);

[0013] (4) obtaining a new baseline region and a new peak region in the nuclear magnetic resonance spectrum data y obtained in step (1) according to the spectrum baseline z obtained in step (3);

[0014] (5) Using a filtering method to perform noise reduction processing on the new baseline region obtained in step (4), the baseline region obtained after the noise reduction processing and the peak region obtained in step (4) together constitute the nuclear magnetic resonance spectrum data after noise reduction.

[0015] Preferably, step (2) specifically includes the following sub-steps:

[0016] (2-1) Processing the NMR spectrum data y obtained in step (1) using a continuous wavelet transform method to obtain a first-order derivative spectrum, dividing the obtained first-order derivative spectrum into 32 regions, obtaining the noise standard deviation of each region, and taking the minimum value of all noise standard deviations as the noise level σ;

[0017] (2-2) Generate a threshold n*σ based on the noise level σ obtained in step (2-1), where n is a noise factor parameter (its value is preferably 8);

[0018] (2-3) Based on the sliding window method, a window with a length of β (β can be set to between 0.05% and 0.15% of the width of the nuclear magnetic resonance spectrum, preferably 0.15%) is continuously slid on the first-order derivative spectrum obtained in step (2-1), and the height difference between the maximum and minimum values ​​in the window is obtained each time the window is slid. If the height difference is greater than or equal to the threshold value n*σ set in step (2-2), the center point of the window is set as the peak region point; if it is less than the threshold value, the center point of the window is set as the baseline region point, thereby obtaining multiple peak region points and multiple baseline region points;

[0019] (2-4) Set counter m=1;

[0020] (2-5) Determine whether m is greater than the total number of peak region points obtained in step (2-3), if so, proceed to step (2-14), otherwise proceed to step (2-6);

[0021] (2-6) obtaining the mth peak region point from all peak region points obtained in step (2-3) as the current starting peak region point;

[0022] (2-7) Set counter M=1;

[0023] (2-8) Determine whether the first adjacent point of the current starting peak region point in the first-order derivative spectrum is a baseline region point or a peak region point. If it is a baseline region point, set m = m + 1 and return to step (2-5). If it is a peak region point, classify the current starting peak region point and its first adjacent point in the first-order derivative spectrum as the Mth peak region, and then proceed to step (2-9);

[0024] (2-9) Set counter j = 2;

[0025] (2-10) Determine whether m+j is greater than the total number of peak region points obtained in step (2-3). If so, proceed to step (2-14); otherwise, proceed to step (2-11);

[0026] (2-11) Determine whether the j-th adjacent point of the current starting peak region point in the first-order derivative spectrum is a baseline region point or a peak region point. If it is a baseline region point, set the j-1-th adjacent point of the current starting peak region point in the first-order derivative spectrum as the ending peak region point, and then proceed to step (2-12). If it is a peak region point, set the j-th adjacent point of the current starting peak region point in the first-order derivative spectrum as the M-th peak region, and then proceed to step (2-13);

[0027] (2-12) Set m=m+j, M=M+1, and return to step (2-5);

[0028] (2-13) Set j = j + 1 and return to step (2-9);

[0029] (2-14) Set counter count = 1;

[0030] (2-15) Determine whether count is less than or equal to the total number of peak regions M. If so, proceed to step (2-16), otherwise the process ends;

[0031] (2-16) Get the starting point y of the count-th peak area startIdx and the end point y endIdx , connect the starting point and the ending point into a straight line, obtain the value of each point in the count-th peak region on the straight line as the linear interpolation value corresponding to the point, obtain the difference between the value of the point on the nuclear magnetic resonance spectrum data y and the linear interpolation value as the peak height h corresponding to the point, and obtain the maximum peak height and the minimum peak height among the peak heights h corresponding to all points in the count-th peak region;

[0032] (2-17) Determine whether the maximum peak height in the count-th peak region obtained in step (2-16) is greater than 0 and the minimum peak height is greater than -n*σ, or the maximum peak height is less than 0 and the minimum peak height is less than -n*σ. If the maximum peak height is greater than 0 and the minimum peak height is greater than -n*σ, mark the count-th peak region as a positive peak region, and then proceed to step (2-18). If the maximum peak height is less than 0 and the minimum peak height is less than -n*σ, mark the count-th peak region as a negative peak region, and then proceed to step (2-18); otherwise, mark the count-th peak region as a mixed region (i.e., the peak region contains both positive peaks and negative peaks), and then proceed to step (2-19);

[0033] (2-18) Set count = count + 1 and return to step (2-15);

[0034] (2-19) Obtain the horizontal coordinate startIdx of the starting point and the horizontal coordinate endIdx of the ending point of the count-th peak region, and set the horizontal coordinate Index of the scanning point = startIdx (i.e., start scanning from the starting point of the peak region);

[0035] (2-20) Determine whether the horizontal coordinate Index of the scanning point is less than or equal to the horizontal coordinate endIdx of the end point. If so, proceed to step (2-21), otherwise return to step (2-18);

[0036] (2-21) Get the point with the horizontal coordinate Index in the count-th peak area, and get the peak height h of the point obtained in step (2-16) Index , and determine the peak height h Index Is it less than -n*σ, and is the point a minimum point in the nuclear magnetic resonance spectrum y? If so, set the horizontal coordinate of the left extension starting point LeftIndex = Index-1, the horizontal coordinate of the right extension starting point RightIndex = Index+1, and then go to step (2-22); otherwise, set the horizontal coordinate of the scanning point Index = Index+1, and return to step (2-20);

[0037] (2-22) Determine whether the horizontal coordinate LeftIndex of the left extension starting point is greater than or equal to the horizontal coordinate startIdx of the starting point. If so, proceed to step (2-23). ​​Otherwise, set the point with the horizontal coordinate startIdx in the count-th peak area as the negative peak starting point, and then proceed to step (2-24).

[0038] (2-23) Get the point with the horizontal coordinate LeftIndex in the count-th peak area, and get the peak height h of this point obtained in step (2-16) LeftIndex , and determine the peak height h LeftIndex Is it greater than 0? If so, set the point with the horizontal coordinate LeftIndex in the count-th peak area as the negative peak starting point, and then go to step (2-24); otherwise, set the horizontal coordinate LeftIndex of the left extension starting point = LeftIndex-1, and then return to step (2-22);

[0039] (2-24) Determine whether the horizontal coordinate RightIndex of the right extension starting point is less than or equal to the horizontal coordinate endIdx of the end point. If so, proceed to step (2-25); otherwise, set the point with the horizontal coordinate endIdx in the count-th peak area as the negative peak end point, and then return to step (2-18);

[0040] (2-25) Get the point with the horizontal coordinate RightIndex in the count-th peak area, and get the peak height h of this point obtained in step (2-16) RightIndex , and determine the peak height h RightIndex Is it greater than 0? If so, set the point with the horizontal coordinate RightIndex in the count-th peak area as the negative peak end point, set the horizontal coordinate Index of the scanning point = RightIndex+1, and return to step (2-20); otherwise, set the horizontal coordinate RightIndex of the right extension starting point = RightIndex+1, and then return to step (2-24).

[0041] Preferably, step (3) specifically includes the following sub-steps:

[0042] (3-1) Initialize the weight vector w to a column vector of length N and element 1, initialize the weight matrix W to a sparse matrix of size N*N with all diagonal elements being the weight vector w, and initialize the second-order difference matrix D to a size of N*N, where N represents the length of the NMR spectrum data y obtained in step (1);

[0043] (3-2) Set the number of iterations t = 1;

[0044] (3-3) Determine whether the number of iterations t is less than or equal to the preset maximum number of iterations maxIter. If so, proceed to step (3-4), otherwise the process ends;

[0045] (3-4) Obtain the spectrum baseline z based on the weight matrix W and second-order difference matrix D obtained in step (3-1):

[0046] z=(W+λD T D) -1 Wy

[0047] Where λ represents the smoothing parameter, and its value range is 1 to 10 13 , preferably 10 12 ;

[0048] (3-5) Create two arrays d and d of length N ― ;

[0049] (3-6) Set counter i=1;

[0050] (3-7) Determine whether the counter i is less than or equal to the length N of the nuclear magnetic resonance spectrum data y, if so, proceed to step (3-8), otherwise proceed to step (3-13);

[0051] (3-8) Get the point with horizontal coordinate i in the NMR spectrum and get its value y on the NMR spectrum data y i , get the point with horizontal coordinate i in the spectrum baseline z, and get its value z on the spectrum baseline z i , determine whether the point with abscissa i in the nuclear magnetic resonance spectrum is located in the positive peak region or the negative peak region obtained in step (2); if the point is located in the positive peak region, proceed to step (3-9); if the point is located in the negative peak region, proceed to step (3-10);

[0052] (3-9) Calculate the i-th element d in array d i =y i ―z i , the i-th element d in the array d iStore in array d and go to step (3-11);

[0053] (3-10) Calculate the i-th element d in array d i =z i ―y i , the i-th element d in the array d i Store in array d and go to step (3-11);

[0054] (3-11) Determine d i Is it less than 0? If so, set array d ― The i-th element in The array d ― The i-th element in Stored in array d ― and go to step (3-12), otherwise set array d ― The i-th element in The array d ― The i-th element in Store in array d ― and go to step (3-12);

[0055] (3-12) Set counter i=i+1 and return to step (3-7);

[0056] (3-13) Get array d ― Standard deviation

[0057] (3-14) Set counter j = 1;

[0058] (3-15) Determine whether the variable j is less than or equal to the length N of the nuclear magnetic resonance spectrum data y. If so, proceed to step (3-16); otherwise, proceed to step (3-19);

[0059] (3-16) Get the jth element d in array d j ;

[0060] (3-17) Use the j-th element d obtained in step (3-16) j and the array d obtained in step (3-13) ― Standard deviation Update the jth element value w in the weight vector w j :

[0061] (3-18) Set counter j = j + 1 and return to step (3-15);

[0062] (3-19)Judgment Is it less than the preset threshold coefficient ε? If so, set the number of iterations t=t+1 and return to step (3-3), otherwise the process ends; where norm represents the norm.

[0063] Preferably, in step (3-17), the jth element value w in the weight vector w is updated j The following formula is used:

[0064]

[0065] Preferably, step (4) is specifically to use a window with a length of γ to continuously slide on the spectrum baseline z obtained in step (3). If the baseline value z at the center point j of the window is j Located in [y min ,y max ], it indicates that the center point is located in the new baseline area, otherwise it indicates that the center point is located in the new peak area, where γ is between 0.05% and 0.15% of the width of the NMR spectrum, and y min and y max They represent the minimum and maximum values ​​of the NMR spectrum data y within the window respectively.

[0066] Preferably, the filtering method in step (5) can be a Wiener filtering method, a Gaussian filtering method, or a median filtering method.

[0067] Preferably, when the Gaussian filtering method is used, a Gaussian filter kernel with a window size of k and a central symmetry is used to calculate the Gaussian function x represents the distance between a point in the NMR spectrum within the filter window and the center of the Gaussian filter core. The standard deviation σ is a key parameter affecting the strength of the filter noise reduction. k is between 0.01% and 0.15% of the width of the NMR spectrum, preferably 0.03%.

[0068] According to another aspect of the present invention, a nuclear magnetic resonance spectrum denoising system based on baseline recognition is provided, comprising:

[0069] The first module is used to obtain nuclear magnetic resonance spectrum data y through a nuclear magnetic resonance instrument;

[0070] The second module is used to process the nuclear magnetic resonance spectrum data y obtained by the first module using a continuous wavelet transform method and a sliding window method to obtain the peak area and baseline area of ​​the nuclear magnetic resonance spectrum data, and divide the identified peak area into a positive peak area and a negative peak area;

[0071] The third module is used to fit the spectrum baseline z according to the positive peak area and the negative peak area obtained in the second module;

[0072] A fourth module is configured to obtain a new baseline region and a new peak region in the nuclear magnetic resonance spectrum data y obtained in the first module according to the spectrum baseline z obtained in the third module;

[0073] The fifth module is used to perform noise reduction processing on the new baseline area obtained by the fourth module using a filtering method. The baseline area obtained after the noise reduction processing and the peak area obtained by the fourth module together constitute the noise-reduced nuclear magnetic resonance spectrum data.

[0074] In general, the above technical solutions conceived by the present invention can achieve the following beneficial effects compared with the prior art:

[0075] (1) The present invention can reduce the dependence on a large amount of real training data and process spectra containing various peak types such as positive peaks and negative peaks: due to the use of (2-17) to step (2-25), it solves the problem of accurate identification of negative peak areas by setting appropriate peak type division conditions and adopting a forward and backward bidirectional scanning strategy to locate the starting and ending points of the negative peak area. It does not need to rely on a large amount of actual nuclear magnetic resonance experimental data for training, effectively reduces the stringent requirements on the quality and quantity of training data, and is more feasible in practical applications;

[0076] (2) The present invention can effectively avoid misjudgment of weak peaks and complex robustness verification: due to the use of steps (3-8) to step (4), it can identify new baseline regions and new peak regions in the nuclear magnetic resonance spectrum by differentially processing the positive peak region and the negative peak region, dynamically adjusting the weights, and fitting a more accurate spectrum baseline. This can more accurately distinguish between signals and noise, reduce the risk of misjudging weak peaks as noise, and reduce the need for extensive verification of the algorithm robustness under different experimental conditions;

[0077] (3) The present invention can maintain signal integrity and high resolution while effectively reducing noise: due to the adoption of step (5), by filtering and reducing noise only in the new baseline region, the rich information in the peak region is retained to the greatest extent, and both signal integrity and noise suppression can be taken into account at the same time, thereby achieving a good noise reduction effect while maintaining the original resolution of the nuclear magnetic resonance spectrum;

[0078] (4) The present invention has a small amount of calculation, a fast running speed, a good noise reduction effect on spectra containing various types of spectral peaks, and a wide adaptability. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 It is a flow chart of the nuclear magnetic resonance spectrum denoising method based on baseline identification of the present invention;

[0080] Figure 2 is the input NMR spectrum containing broad peaks;

[0081] Figure 3It is the baseline of the broad peak spectrum fitting in the present invention;

[0082] Figure 4 It is the output of the NMR spectrum with broad peaks after noise reduction;

[0083] Figure 5 is the input NMR spectrum containing negative peaks;

[0084] Figure 6 It is the baseline of the spectrum fitting containing negative peaks in the present invention;

[0085] Figure 7 It is the NMR spectrum with negative peaks after output noise reduction. DETAILED DESCRIPTION

[0086] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely for the purpose of explaining the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

[0087] The basic idea of ​​the present invention is as follows: first, nuclear magnetic resonance spectrum data is processed by continuous wavelet transform and sliding window method, so as to preliminarily identify the peak area and baseline area of ​​the spectrum, and divide the peak area into positive peak area and negative peak area; second, the baseline of the spectrum is fitted according to the divided positive and negative peak areas, so as to obtain an estimation of the baseline morphology, providing a basis for subsequent more accurate identification of the baseline area and peak area; then, the fitted baseline is used to re-divide the new baseline area and peak area through the sliding window method, so as to more accurately define the area in the nuclear magnetic resonance spectrum that needs noise reduction; finally, a filtering method is used to perform noise reduction on the new baseline area in the nuclear magnetic resonance spectrum, so as to finally obtain the nuclear magnetic resonance spectrum data after noise reduction.

[0088] like Figure 1 As shown, the present invention provides a method for reducing noise in nuclear magnetic resonance spectra based on baseline recognition, comprising the following steps:

[0089] (1) obtaining nuclear magnetic resonance spectrum data y through a nuclear magnetic resonance instrument;

[0090] (2) The nuclear magnetic resonance spectrum data y obtained in step (1) is processed by using a continuous wavelet transform method and a sliding window method to obtain a peak region and a baseline region of the nuclear magnetic resonance spectrum data, and the identified peak region is divided into a positive peak region and a negative peak region.

[0091] (3) fitting the spectrum baseline z according to the positive peak area and the negative peak area obtained in step (2);

[0092] (4) obtaining a new baseline region and a new peak region in the nuclear magnetic resonance spectrum data y obtained in step (1) according to the spectrum baseline z obtained in step (3);

[0093] Specifically, a window with a length of γ (γ can be set between 0.05% and 0.15% of the width of the NMR spectrum, preferably 0.15%) slides continuously on the spectrum baseline z obtained in step (3). If the baseline value z at the center point j of the window is j Located in [y min ,y max ] between (where y min and y max represent the minimum and maximum values ​​of the NMR spectrum data y in the window, respectively), then it means that the center point is located in the new baseline area, otherwise it means that the center point is located in the new peak area.

[0094] The advantage of this step is that the sliding window method is used to identify new baseline regions and new peak regions in the NMR spectrum based on the spectrum baseline, which can more accurately distinguish between signals and noise and reduce the risk of misjudging weak peaks as noise.

[0095] (5) Using a filtering method to perform noise reduction processing on the new baseline region obtained in step (4), the baseline region obtained after the noise reduction processing and the peak region obtained in step (4) together constitute the nuclear magnetic resonance spectrum data after noise reduction.

[0096] Specifically, the filtering method in this step can be Wiener filtering, Gaussian filtering, median filtering, etc.

[0097] When the Gaussian filter method is used, a Gaussian filter kernel with a window size of k (k can be set between 0.01% and 0.15% of the NMR spectrum width, preferably 0.03%) and a centrally symmetric Gaussian filter kernel is used to calculate the Gaussian function. x represents the distance between the point of the NMR spectrum within the filter window and the center of the Gaussian filter kernel, and the standard deviation σ is a key parameter affecting the strength of the filtering noise reduction.

[0098] The advantage of this step is that only the new baseline region is filtered and denoised, which maximizes the retention of rich information in the peak region and takes into account both signal integrity and noise suppression. This allows for good noise reduction while maintaining the original resolution of the NMR spectrum.

[0099] Step (2) specifically includes the following sub-steps:

[0100] (2-1) Processing the NMR spectrum data y obtained in step (1) using a continuous wavelet transform method to obtain a first-order derivative spectrum, dividing the obtained first-order derivative spectrum into 32 regions, obtaining the noise standard deviation of each region, and taking the minimum value of all noise standard deviations as the noise level σ;

[0101] (2-2) Generate a threshold n*σ based on the noise level σ obtained in step (2-1), where n is a noise factor parameter (its value is preferably 8);

[0102] (2-3) Based on the sliding window method, a window with a length of β (β can be set to between 0.05% and 0.15% of the width of the nuclear magnetic resonance spectrum, preferably 0.15%) is continuously slid on the first-order derivative spectrum obtained in step (2-1), and the height difference between the maximum and minimum values ​​in the window is obtained each time the window is slid. If the height difference is greater than or equal to the threshold value n*σ set in step (2-2), the center point of the window is set as the peak region point; if it is less than the threshold value, the center point of the window is set as the baseline region point, thereby obtaining multiple peak region points and multiple baseline region points;

[0103] (2-4) Set counter m=1;

[0104] (2-5) Determine whether m is greater than the total number of peak region points obtained in step (2-3), if so, proceed to step (2-14), otherwise proceed to step (2-6);

[0105] (2-6) obtaining the mth peak region point from all peak region points obtained in step (2-3) as the current starting peak region point;

[0106] (2-7) Set counter M=1;

[0107] (2-8) Determine whether the first adjacent point of the current starting peak region point in the first-order derivative spectrum is a baseline region point or a peak region point. If it is a baseline region point, set m = m + 1 and return to step (2-5). If it is a peak region point, classify the current starting peak region point and its first adjacent point in the first-order derivative spectrum as the Mth peak region, and then proceed to step (2-9);

[0108] (2-9) Set counter j = 2;

[0109] (2-10) Determine whether m+j is greater than the total number of peak region points obtained in step (2-3). If so, proceed to step (2-14); otherwise, proceed to step (2-11);

[0110] (2-11) Determine whether the j-th adjacent point of the current starting peak region point in the first-order derivative spectrum is a baseline region point or a peak region point. If it is a baseline region point, set the j-1-th adjacent point of the current starting peak region point in the first-order derivative spectrum as the ending peak region point, and then proceed to step (2-12). If it is a peak region point, set the j-th adjacent point of the current starting peak region point in the first-order derivative spectrum as the M-th peak region, and then proceed to step (2-13);

[0111] (2-12) Set m=m+j, M=M+1, and return to step (2-5);

[0112] (2-13) Set j = j + 1 and return to step (2-9);

[0113] (2-14) Set counter count = 1;

[0114] (2-15) Determine whether count is less than or equal to the total number M of peak regions. If so, proceed to step (2-16), otherwise the process ends.

[0115] (2-16) Get the starting point y of the count-th peak area startIdx and the end point y endIdx , connect the starting point and the ending point into a straight line, obtain the value of each point in the count-th peak region on the straight line as the linear interpolation corresponding to the point, obtain the difference between the value of the point on the nuclear magnetic resonance spectrum data y and the linear interpolation as the peak height h corresponding to the point, and obtain the maximum peak height and the minimum peak height among the peak heights h corresponding to all points in the count-th peak region.

[0116] (2-17) Determine whether the maximum peak height in the count-th peak region obtained in step (2-16) is greater than 0 and the minimum peak height is greater than -n*σ, or the maximum peak height is less than 0 and the minimum peak height is less than -n*σ. If the maximum peak height is greater than 0 and the minimum peak height is greater than -n*σ, mark the count-th peak region as a positive peak region, and then proceed to step (2-18). If the maximum peak height is less than 0 and the minimum peak height is less than -n*σ, mark the count-th peak region as a negative peak region, and then proceed to step (2-18); otherwise, mark the count-th peak region as a mixed region (i.e., the peak region contains both positive peaks and negative peaks), and then proceed to step (2-19).

[0117] (2-18) Set count = count + 1 and return to step (2-15).

[0118] (2-19) Get the horizontal coordinate startIdx of the starting point and the horizontal coordinate endIdx of the ending point of the count-th peak region, and set the horizontal coordinate Index of the scanning point = startIdx (ie, start scanning from the starting point of the peak region).

[0119] (2-20) Determine whether the horizontal coordinate Index of the scanning point is less than or equal to the horizontal coordinate endIdx of the end point. If so, proceed to step (2-21), otherwise return to step (2-18).

[0120] (2-21) Get the point with the horizontal coordinate Index in the count-th peak area, and get the peak height h of the point obtained in step (2-16) Index , and determine the peak height h Index Is it less than -n*σ, and is the point a minimum point in the nuclear magnetic resonance spectrum y? If so, set the horizontal coordinate of the left extension starting point LeftIndex = Index-1, the horizontal coordinate of the right extension starting point RightIndex = Index+1, and then go to step (2-22); otherwise, set the horizontal coordinate of the scanning point Index = Index+1, and return to step (2-20).

[0121] (2-22) Determine whether the horizontal coordinate LeftIndex of the left extension starting point is greater than or equal to the horizontal coordinate startIdx of the starting point. If so, proceed to step (2-23); otherwise, set the point with the horizontal coordinate startIdx in the count-th peak area as the negative peak starting point, and then proceed to step (2-24).

[0122] (2-23) Get the point with the horizontal coordinate LeftIndex in the count-th peak area, and get the peak height h of this point obtained in step (2-16) LeftIndex , and determine the peak height h LeftIndex Is it greater than 0? If so, set the point with the horizontal coordinate LeftIndex in the count-th peak area as the negative peak starting point, and then go to step (2-24); otherwise, set the horizontal coordinate LeftIndex of the left extension starting point to LeftIndex = LeftIndex-1, and then return to step (2-22).

[0123] (2-24) Determine whether the horizontal coordinate RightIndex of the right extension starting point is less than or equal to the horizontal coordinate endIdx of the end point. If so, proceed to step (2-25). Otherwise, set the point with the horizontal coordinate endIdx in the count-th peak region as the negative peak end point, and then return to step (2-18).

[0124] (2-25) Get the point with the horizontal coordinate RightIndex in the count-th peak area, and get the peak height h of this point obtained in step (2-16) RightIndex , and determine the peak height h RightIndex Is it greater than 0? If so, set the point with the horizontal coordinate RightIndex in the count-th peak area as the negative peak end point, set the horizontal coordinate Index of the scanning point = RightIndex+1, and return to step (2-20); otherwise, set the horizontal coordinate RightIndex of the right extension starting point = RightIndex+1, and then return to step (2-24).

[0125] The advantage of the above steps (2-17) to (2-25) is that there is no need to rely on a large amount of actual nuclear magnetic resonance experimental data for training, and the problem of accurate identification of negative peak regions is solved.

[0126] Step (3) specifically includes the following sub-steps:

[0127] (3-1) Initialize the weight vector w to a column vector of length N and element 1, initialize the weight matrix W to a sparse matrix of size N*N with all diagonal elements being the weight vector w, and initialize the second-order difference matrix D to a size of N*N, where N represents the length of the NMR spectrum data y obtained in step (1);

[0128] (3-2) Set the number of iterations t = 1;

[0129] (3-3) Determine whether the number of iterations t is less than or equal to the preset maximum number of iterations maxIter. If so, proceed to step (3-4), otherwise the process ends;

[0130] Specifically, the preset maximum number of iterations maxIter ranges from 20 to 200, preferably 200.

[0131] (3-4) Obtain the spectrum baseline z based on the weight matrix W and second-order difference matrix D obtained in step (3-1):

[0132] z=(W+λD T D) ―1 Wy

[0133] Where λ represents the smoothing parameter, and its value range is 1 to 10 13 , preferably 10 12 .

[0134] (3-5) Create two arrays d and d of length N ― ;

[0135] (3-6) Set counter i=1;

[0136] (3-7) Determine whether the counter i is less than or equal to the length N of the nuclear magnetic resonance spectrum data y, if so, proceed to step (3-8), otherwise proceed to step (3-13);

[0137] (3-8) Get the point with horizontal coordinate i in the NMR spectrum and get its value y on the NMR spectrum data y i , get the point with horizontal coordinate i in the spectrum baseline z, and get its value z on the spectrum baseline z i , determine whether the point with the horizontal coordinate i in the nuclear magnetic resonance spectrum is located in the positive peak area or the negative peak area obtained in step (2). If the point is located in the positive peak area, proceed to step (3-9); if the point is located in the negative peak area, proceed to step (3-10).

[0138] (3-9) Calculate the i-th element d in array d i =y i ―z i , the i-th element d in the array d i Store in array d and go to step (3-11);

[0139] (3-10) Calculate the i-th element d in array d i =z i ―y i , the i-th element d in the array d i Store in array d and go to step (3-11);

[0140] (3-11) Determine d i Is it less than 0? If so, set array d ― The i-th element in The array d ― The i-th element in Stored in array d ― and go to step (3-12), otherwise set array d ― The i-th element in The array d ― The i-th element in Store in array d ― and go to step (3-12);

[0141] (3-12) Set counter i=i+1 and return to step (3-7);

[0142] (3-13) Get array d ― Standard deviation

[0143] (3-14) Set counter j = 1;

[0144] (3-15) Determine whether the variable j is less than or equal to the length N of the nuclear magnetic resonance spectrum data y. If so, proceed to step (3-16); otherwise, proceed to step (3-19);

[0145] (3-16) Get the jth element d in array d j ;

[0146] (3-17) Use the j-th element d obtained in step (3-16) j and the array d obtained in step (3-13) ― Standard deviation Update the jth element value w in the weight vector w j :

[0147]

[0148] (3-18) Set counter j = j + 1 and return to step (3-15);

[0149] (3-19)Judgment Is it less than the preset threshold coefficient ε? If so, set the number of iterations t=t+1 and return to step (3-3), otherwise the process ends; where norm represents the norm;

[0150] Specifically, the preset threshold coefficient ε ranges from 0.0001 to 0.01, preferably 0.001.

[0151] The advantage of the above steps (3-8) to (3-17) is that by differentially processing the positive peak region and the negative peak region, the weight is dynamically adjusted to fit a more accurate nuclear magnetic resonance spectrum baseline, thereby providing a basis for more accurate identification of the baseline region and the peak region.

[0152] In summary, the NMR spectrum denoising method provided by the present invention accurately identifies peak and baseline regions by fitting the spectrum baseline, and then uses a filtering algorithm to denoise the resulting baseline region. Most algorithms tend to identify broad peak signals as baselines, but the present invention accurately identifies broad peak regions and can process spectra containing a variety of peak types, including positive and negative peaks, broad and narrow peaks. By filtering only the baseline region, the rich information in the peak region is retained while achieving excellent noise suppression.

[0153] It will be easily understood by those skilled in the art that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for reducing noise in nuclear magnetic resonance spectra based on baseline recognition, characterized in that: The steps include: (1) obtaining nuclear magnetic resonance spectrum data y through a nuclear magnetic resonance instrument; (2) processing the nuclear magnetic resonance spectrum data y obtained in step (1) by using a continuous wavelet transform method and a sliding window method to obtain a peak region and a baseline region of the nuclear magnetic resonance spectrum data, and dividing the identified peak region into a positive peak region and a negative peak region; (3) fitting the spectrum baseline z according to the positive peak area and the negative peak area obtained in step (2); (4) obtaining a new baseline region and a new peak region in the nuclear magnetic resonance spectrum data y obtained in step (1) according to the spectrum baseline z obtained in step (3); (5) Using a filtering method to perform noise reduction processing on the new baseline region obtained in step (4), the baseline region obtained after the noise reduction processing and the peak region obtained in step (4) together constitute the nuclear magnetic resonance spectrum data after noise reduction.

2. The method for reducing noise in nuclear magnetic resonance spectra based on baseline identification according to claim 1, characterized in that: Step (2) specifically includes the following sub-steps: (2-1) Processing the NMR spectrum data y obtained in step (1) using a continuous wavelet transform method to obtain a first-order derivative spectrum, dividing the obtained first-order derivative spectrum into 32 regions, obtaining the noise standard deviation of each region, and taking the minimum value of all noise standard deviations as the noise level σ; (2-2) Generate a threshold n*σ based on the noise level σ obtained in step (2-1), where n is a noise factor parameter (its value is preferably 8); (2-3) Based on the sliding window method, a window with a length of β (β can be set to between 0.05% and 0.15% of the width of the nuclear magnetic resonance spectrum, preferably 0.15%) is continuously slid on the first-order derivative spectrum obtained in step (2-1), and the height difference between the maximum and minimum values ​​in the window is obtained each time the window is slid. If the height difference is greater than or equal to the threshold value n*σ set in step (2-2), the center point of the window is set as the peak region point; if it is less than the threshold value, the center point of the window is set as the baseline region point, thereby obtaining multiple peak region points and multiple baseline region points; (2-4) Set counter m=1; (2-5) Determine whether m is greater than the total number of peak region points obtained in step (2-3), if so, proceed to step (2-14), otherwise proceed to step (2-6); (2-6) obtaining the mth peak region point from all peak region points obtained in step (2-3) as the current starting peak region point; (2-7) Set counter M=1; (2-8) Determine whether the first adjacent point of the current starting peak region point in the first-order derivative spectrum is a baseline region point or a peak region point. If it is a baseline region point, set m = m + 1 and return to step (2-5). If it is a peak region point, classify the current starting peak region point and its first adjacent point in the first-order derivative spectrum as the Mth peak region, and then proceed to step (2-9); (2-9) Set counter j = 2; (2-10) Determine whether m+j is greater than the total number of peak region points obtained in step (2-3). If so, proceed to step (2-14); otherwise, proceed to step (2-11); (2-11) Determine whether the j-th adjacent point of the current starting peak region point in the first-order derivative spectrum is a baseline region point or a peak region point. If it is a baseline region point, set the j-1-th adjacent point of the current starting peak region point in the first-order derivative spectrum as the ending peak region point, and then proceed to step (2-12). If it is a peak region point, set the j-th adjacent point of the current starting peak region point in the first-order derivative spectrum as the M-th peak region, and then proceed to step (2-13); (2-12) Set m=m+j, M=M+1, and return to step (2-5); (2-13) Set j = j + 1 and return to step (2-9); (2-14) Set counter count = 1; (2-15) Determine whether count is less than or equal to the total number M of peak regions. If so, proceed to step (2-16), otherwise the process ends. (2-16) Get the starting point y of the count-th peak area startIdx and the end point y endIdx , connect the starting point and the ending point into a straight line, obtain the value of each point in the count-th peak region on the straight line as the linear interpolation value corresponding to the point, obtain the difference between the value of the point on the nuclear magnetic resonance spectrum data y and the linear interpolation value as the peak height h corresponding to the point, and obtain the maximum peak height and the minimum peak height among the peak heights h corresponding to all points in the count-th peak region; (2-17) Determine whether the maximum peak height in the count-th peak region obtained in step (2-16) is greater than 0 and the minimum peak height is greater than -n*σ, or the maximum peak height is less than 0 and the minimum peak height is less than -n*σ. If the maximum peak height is greater than 0 and the minimum peak height is greater than -n*σ, mark the count-th peak region as a positive peak region, and then proceed to step (2-18). If the maximum peak height is less than 0 and the minimum peak height is less than -n*σ, mark the count-th peak region as a negative peak region, and then proceed to step (2-18); otherwise, mark the count-th peak region as a mixed region (i.e., the peak region contains both positive peaks and negative peaks), and then proceed to step (2-19); (2-18) Set count = count + 1 and return to step (2-15); (2-19) Obtain the horizontal coordinate startIdx of the starting point and the horizontal coordinate endIdx of the ending point of the count-th peak region, and set the horizontal coordinate Index of the scanning point = startIdx (i.e., start scanning from the starting point of the peak region); (2-20) Determine whether the horizontal coordinate Index of the scanning point is less than or equal to the horizontal coordinate endIdx of the end point. If so, proceed to step (2-21), otherwise return to step (2-18); (2-21) Get the point with the horizontal coordinate Index in the count-th peak area, and get the peak height h of the point obtained in step (2-16) Index , and determine the peak height h Index Is it less than -n*σ, and is the point a minimum point in the nuclear magnetic resonance spectrum y? If so, set the horizontal coordinate of the left extension starting point LeftIndex = Index-1, the horizontal coordinate of the right extension starting point RightIndex = Index+1, and then go to step (2-22); otherwise, set the horizontal coordinate of the scanning point Index = Index+1, and return to step (2-20); (2-22) Determine whether the horizontal coordinate LeftIndex of the left extension starting point is greater than or equal to the horizontal coordinate startIdx of the starting point. If so, proceed to step (2-23). ​​Otherwise, set the point with the horizontal coordinate startIdx in the count-th peak area as the negative peak starting point, and then proceed to step (2-24). (2-23) Get the point with the horizontal coordinate LeftIndex in the count-th peak area, and get the peak height h of this point obtained in step (2-16) LeftIndex , and determine the peak height h LeftIndex Is it greater than 0? If so, set the point with the horizontal coordinate LeftIndex in the count-th peak area as the negative peak starting point, and then go to step (2-24); otherwise, set the horizontal coordinate LeftIndex of the left extension starting point = LeftIndex-1, and then return to step (2-22); (2-24) Determine whether the horizontal coordinate RightIndex of the right extension starting point is less than or equal to the horizontal coordinate endIdx of the end point. If so, proceed to step (2-25); otherwise, set the point with the horizontal coordinate endIdx in the count-th peak area as the negative peak end point, and then return to step (2-18); (2-25) Get the point with the horizontal coordinate RightIndex in the count-th peak area, and get the peak height h of this point obtained in step (2-16) RightIndex , and determine the peak height h RightIndex Is it greater than 0? If so, set the point with the horizontal coordinate RightIndex in the count-th peak area as the negative peak end point, set the horizontal coordinate Index of the scanning point = RightIndex+1, and return to step (2-20); otherwise, set the horizontal coordinate RightIndex of the right extension starting point = RightIndex+1, and then return to step (2-24).

3. The nuclear magnetic resonance spectrum denoising method based on baseline identification according to claim 1 or 2, characterized in that: Step (3) specifically includes the following sub-steps: (3-1) Initialize the weight vector w to a column vector of length N and element 1, initialize the weight matrix W to a sparse matrix of size N*N with all diagonal elements being the weight vector w, and initialize the second-order difference matrix D to a size of N*N, where N represents the length of the NMR spectrum data y obtained in step (1); (3-2) Set the number of iterations t = 1; (3-3) Determine whether the number of iterations t is less than or equal to the preset maximum number of iterations maxIter. If so, proceed to step (3-4), otherwise the process ends; (3-4) Obtain the spectrum baseline z based on the weight matrix W and second-order difference matrix D obtained in step (3-1): z=(W+λD T D) -1 Wy Where λ represents the smoothing parameter, and its value range is 1 to 10 13 , preferably 10 12 ; (3-5) Create two arrays d and d of length N - ; (3-6) Set counter i=1; (3-7) Determine whether the counter i is less than or equal to the length N of the nuclear magnetic resonance spectrum data y, if so, proceed to step (3-8), otherwise proceed to step (3-13); (3-8) Get the point with horizontal coordinate i in the NMR spectrum and get its value y on the NMR spectrum data y i , get the point with horizontal coordinate i in the spectrum baseline z, and get its value z on the spectrum baseline z i , determine whether the point with abscissa i in the nuclear magnetic resonance spectrum is located in the positive peak region or the negative peak region obtained in step (2); if the point is located in the positive peak region, proceed to step (3-9); if the point is located in the negative peak region, proceed to step (3-10); (3-9) Calculate the i-th element d in array d i =y i -z i , the i-th element d in the array d i Store in array d and go to step (3-11); (3-10) Calculate the i-th element d in array d i =z i -y i , the i-th element d in the array d i Store in array d and go to step (3-11); (3-11) Determine d i Is it less than 0? If so, set array d - The i-th element in The array d - The i-th element in Stored in array d - and go to step (3-12), otherwise set array d - The i-th element in The array d - The i-th element in Store in array d - and go to step (3-12); (3-12) Set counter i=i+1 and return to step (3-7); (3-13) Get the standard deviation σ of array d- d -; (3-14) Set counter j = 1; (3-15) Determine whether the variable j is less than or equal to the length N of the nuclear magnetic resonance spectrum data y. If so, proceed to step (3-16); otherwise, proceed to step (3-19); (3-16) Get the jth element d in array d j ; (3-17) Use the j-th element d obtained in step (3-16) j and the standard deviation σ of the array d- obtained in step (3-13) d -Update the jth element value w in the weight vector w j : (3-18) Set counter j = j + 1 and return to step (3-15); (3-19)Judgment Is it less than the preset threshold coefficient ε? If so, set the number of iterations t=t+1 and return to step (3-3), otherwise the process ends; where norm represents the norm.

4. The method for reducing noise in nuclear magnetic resonance spectra based on baseline identification according to any one of claims 1 to 3, characterized in that: In step (3-17), update the j-th element value w in the weight vector w j The following formula is used:

5. The method for reducing noise in nuclear magnetic resonance spectra based on baseline identification according to claim 4, characterized in that: Step (4) is to use a window with a length of γ to slide continuously on the spectrum baseline z obtained in step (3). If the baseline value z at the center point j of the window is j Located in [y min ,y max ], it indicates that the center point is located in the new baseline area, otherwise it indicates that the center point is located in the new peak area, where γ is between 0.05% and 0.15% of the width of the NMR spectrum, and y min and y max They represent the minimum and maximum values ​​of the NMR spectrum data y within the window respectively.

6. The method for reducing noise in nuclear magnetic resonance spectra based on baseline identification according to claim 5, characterized in that: The filtering method in step (5) may be a Wiener filtering method, a Gaussian filtering method, or a median filtering method.

7. The method for reducing noise in nuclear magnetic resonance spectra based on baseline identification according to claim 6, characterized in that: When the Gaussian filtering method is used, a Gaussian filter kernel with a window size of k and a central symmetry is used to calculate the Gaussian function x represents the distance between a point in the NMR spectrum within the filter window and the center of the Gaussian filter core. The standard deviation σ is a key parameter affecting the strength of the filter noise reduction. k is between 0.01% and 0.15% of the width of the NMR spectrum, preferably 0.03%.

8. A nuclear magnetic resonance spectrum denoising system based on baseline recognition, characterized in that: include: The first module is used to obtain nuclear magnetic resonance spectrum data y through a nuclear magnetic resonance instrument; The second module is used to process the nuclear magnetic resonance spectrum data y obtained by the first module using a continuous wavelet transform method and a sliding window method to obtain the peak area and baseline area of ​​the nuclear magnetic resonance spectrum data, and divide the identified peak area into a positive peak area and a negative peak area; The third module is used to fit the spectrum baseline z according to the positive peak area and the negative peak area obtained in the second module; A fourth module is configured to obtain a new baseline region and a new peak region in the nuclear magnetic resonance spectrum data y obtained in the first module according to the spectrum baseline z obtained in the third module; The fifth module is used to perform noise reduction processing on the new baseline area obtained by the fourth module using a filtering method. The baseline area obtained after the noise reduction processing and the peak area obtained by the fourth module together constitute the noise-reduced nuclear magnetic resonance spectrum data.

Images