Electromagnetic Noise Suppression of Magnetic Resonance Sounding Combined with Data Acquisition and Multi-Frame Spectral Subtraction in the Frequency Domain

Tingting Lin 1,2,* , Xiaokang Yao 1,2 , Sijia Yu 1,2 and Yang Zhang 1,2,* 1 College of Instrumentation and Electrical Engineering, Jilin University, Changchun 130061, China; mexk16@mails.jlu.edu.cn (X.Y.); yusj18@mails.jlu.edu.cn (S.Y.) 2 Key Lab of Geo-Exploration Instrumentation of Ministry of Education, Jilin University, Changchun 130061, China * Correspondence: ttlin@jlu.edu.cn (T.L.); yang_zhang@jlu.edu.cn (Y.Z.)


Introduction
Magnetic resonance sounding (MRS) is a direct and nondestructive method in underground water detection. It provides water content distribution and estimates hydrogeological properties and has been widely applied in hydrogeological investigation [1][2][3]. However, the low signal-to-noise ratio (SNR) is always the main reason to restrict the application of this method [4,5]. The weak signal, which is only tens to thousands of nano-volts, is difficult to be detect. Moreover, the relaxation time of MRS signal is relatively short, i.e., a few hundred milliseconds, which requires higher sensitivity of the MRS receiver [6,7]. Therefore, obtaining the reliable MRS signal has become a key point in MRS measurement.
In order to improve the SNR of MRS signal, both the MRS instrument hardware [8,9] and the post-processing algorithm [10][11][12] have been researched. In terms of acquiring more accurate MRS signal, a series of de-noising methods are developed, such as the notch filter method [5], adaptive notch

Noise recording
Noisy MRS signal recording N1 S1 N2 S2 Ni Si The first cycle The i-th cycle The second cycle From Figure 1, we can see that the ambient noise is first collected, then the noisy MRS signal is recorded. For example, the MRS receiver first collects the noise at the measurement site (N1 in Figure 1), then the MRS transmitter emits excitation current, and then the MRS receiver works again to receive the noisy MRS signal (S1 in Figure 1). By now, a cycle of the recording is completed. Then, the second cycle of the recording and the third cycle of the recording are repeated one by one. The left part of Figure 1 shows the pure noise. These noises will provide a noise spectrum during spectrum subtraction. The right part of Figure 1 illustrates the noisy MRS signals, which are used to provide the spectrum of noisy MRS signals.
After the acquisition of the pure noise and the noisy MRS signal, the multi-frame spectrum subtraction (MFSS) for noise suppression can be further performed. There are two methods to be selected for noise reduction. One is to process each data samples with MFSS separately. The other is using MFSS to process the stacked data. We know that the statistical average method is effective in SNR improvement. Besides, for the noise data collected in a different time, the spectrum of the stacked noise data may better represent the spectrum distribution of the ambient noise. Therefore, we prefer to utilize the stacked data using in MFSS.
According to the above analysis, when using the DA-MFSS method for noise suppression, the receiver of the MRS instrument first collects the pure noise at measurement site, then the transmitter transmits pulse current for 40 ms, and then the receiver collects the noisy MRS signal. The work procedure is shown in Figure 2. From Figure 1, we can see that the ambient noise is first collected, then the noisy MRS signal is recorded. For example, the MRS receiver first collects the noise at the measurement site (N1 in Figure 1), then the MRS transmitter emits excitation current, and then the MRS receiver works again to receive the noisy MRS signal (S1 in Figure 1). By now, a cycle of the recording is completed. Then, the second cycle of the recording and the third cycle of the recording are repeated one by one. The left part of Figure 1 shows the pure noise. These noises will provide a noise spectrum during spectrum subtraction. The right part of Figure 1 illustrates the noisy MRS signals, which are used to provide the spectrum of noisy MRS signals.
After the acquisition of the pure noise and the noisy MRS signal, the multi-frame spectrum subtraction (MFSS) for noise suppression can be further performed. There are two methods to be selected for noise reduction. One is to process each data samples with MFSS separately. The other is using MFSS to process the stacked data. We know that the statistical average method is effective in SNR improvement. Besides, for the noise data collected in a different time, the spectrum of the stacked noise data may better represent the spectrum distribution of the ambient noise. Therefore, we prefer to utilize the stacked data using in MFSS.
According to the above analysis, when using the DA-MFSS method for noise suppression, the receiver of the MRS instrument first collects the pure noise at measurement site, then the transmitter transmits pulse current for 40 ms, and then the receiver collects the noisy MRS signal. The work procedure is shown in Figure 2.

