CN105760700A - Adaptive fixed-point IVA algorithm applicable to analysis on multi-subject complex fMRI data - Google Patents

Abstract

An adaptive fixed-point IVA algorithm suitable for multi-subject complex fMRI data analysis belongs to the field of biomedical signal processing. Using nonlinear functions based on multidimensional generalized Gaussian distribution (MGGD) to estimate the source vector component (SCV) distribution of complex fMRI data; using the maximum likelihood method to adaptively estimate the shape parameters of MGGD, automatically matching with the changing SCV distribution; in SCV-dominated The subspace update is based on the nonlinear function of MGGD to realize the denoising of the complex fMRI data; the pseudo-covariance matrix of the input data is added in the algorithm update process, and the non-circular characteristics of the complex fMRI data are directly used to further improve the IVA analysis of the complex fMRI data pertinence. The invention can effectively analyze complex fMRI data of multiple subjects with high noise level but the most comprehensive brain function information, and can provide more information for brain function research and brain disease diagnosis under the unfavorable conditions of large differences among subjects and low signal-to-noise ratio good basis.

Description

An Adaptive Fixed-Point IVA Algorithm Suitable for Multi-subject Complex fMRI Data Analysis

technical field

The invention relates to the field of biomedical signal processing, in particular to an analysis method for multi-subject complex functional magnetic resonance imaging (functionalmagneticresonanceimaging, fMRI) data.

Background technique

fMRI is an indispensable and powerful tool for current brain science research. Through the analysis of fMRI data, people can more comprehensively and deeply reveal brain function and brain mechanism. Raw acquired fMRI data is complex, including magnitude data and phase data. Among them, the noise of the phase data is relatively large, causing the noise of the complex fMRI data to be larger than that of the amplitude fMRI data. For this reason, phase fMRI data are discarded directly in most fMRI studies, and only magnitude fMRI data analysis is performed. However, more and more studies have shown that phase data contains a lot of unique and physiologically meaningful information, such as blood oxygenation levels during functional activation, the influence of large and small blood vessels, identification of different tissue types, and so on. Therefore, the analysis of complex fMRI data is helpful to advance the research of brain science.

Independent vector analysis (IVA) is a joint independent component analysis (ICA) method, which can solve the order ambiguity problem of single ICA method in multi-subject fMRI data analysis. Through the analysis of the IVA algorithm, multiple (usually 20-50) brain spatial activation area (spatial map, SM) components and their corresponding time course (time course, TC) components of each subject can be obtained from the fMRI data of multiple subjects. In this way, the difference information of the space-time response of each subject can be obtained, which can provide an important basis for the follow-up brain function research and clinical diagnosis.

IVA achieves SM signal separation by minimizing the mutual information between source vector components (source vector component, SCV, taking fMRI data analysis as an example, SCV contains the same SM components of each subject), while maximizing the correlation of SCV internal components . IVA algorithm usually selects multi-dimensional probability density function or nonlinear function to realize the estimation of SCV probability density distribution. Currently, IVA has real arithmetic and complex arithmetic. Among them, the real IVA algorithm has been used in the analysis of multi-subject amplitude fMRI data; the complex IVA algorithm has been used in the blind separation of speech and audio domains, but it has not been applied to the multi-subject complex fMRI data analysis. There are three main problems:

First, complex fMRI data are noisy. Due to the lack of denoising measures, the existing complex IVA algorithm cannot perform direct and effective analysis on multi-subject complex fMRI data.

Second, the SCV distribution of multi-subject complex fMRI data differs from that of frequency-domain speech. As a result, the multidimensional probability density function or nonlinear function in the existing complex IVA algorithm is only suitable for speech signals, not suitable for multi-subject complex fMRI data.

Third, the SCV distributions of multi-subject complex fMRI data varied widely. The reason for this is that complex fMRI data have a wide variety of SCV components, including task-related components, instantaneous task-related components, default network components, auditory components, visual components, motion components, scanner noise components, breathing noise components, head movement noise components, etc. Dozens of them. In this case, the distribution of each SCV component cannot be accurately estimated by using a fixed multidimensional probability density function or a nonlinear function.

Contents of the invention

The object of the present invention is to provide an adaptive fixed-point IVA algorithm suitable for multi-subject complex fMRI data analysis, by introducing fMRI data denoising, non-circular characteristics of complex fMRI data, and adaptively learning SCV distribution to solve the above-mentioned problem. Three aspects of the problem, in the IVA analysis of multi-subject complex fMRI data, the results are significantly better than the existing algorithms.