Noise Suppression Procedure for DA-MFSS
After data acquisition, noise suppression by MFSS can be carried out, which mainly includes five processing stages. In order to illustrate the process procedure clearly, we conclude it and illustrate it in Figure 3.

Noise Suppression Procedure for DA-MFSS
After data acquisition, noise suppression by MFSS can be carried out, which mainly includes five processing stages. In order to illustrate the process procedure clearly, we conclude it and illustrate it in Figure 3.

Noise Suppression Procedure for DA-MFSS
After data acquisition, noise suppression by MFSS can be carried out, which mainly includes five processing stages. In order to illustrate the process procedure clearly, we conclude it and illustrate it in Figure 3.   In part A, the data collection and statistical average; In part B, calculating the spectrum of the stacked noise; In part C, calculating the spectrum of the stacked noisy MRS signal; In part D, multi-frame spectral subtraction of the above two spectrum; In part E, calculating the MRS signal in the time domain through the inverse transformation and the overlap-add method. The DA-MFSS mainly includes the following five steps.

1.
Stage 1, the pure noise and the noisy MRS signals are calculated by statistical average, as shown in part A of Figure 3.

2.
Stage 2, the stacked noise is used to calculate the spectrum value at each frequency, as shown in part B of Figure 3.

3.
Stage 3, the stacked noisy MRS signal is applied to calculate the spectrum value at each frequency, as shown in part C of Figure 3.

4.
Stage 4, according to the MFSS theory, the above two spectral values are brought into the calculation of the energy value at different frequencies, as shown in part D of Figure 3.

5.
Stage 5, based on inverse transformation formula and the overlap-add method, the MRS signal in the time domain after spectral subtraction is obtained, as shown in part E of Figure 3.

Discrete Fourier Transform of the Segmented Noisy MRS Signal
The stacked noisy MRS signal is expressed as discrete form s[n]. After adding window and framing, the i-th frame noisy signal is s i [m].
where i is the frame index, w[·] represents window function, H is the sliding size of the window function, m denotes the time index within each frame, m = 1, · · · , M, and M is the frame length. In order to show an excellent effect of noise suppression, here we list several common window functions for selection, as shown in Table 1.

Exponential function type
In order to compare the window shape in the time domain and the size of the main lobe and side lobe in the frequency domain more clearly, we compared the rectangular window, the Hamming window, the Blackman window and the Gaussian window, as shown in Figure 4. In order to compare the window shape in the time domain and the size of the main lobe and side lobe in the frequency domain more clearly, we compared the rectangular window, the Hamming window, the Blackman window and the Gaussian window, as shown in Figure 4.  From the comparison of the five windows in Figure 4, we can conclude that the main lobe of the rectangular window is narrowest, and its frequency resolution is the highest. The main lobe attenuation of the rectangular window is the smallest, but its side lobe attenuation is also small, which leads to serious energy leakage of the side lobe. The main lobe of the Hamming window is slightly wider than that of the rectangular window, so its frequency resolution is not as good as that of the rectangular window, but the side lobe attenuation of the Hamming window is larger than that of the rectangular window. The main lobe of the Hanning window is approximately equal to that of the Hamming window. The first side lobe of the Hanning window is bigger than the first side lobe of the Hamming window. The second side lobe attenuation of the Hanning window is close to that of the Hamming window. The third side lobe attenuation of the Hanning window is larger than that of the rectangular window, Hamming window and Gaussian window, but smaller than that of the Blackman window. The main lobe of the Blackman window is slightly wider than that of the Hamming window, that is, its frequency resolution is lower than that of the Hamming window. However, the side lobe attenuation of the Blackman window is the largest among the five kinds of windows, so its spectrum tailing is the smallest. In the figure, the main lobe of the Gaussian window is wider than that of the other four windows, and its frequency resolution is the lowest among the five windows. Nevertheless, the side lobe attenuation of the Gaussian window is larger than that of rectangular window and Hamming window, and smaller than that of the Blackman window and the Hanning window. Note that, when choosing a window function, we expect to choose a window function with a high-frequency resolution and small side lobe. However, due to the contradiction between the main lobe width and the attenuation of the side lobe, it is impossible to find the window function that can achieve both high-frequency resolution and very small side lobe. Therefore, when choosing the window function, according to the characteristics of actual data, we should choose an appropriate window function.
Based on the above discussion, we consider the tradeoff between the frequency resolution and the side lobe attenuation of window function, and choose the window function [] wm as a Hamming window. From the comparison of the five windows in Figure 4, we can conclude that the main lobe of the rectangular window is narrowest, and its frequency resolution is the highest. The main lobe attenuation of the rectangular window is the smallest, but its side lobe attenuation is also small, which leads to serious energy leakage of the side lobe. The main lobe of the Hamming window is slightly wider than that of the rectangular window, so its frequency resolution is not as good as that of the rectangular window, but the side lobe attenuation of the Hamming window is larger than that of the rectangular window. The main lobe of the Hanning window is approximately equal to that of the Hamming window. The first side lobe of the Hanning window is bigger than the first side lobe of the Hamming window. The second side lobe attenuation of the Hanning window is close to that of the Hamming window. The third side lobe attenuation of the Hanning window is larger than that of the rectangular window, Hamming window and Gaussian window, but smaller than that of the Blackman window. The main lobe of the Blackman window is slightly wider than that of the Hamming window, that is, its frequency resolution is lower than that of the Hamming window. However, the side lobe attenuation of the Blackman window is the largest among the five kinds of windows, so its spectrum tailing is the smallest. In the figure, the main lobe of the Gaussian window is wider than that of the other four windows, and its frequency resolution is the lowest among the five windows. Nevertheless, the side lobe attenuation of the Gaussian window is larger than that of rectangular window and Hamming window, and smaller than that of the Blackman window and the Hanning window. Note that, when choosing a window function, we expect to choose a window function with a high-frequency resolution and small side lobe. However, due to the contradiction between the main lobe width and the attenuation of the side lobe, it is impossible to find the window function that can achieve both high-frequency resolution and very small side lobe. Therefore, when choosing the window function, according to the characteristics of actual data, we should choose an appropriate window function.
Based on the above discussion, we consider the tradeoff between the frequency resolution and the side lobe attenuation of window function, and choose the window function w[m] as a Hamming window. where m is the point index in the window function and M is the window length. Furthermore, in the frequency domain, the spectrum S i (k) of one frame s i [m] of noisy MRS signal s[n] is obtained by discrete Fourier transform.
Here, K is the total number of frequency points of one frame s i [m] in discrete Fourier transform, and k is the frequency point index. From Equation (5), we can get the amplitude of the spectrum S i (k) at the k-th frequency point as follows: Simultaneously, according to Equation (5), the phase A i (k) of the spectrum S i (k) at the frequency point k in the i-th frame of noisy MRS data can be obtained.

Energy Calculation of Each Frame Noise Data in the Frequency Domain
After the noise dataset is stacked, the representative noise r[n] can be obtained, which can be expressed as follows: where h[n], p[n] and v[n] represent power-line harmonic noise, spike noise and random noise, respectively. Similarly, the typical noise r[n] is windowed and framed.
Furthermore, the spectrum R i (k) of the i-th frame noise r i [m] is obtained by discrete Fourier transform.
where K is the total number of frequency points of one frame r i [m] in discrete Fourier transform, and k denotes the frequency point index.
Assuming that the time length of the representative noise r[n] obtained by statistical average is IS, then the corresponding frame number is NIS: where the operator · describes solving the maximum integer corresponding to a real number. In order to explain the calculation principle of the number of frames when the representative noise data r[n] is divided into frames, we give Figure 5 for description. In this figure, we take the noise data with 25 sampling points as an example, in which the small rectangular block represents a sampling point, such as sampling point 0, sampling point 4 and sampling point 9, etc. If the window length M = 7 and the sliding step H = 4, then the noise data can be divided into five frames. On the other hand, we calculate NIS = 24−7 4 + 1 = 5 by Equation (11). According to Figure 5, we conclude that the actual total number of frames is equal to the calculated total number of frames. Therefore, when calculating the total number of frames, the total number of frames of the whole data can be calculated quickly through Equation (11). Electronics 2020, 9, x FOR PEER REVIEW 8 of 17 5, we conclude that the actual total number of frames is equal to the calculated total number of frames. Therefore, when calculating the total number of frames, the total number of frames of the whole data can be calculated quickly through Equation (11). According to Equations (10) and (11), the average energy value () Dk of the spectrum at the frequency point k of the representative noise [] rn is as follows:

MRS Signal Estimation by MFSS
where  is the judgment factor of spectral subtraction,  is the gain compensation factor.
According to Equation (13), the estimated amplitude | ( ) | i NMR k of MRS signal is obtained, and then a new spectrum can be constructed by combining the phase angle () i Ak estimated by Equation (7), for example, the spectrum of the i-th frame at frequency point k is as follows.
Furthermore, in the light of Equation (14), the spectral values of each frequency point of the i-th frame can be obtained, so that the estimated spectrum According to Equations (10) and (11), the average energy value D(k) of the spectrum at the frequency point k of the representative noise r[n] is as follows:

MRS Signal Estimation by MFSS
After obtaining the spectrum S i (k) of the noisy MRS signal s[n] and the average energy D(k) of the representative noise r[n], we can use MFSS to estimate the energy EN i (k) of the i-th frame data of MRS signal at the frequency point k.
where α is the judgment factor of spectral subtraction, β is the gain compensation factor.
According to Equation (13), the estimated amplitude NMR i (k) of MRS signal is obtained, and then a new spectrum can be constructed by combining the phase angle A i (k) estimated by Equation (7), for example, the spectrum of the i-th frame at frequency point k is as follows.
Furthermore, in the light of Equation (14), the spectral values of each frequency point of the i-th frame can be obtained, so that the estimated spectrum NMR i of the i-th frame can be acquired.
Then, through the fast Fourier inverse transform, we can get the i-th frame estimated signal nmr i with the length M in the time domain  (15), the estimated MRS signals are segmented, and they cannot be directly and simply spliced into a complete MRS signal. In order to obtain a complete MRS signal, we apply the overlap-add method to construct the complete MRS signal. To illustrate the overlap-add method to estimate the complete MRS signal, we draw Figure 6 to describe the signal construction. In the figure, according to Equations (13) and (14), the estimated spectrum NMR i of the i-th frame can be obtained, for example, the estimated spectrum NMR 1 , NMR 2 , NMR 3 , NMR 4 , etc. Then, the i-th frame estimated signal nmr i with the length M in the time domain can be further obtained through the fast Fourier inverse transform, such as the first frame estimated signal nmr 1 , the second frame estimated signal nmr 2 , the third frame estimated signal nmr 3 , the fourth frame estimated signal nmr 4 , etc. According to Equation (15), the estimated MRS signals are segmented, and they cannot be directly and simply spliced into a complete MRS signal. In order to obtain a complete MRS signal, we apply the overlap-add method to construct the complete MRS signal. To illustrate the overlap-add method to estimate the complete MRS signal, we draw Figure 6 to describe the signal construction. In the figure, according to Equations (13) and (14), the estimated spectrum i NMR of the i-th frame can be obtained, for example, the estimated spectrum In the same way, we can get Part C, Part D, Part E, Part F and Part G of the complete MRS signal.  Figure 6. Estimated MRS signal construction via the overlap-add method after multi-frame spectral subtraction and fast Fourier inverse transform.
From Figure 6, we can conclude that according to the calculated  From Figure 6, we can conclude that according to the calculated nmr 1 , nmr 2 , nmr 3 , · · · nmr i by Equation (15), we can construct the estimated complete MRS signal nmr through the overlap-add method.