The technical solution of the present invention is to adopt a non-linear function based on multivariate generalized Gaussian distribution (MGGD) to estimate the SCV distribution of complex fMRI data; adopt the maximum likelihood estimation (maximumlikelihoodestimation, MLE) method to adaptively estimate the shape parameter of MGGD, In order to automatically match with the changing SCV distribution; update the nonlinear function based on MGGD in the SCV dominant subspace to realize the denoising of complex fMRI data; add the pseudo-covariance matrix of the input data in the algorithm update process to directly use the complex number The acyclic nature of fMRI data further improves the pertinence of IVA for complex fMRI data analysis. The specific implementation steps are as follows:

Step 1: Input multi-subject complex fMRI data k=1,...,K. Among them, K represents the number of subjects; J represents the number of whole brain scans in the time dimension; M represents the number of voxels in the brain in the space dimension.

The second step: perform PCA compression and whitening on the complex fMRI data X ^(k) of each subject. Let N be the number of SM and TC components of each subject. Compress and whiten the complex fMRI data X ^{(k) of} subject k as As a compressed matrix, N<J, reduce the time dimension from J to N, is a whitening matrix, such that has a variance of 1.

The third step: initialization. Randomly initialize the unmixing matrix k=1,...,K, acting on Calculate the initial value of the nth SCV component is the nth column of W ^(k) , n=1,...,N, yes In the mth column of , m=1,...,M, superscript T means transpose, and superscript H means conjugate transpose. For simplicity, m is ignored in the following, and y _n (m), and x ^(k) (m) are abbreviated as y _n , and x ^(k) . Using formula (1) to calculate the initial value of the cost function of the adaptive fixed-point IVA algorithm

Among them, E represents mathematical expectation, p(y _n ) is the probability density function of y _n ; G(·) is a real-valued nonlinear function; |·| represents the modulus value of complex numbers. Let _λn be the covariance matrix The largest eigenvalue of , is the eigenvector corresponding to λ _n , let Then G( ) in formula (1) adopts the nonlinear function based on MGGD shown in formula (2):

$\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&lsqb;</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; (</span>{<span\; class="notranslate">\; \&Sigma;</span>}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; K</span>}}{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&rsqb;</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

Wherein, β _n is the shape parameter of MGGD, and is also the shape parameter of each SCV distribution, the initialization β _n = β ₀ , β ₀ can be 0.4-0.5.

Step 4: Update the unmixing matrix W ^(k) . For each column of W ^(k) n=1,...,N, update according to formula (3):

$\begin{array}{c}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&LeftArrow;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>}{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \}</span>\\ <span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \}</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}\\ <span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; \}</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \}</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>}\end{array}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

The superscript * indicates conjugation, and G′(·) and G″(·) are the first and second derivatives of G(·) respectively:

${\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$

${\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$

After all the columns of W ^(k) are updated, the de-correlation operation shown in formula (6) is performed on W ^(k) :

W ^(k) ←(W ^(k) (W ^(k) ) ^H ) ^-12 W ^(k) (6)

Step 5: Estimate and update the shape parameter β _n of each SCV distribution by using the maximum likelihood estimation method based on Newton-Raphson optimization:

${\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; \&LeftArrow;</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; g</span>}\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; ;</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; \&part;</span>\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; g</span>}\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; ;</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; /</span><span\; class="notranslate">\; \&part;</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

where the maximum likelihood estimate of β _n is

In the formula, Γ(·) is the gamma function; the inverse of Σ _n is

Step 6: Calculate the cost function of this iteration. Record the number of iterations this time as iter, and use formula (1) to calculate the cost function of this iteration

Step 7: Judging the iteration termination condition. Calculate the cost function of this iteration The difference from the cost function of the previous iteration when When the absolute value of is less than the preset threshold ε or reaches the maximum number of iterations iter _max , the adaptive fixed-point IVA algorithm ends the iteration, otherwise returns to the fourth step. The value range of the preset threshold ε can be set as 10 ^-5 to 10 ^-6 , and the maximum number of iterations iter _max can be set as 500 to 2000.

Step 8: Calculate the N SM components of each subject and N TC components n=1,...,N, k=1,...,K, as follows:

${\mathrm{<span\; class="notranslate">\; Y</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; W</span>}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span>$

${\mathrm{<span\; class="notranslate">\; R</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; W</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$

in,

Step 9: Perform phase correction and SM phase denoising. Using the phase correction and denoising method in the article "Yu, M.C., Lin, Q.H., Kuang, L.D., Gong, X.F., Cong, F., Calhoun, V.D., 2015. ICAoffullcomplex-valued fMRI data using phase information of spatial maps. Journal of Neuroscience Methods 24, 75-91", first Apply each TC component for phase correction, and then remove voxels in the corresponding SM component whose phase range is outside -π/4 to π/4, to obtain the phase-denoised SM component.

Step 10: Output N SM components after phase denoising and N TC components after phase correction for each subject.

The effect and the benefit that the present invention reaches are, and existing fast fixed-point IVA algorithm (fastfixed-pointIVA, FIVA), acyclic fast fixed-point IVA algorithm (non-circularfastfixed-pointIVA, non-FIVA), IVA-GL algorithm (the algorithm First use the IVA-G algorithm based on multivariate Gaussian distribution (multivariateGaussiandistribution) for initialization, and then use the IVA_L algorithm based on multivariate Laplace distribution (multivariateLaplacedistribution) for separation), and further add SCV dominant subspace to FIVA and non-FIVA Compared with five IVA algorithms including FIVAs for denoising and non-FIVAs, the present invention has obvious advantages in the unfavorable situation of large differences among subjects and low signal-to-noise ratio. Among them, in the analysis of complex fMRI data collected under the task of knocking fingers by 16 subjects, the correlation coefficient of the SM component amplitude estimated by each algorithm and the prior reference signal is used as the performance index. Compared with the above five IVA algorithms, the present invention The performance of the estimated task-related components is improved by 9%-49%; the performance of the estimated default network DMN (default mode network) component is improved by 12%-21%. In addition, the SCV internal error rate (proportion of inconsistent components in SCV) of task-related components and DMN components estimated by the present invention is respectively reduced by 40%-65% and 71%-84% compared with the above five IVA algorithms. Therefore, the present invention can effectively analyze multi-subject complex fMRI data with high noise level but the most comprehensive brain function information, thereby providing better basis for brain function research and brain disease diagnosis.

Description of drawings

Accompanying drawing is the working flow chart of analyzing multi-subject complex fMRI data in the present invention.

detailed description

A specific embodiment of the present invention will be described in detail below in conjunction with the technical scheme and accompanying drawings.

There are 16 complex fMRI data collected under the finger-tapping task, that is, K=16. Each subject underwent J=165 scans, each scan obtained whole brain data of 53×63×46, and the number of voxels in the brain M=59610. Assuming that the number of SM and TC components of each subject is N=50, the steps of the multi-subject complex fMRI data analysis using the present invention are shown in the accompanying drawings.

Step 1: Input multi-subject complex fMRI data k=1,...,16.

The second step: perform PCA compression and whitening on the complex fMRI data X ^(k) of each subject. Compress and whiten the complex fMRI data X ^(k) for each subject as compressed array whitening array

The third step: initialization. Randomly initialize the unmixing matrix k=1,...,16, set the initial value of the shape parameter β _n of each SCV distribution β ₀ =0.4, n=1,...,50. Using formula (1) to calculate the initial value of the cost function of the adaptive fixed-point IVA algorithm

Step 4: Update the unmixing matrix W ^(k) . Use equation (3) to update all columns of the unmixing matrix W ^(k) n=1,...,50; using formula (6) to decorrelate W ^(k) .

Step 5: Use formula (7) to update the shape parameter β _n of each SCV distribution, where n=1, . . . , 50.

Step 6: Use formula (1) to calculate the cost function of this iteration

Step 7: Judging the iteration termination condition. The preset threshold ε=10- ⁶ , the maximum number of iterations iter _max =1000, calculate the cost function of this iteration The difference from the cost function of the previous iteration when When the absolute value of is less than 10- ⁶ or reaches the maximum number of iterations 1000, the algorithm iteration ends, otherwise return to the fourth step.

Step 8: Calculate the 50 SM components of each subject using formula (10) and formula (11) and 50 TC ingredients

Step 9: Perform phase correction and SM phase denoising. Firstly, each TC component is applied for phase correction, and then voxels in the corresponding SM component whose phase range is outside -π/4 to π/4 are removed to obtain a phase-denoised SM.

Step 10: Output 50 SM components after phase denoising and 50 TC components after phase correction for each subject.

Claims (1)

1. An adaptive fixed-point IVA algorithm suitable for multi-subject complex fMRI data analysis, using a nonlinear function based on multidimensional generalized Gaussian distribution MGGD to estimate the SCV distribution of complex fMRI data; using the maximum likelihood estimation MLE method to adaptively estimate The shape parameters of MGGD are automatically matched with the changing SCV distribution; the nonlinear function based on MGGD is updated in the SCV dominant subspace to realize the denoising of complex fMRI data; the pseudo-covariance matrix of the input data is added during the algorithm update process, Directly use the non-circular characteristics of complex fMRI data to further improve the pertinence of IVA for complex fMRI data analysis; including the following steps: Step 1: Input multi-subject complex fMRI data K represents the number of subjects; J represents the number of whole brain scans in the time dimension; M represents the number of voxels in the brain in the space dimension; Step 2: Perform PCA compression and whitening on the complex fMRI data X ^(k) of each subject; assuming that the number of SM and TC components of each subject is N, compress the complex fMRI data X ^{(k) of} subject k and whitened as As a compressed matrix, N<J, reduce the time dimension from J to N, is a whitening matrix, such that The variance of is 1; Step 3: Initialization; random initialization of the unmixing matrix Acting on Calculate the initial value of the nth SCV component is the nth column of W ^(k) , n=1,...,N, yes In the mth column of , m=1,...,M, superscript T means transpose, superscript H means conjugate transpose; ignore m, set y _n (m), and x ^(k) (m) are abbreviated as y _n , and x ^(k) ; Utilize formula (1) to calculate the initial value of the cost function of adaptive fixed-point IVA algorithm Among them, E represents mathematical expectation, p(y _n ) is the probability density function of y _n ; G( ) is a real-valued nonlinear function; |·| represents the modulus of complex numbers; let λ _n be the covariance matrix The largest eigenvalue of , is the eigenvector corresponding to λ _n , let In formula (1), G( ) adopts the nonlinear function based on MGGD shown in formula (2): $\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&lsqb;</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; K</span>}}{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&rsqb;</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$ In the formula, β _n is the shape parameter of MGGD, which is also the shape parameter of each SCV distribution, and the initialization β _n = β ₀ ; Step 4: Update the unmixing matrix W ^(k) ; for each column of W ^(k) Update according to formula (3): $\begin{array}{c}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&LeftArrow;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>}{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \}</span>\\ <span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \}</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}\\ <span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; \}</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; \{</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \}</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>}\end{array}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$ The superscript * indicates conjugation, and G′(·) and G″(·) are the first and second derivatives of G(·) respectively: ${\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$ ${\mathrm{<span\; class="notranslate">\; G</span>}}^{<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>}^{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$ After all the columns of W ^(k) are updated, the de-correlation operation shown in formula (6) is performed on W ^(k) : W ^(k) ←(W ^(k) (W ^(k) ) ^H ) ^-1/2 W ^(k) (6) Step 5: Estimate and update the shape parameter β _n of each SCV distribution by using the maximum likelihood estimation method based on Newton-Raphson optimization: ${\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; \&LeftArrow;</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; log</span>}\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; ;</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; \&part;</span>\mathrm{<span\; class="notranslate">\; log</span>}\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; ;</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; /</span><span\; class="notranslate">\; \&part;</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$ where the maximum likelihood estimate of β _n is In the formula, Γ(·) is the gamma function; the inverse of Σ _n is Step 6: Calculate the cost function of this iteration; record the number of iterations as iter, and use formula (1) to calculate the cost function of this iteration Step 7: Judging the termination condition of the iteration; calculating the cost function of this iteration The difference from the cost function of the previous iteration when When the absolute value of is less than the preset threshold ε or reaches the maximum number of iterations iter _max , the adaptive fixed-point IVA algorithm ends the iteration, otherwise returns to the fourth step; Step 8: Calculate the N SM components of each subject and N TC components as follows: ${\mathrm{<span\; class="notranslate">\; Y</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; W</span>}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span>$ ${\mathrm{<span\; class="notranslate">\; R</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; W</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$ in, Step 9: Perform phase correction and SM phase denoising; first apply each TC component for phase correction, and then remove voxels whose phase range is outside -π/4 to π/4 in the corresponding SM component to obtain the phase denoising SM components; Step 10: Output N SM components after phase denoising and N TC components after phase correction for each subject.