Simulation Experiment
In order to verify the noise suppression effect of DA-MFSS, we utilize the synthetic noisy MRS data to conduct a simulation experiment.
The initial amplitude nmr 0 (q), the apparent relaxation time T * 2 (q), the Larmor frequency f Larmor and the initial phase angle ϕ 0 (q) of the ideal MRS signal are set as shown in Table 2. In the simulated noisy MRS signal, the ideal MRS signal and the noise is included. The parameters of noise added in the synthesis is shown in Table 3. In the simulation experiment, we set the recording time length to 256 ms and the sampling frequency f s = 25 kHz. In order to facilitate the comparison with the synthesized noisy MRS signal after noise suppression, this section shows the synthesis process of the noisy MRS signal in Figure 7. Figure 7a shows the ideal MRS signal. The ideal signal parameters are shown in Table 2. Figure 7b shows the Gaussian noise 1, and its parameters are shown in Table 3. Figure 7d shows the noisy MRS signal corrupted with noise 1, in which the SNR is −1.61 dB. To make the noise closer to the actual noise, Figure 7c shows Gaussian noise 2 with the mean value of −0.19 nV and the standard deviation of 61.45 nV, but it is different from noise 1 shown in Figure 7b

Simulation Experiment
In order to verify the noise suppression effect of DA-MFSS, we utilize the synthetic noisy MRS data to conduct a simulation experiment.
The initial amplitude 0 () nmr q , the apparent relaxation time * 2 () Tq , the Larmor frequency Larmor f and the initial phase angle 0 () q  of the ideal MRS signal are set as shown in Table 2. In the simulated noisy MRS signal, the ideal MRS signal and the noise is included. The parameters of noise added in the synthesis is shown in Table 3.  In the simulation experiment, we set the recording time length to 256 ms and the sampling frequency s f = 25 kHz. In order to facilitate the comparison with the synthesized noisy MRS signal after noise suppression, this section shows the synthesis process of the noisy MRS signal in Figure 7. Figure 7a shows the ideal MRS signal. The ideal signal parameters are shown in Table 2. Figure 7b shows the Gaussian noise 1, and its parameters are shown in Table 3. Figure 7d shows the noisy MRS signal corrupted with noise 1, in which the SNR is −1.61 dB. To make the noise closer to the actual noise, Figure 7c shows Gaussian noise 2 with the mean value of −0.19 nV and the standard deviation of 61.45 nV, but it is different from noise 1 shown in Figure 7b.  In the noise suppression process of DA-MFSS, we set the window function as Hamming window, the window length as M = 1200 (1200/ f s = 48 ms), and the frame shift as H = 600 (600/ f s = 24 ms). We choose the judgment factor α = 9.8 and the gain compensation factor β = 0.001. Simultaneously, the Gaussian noise 2 is selected as the noise of spectral subtraction. According to the noise suppression steps of MFSS in Figure 3, after noise elimination the SNR of MRS data reaches 15.27 dB. Furthermore, in order to compare the waveforms before and after MRS data denoising, Figure 8 is drawn to illustrate the noise suppression effect of MFSS. From the figure, we can conclude that the de-noised MRS signal is close to the ideal MRS signal, and the noise is eliminated. Moreover, the SNR is improved by 16.89 dB.
In the noise suppression process of DA-MFSS, we set the window function as Hamming window, the window length as M = 1200 (1200/ s f = 48 ms), and the frame shift as H = 600 (600/ s f = 24 ms). We choose the judgment factor  = 9.8 and the gain compensation factor  = 0.001.
Simultaneously, the Gaussian noise 2 is selected as the noise of spectral subtraction. According to the noise suppression steps of MFSS in Figure 3, after noise elimination the SNR of MRS data reaches 15.27 dB. Furthermore, in order to compare the waveforms before and after MRS data denoising, Figure 8 is drawn to illustrate the noise suppression effect of MFSS. From the figure, we can conclude that the de-noised MRS signal is close to the ideal MRS signal, and the noise is eliminated. Moreover, the SNR is improved by 16

Field Data Processing
In order to verify the performance of DA-MFSS, we utilized the JLMRS instrument [13] to carry out a field experiment in Shaoguo Town, Changchun City, Jilin Province, China. In the field experiment, the adopted data acquisition method is as described in Section 2, that is, in each acquisition cycle, the pure noise data are collected first, then the excitation pulse moment is transmitted, and then the noisy magnetic resonance signal is collected. In the field experiment, the adopted data acquisition method is as described in the second section, that is to collect the noise data first, then transmit the excitation pulse moment, and then collect the noisy MRS signal. The transmitting coil and receiving coil are both single turn square coil with 100 m × 100 m, the single transmitting current duration is 40 ms, and the single recording duration is 256 ms. After laying the coil and setting up the MRS instrument, we performed 32 stacks, and then average the recorded data to calculate the pure noise [] rn and the noisy MRS signal [] sn. According to the steps in Figure 3, through the statistical average for the data, the stacked noise is obtained, as is shown in Figure 9. Figure 9a shows the pure noise [] rn and Figure 9b shows the amplitude spectrum of the representative noise [] rn in the frequency domain. From the figure, it can be concluded that the electromagnetic noise from the measurement area is large and complex.

Field Data Processing
In order to verify the performance of DA-MFSS, we utilized the JLMRS instrument [13] to carry out a field experiment in Shaoguo Town, Changchun City, Jilin Province, China. In the field experiment, the adopted data acquisition method is as described in Section 2, that is, in each acquisition cycle, the pure noise data are collected first, then the excitation pulse moment is transmitted, and then the noisy magnetic resonance signal is collected. In the field experiment, the adopted data acquisition method is as described in the second section, that is to collect the noise data first, then transmit the excitation pulse moment, and then collect the noisy MRS signal. The transmitting coil and receiving coil are both single turn square coil with 100 m × 100 m, the single transmitting current duration is 40 ms, and the single recording duration is 256 ms. After laying the coil and setting up the MRS instrument, we performed 32 stacks, and then average the recorded data to calculate the pure noise r[n] and the noisy MRS signal s[n]. According to the steps in Figure 3, through the statistical average for the data, the stacked noise is obtained, as is shown in Figure 9. Figure 9a shows the pure noise r[n] and Figure 9b shows the amplitude spectrum of the representative noise r[n] in the frequency domain. From the figure, it can be concluded that the electromagnetic noise from the measurement area is large and complex. Similarly, the collected noisy MRS data is processed by stacking. Figure 10a shows the noisy MRS signal [] sn and Figure 10b shows the amplitude spectrum of the noisy MRS signal [] sn . According to the figure, we can conclude that a large amount of harmonic noise and random noise completely submerge the MRS signal, and the trend of the MRS signal cannot be found in the time domain. Similarly, the collected noisy MRS data is processed by stacking. Figure 10a shows the noisy MRS signal s[n] and Figure 10b shows the amplitude spectrum of the noisy MRS signal s[n]. According to the figure, we can conclude that a large amount of harmonic noise and random noise completely submerge the MRS signal, and the trend of the MRS signal cannot be found in the time domain. Similarly, the collected noisy MRS data is processed by stacking. Figure 10a shows the noisy MRS signal [] sn and Figure 10b shows the amplitude spectrum of the noisy MRS signal [] sn .
According to the figure, we can conclude that a large amount of harmonic noise and random noise completely submerge the MRS signal, and the trend of the MRS signal cannot be found in the time domain. Next, we analyze the selection of parameters in two cases to obtain a better noise reduction effect of spectral subtraction.

Select Short Frame Length and Small Frame Sliding Distance
In the noise reduction through MFSS, we still set the window function as a Hamming window, the short frame length as M = 1200 and the small frame shift as H = 600. Taking the representative noise [] rn as the noise used for spectral subtraction, we select the judgment factor  = 4.28 and the gain compensation factor  = 0.001. According to the steps of DA-MFSS in Figure 3, the standard deviation of the noisy MRS signal (shown in Figure 10) is reduced to 386.40 nV from 741.37 prior to multi-frame spectral subtraction.
The comparison before and after noise reduction is shown in Figure 11a. It can be seen from the figure that the noise is suppressed though MFSS. From Figure 11b, we can conclude that, through spectral subtraction, in addition to the random noise being suppressed, the harmonic noise at multiple power frequencies is also suppressed. Moreover, the noise around the Larmor frequency is also suppressed. Based on noise suppression by DA-MFSS, we utilize a band-pass filter Next, we analyze the selection of parameters in two cases to obtain a better noise reduction effect of spectral subtraction.

Select Short Frame Length and Small Frame Sliding Distance
In the noise reduction through MFSS, we still set the window function as a Hamming window, the short frame length as M = 1200 and the small frame shift as H = 600. Taking the representative noise r[n] as the noise used for spectral subtraction, we select the judgment factor α = 4.28 and the gain compensation factor β = 0.001. According to the steps of DA-MFSS in Figure 3, the standard deviation of the noisy MRS signal (shown in Figure 10) is reduced to 386.40 nV from 741.37 prior to multi-frame spectral subtraction.
The comparison before and after noise reduction is shown in Figure 11a. It can be seen from the figure that the noise is suppressed though MFSS. From Figure 11b, we can conclude that, through spectral subtraction, in addition to the random noise being suppressed, the harmonic noise at multiple power frequencies is also suppressed. Moreover, the noise around the Larmor frequency is also suppressed. Based on noise suppression by DA-MFSS, we utilize a band-pass filter with band-pass range of 2320 Hz to 2350 Hz to extract the envelope of MRS signal with the initial amplitude of 413.37 nV and the apparent relaxation time of 124.90 ms, as shown in Figure 11a. In the case of long frame length and large sliding distance, we still set the window function as the Hamming window, but the frame length is 2600 and the frame shift is 1300. Taking the pure noise [] rn in Figure 9 as the noise used for spectral subtraction, we select the judgment factor  =

Select Long Frame Length and Large Frame Sliding Distance
In the case of long frame length and large sliding distance, we still set the window function as the Hamming window, but the frame length is 2600 and the frame shift is 1300. Taking the pure noise r[n] in Figure 9 as the noise used for spectral subtraction, we select the judgment factor α = 9.29 and the gain compensation factor β = 0.0001. In the light of the steps in Figure 3, the standard deviation of the noisy MRS signal (shown in Figure 10) is reduced to 157.57 nV from 741.37 nV prior to MFSS.
The noise reduction result is shown in Figure 12 under the condition of long frame length and large frame sliding distance. It can be seen from the figure that the noise is greatly suppressed though MFSS. From Figure 12b, we can conclude that by the spectral subtraction, in addition to the random noise being suppressed, the harmonic noise at multiple power frequencies is also effectively suppressed. Moreover, the noise around the Larmor frequency is also suppressed. Based on noise suppression by DA-MFSS, we utilize a band-pass filter with a band-pass range of 2320 to 2350 Hz to extract the envelope of MRS signal with the initial amplitude of 408.31 nV and the apparent relaxation time of 141.97 ms, as shown in Figure 12a.

Select Long Frame Length and Large Frame Sliding Distance
In the case of long frame length and large sliding distance, we still set the window function as the Hamming window, but the frame length is 2600 and the frame shift is 1300. Taking the pure noise [] rn in Figure 9 as the noise used for spectral subtraction, we select the judgment factor  = 9.29 and the gain compensation factor  = 0.0001. In the light of the steps in Figure 3, the standard deviation of the noisy MRS signal (shown in Figure 10) is reduced to 157.57 nV from 741.37 nV prior to MFSS.
The noise reduction result is shown in Figure 12 under the condition of long frame length and large frame sliding distance. It can be seen from the figure that the noise is greatly suppressed though MFSS. From Figure 12b, we can conclude that by the spectral subtraction, in addition to the random noise being suppressed, the harmonic noise at multiple power frequencies is also effectively suppressed. Moreover, the noise around the Larmor frequency is also suppressed. Based on noise suppression by DA-MFSS, we utilize a band-pass filter with a band-pass range of 2320 to 2350 Hz to extract the envelope of MRS signal with the initial amplitude of 408.31 nV and the apparent relaxation time of 141.97 ms, as shown in Figure 12a. Compared with Figure 11 and Figure 12, it is concluded that, in the case of long window length, large frame shift and large the judgment factor, the better noise suppression effect can be obtained. For the convenience of observing the residual noise, we plot the amplitude spectrum after DA-MFSS in Figure 13.  Compared with Figures 11 and 12, it is concluded that, in the case of long window length, large frame shift and large the judgment factor, the better noise suppression effect can be obtained. For the convenience of observing the residual noise, we plot the amplitude spectrum after DA-MFSS in Figure 13. Compared with Figure 11 and Figure 12, it is concluded that, in the case of long window length, large frame shift and large the judgment factor, the better noise suppression effect can be obtained. For the convenience of observing the residual noise, we plot the amplitude spectrum after DA-MFSS in Figure 13.   According to Figure 13, by DA-MFSS, the random noise is effectively suppressed. Most of the power frequency harmonic noise is also suppressed; however, there are still some power frequency harmonic noise residues (such as frequency points p'3, p'7, p'13 and p '15). Specifically, the residual power frequency harmonic noise at 1751, 2402, 3155 and 3250 Hz has not been completely eliminated, but the residual harmonic noise is much smaller than the power frequency harmonic in Figure 12. According to Figure 13, by DA-MFSS, the random noise is effectively suppressed. Most of the power frequency harmonic noise is also suppressed; however, there are still some power frequency harmonic noise residues (such as frequency points p'3, p'7, p'13 and p '15). Specifically, the residual power frequency harmonic noise at 1751, 2402, 3155 and 3250 Hz has not been completely eliminated, but the residual harmonic noise is much smaller than the power frequency harmonic in Figure 12.

Discussion
According to the previous parts, we can conclude that the DA-MFSS can not only suppress the random noise, but also eliminate the interference of harmonic noise. This method provides a new idea to deal with the noise embedded in the MRS data, which can process multiple types of noise at one time in the noise processing process. It can be seen from the results of simulation experiment and field experiment that the proposed method based on DA-MFSS has an excellent de-noising effect on MRS data containing complex noise.
The data used in DA-MFSS are the stacking data. Since the noise used in stacking comes from different time periods, the spectral value obtained from the noise data by stacking at each frequency can better reflect and represent the noise spectrum distribution around the measurement site. Therefore, the statistical averaging is not only conducive to the elimination of random noise, but also conducive to multi-frame spectral subtraction.
In DA-MFSS, the difference of denoising effect is related to the frame length, the judgment factor α and the gain compensation factor β. When setting the long window length, the large judgment factor and the small gain compensation factor (for example, the window length M = 2600, the frame shift H = 1300, the judgment factor α = 9.29, the gain compensation factor β = 0.0001), the good noise elimination effect can be obtained, as shown in Figure 12. However, when setting short window length and small judgment factor (for example, the window length M = 1200, the judgment factor α = 4.28), the noise suppression effect is somewhat reduced, as shown in Figure 11. Therefore, when using DA-MFSS to suppress noise, it is necessary to reasonably choose a relatively large window length, large judgment factor and small gain compensation factor.

Conclusions
In this paper, focusing on the noise interference of noisy MRS signals, the noise and noisy data acquisition and the theory of spectral subtraction are studied, respectively. The noise-suppression method based on data acquisition and multi-frame spectral subtraction (DA-MFSS) is proposed to process the noisy MRS data.
The method first collects pure noise through the receiving coil in the measurement area, and then collects the noisy MRS signal after transmitting pulse moments. The pure noise and the noisy MRS signal acquisition are repeated for several times, and then the pure noise and the noisy signal are statistically averaged to obtain the representative data prior to MFSS. Then, the MFSS is used to suppress noise of the segmented noisy MRS data. Furthermore, the overlap-add method is utilized to synthesize the segmented estimated signals into the complete MRS signal. Finally, through the numerical simulation and noise suppression verification of on-site noisy MRS data, the SNR of the MRS data can be improved by 16.89 dB, and the method proposed in this paper can well suppress random noise and harmonic noise. This method has the advantage of suppressing multiple types of noise at the same time, making up for the shortcomings of the traditional noise-suppression method in MRS data processing.