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

Abstract

The invention discloses an adaptive fixed-point IVA algorithm applicable to analysis on multi-subject complex fMRI data, and belongs to the field of biomedical signal processing. The algorithm comprises the following steps: estimating an SCV distribution of complex fMRI data by adopting an MGGD-based nonlinear function; adaptively estimating a shape parameter of an MGGD by adopting a maximum likelihood estimation method, and automatically matching the shape parameter and the variable SCV distribution; updating the MGGD-based nonlinear function in an SCV-dominated subspace to implement noise elimination of the complex fMRI data; adding a pseudo-covariance matrix of input data in an algorithm updating process, and further improving the pertinence of IVA on the complex fMRI data by directly utilizing a non-circular characteristic of the complex fMRI data. According to the algorithm, multi-subject complex fMRI data of which the noise level is high but the brain function information is most comprehensive can be effectively analyzed, and under the unfavorable conditions of great differences among subjects and low signal to noise ratio, better bases can be provided for brain function researches and brain disease diagnosis.

Description

Self-adaptive fixed-point IVA algorithm suitable for multi-test complex fMRI data analysis

Technical Field

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

Background

fMRI is currently an indispensable powerful tool for brain science research. Through analysis of fMRI data, one can more fully and deeply reveal brain function and brain mechanics. The originally acquired fMRI data is complex, including amplitude data and phase data. In this case, the noise of the phase data is larger, resulting in a larger noise of the complex fMRI data than the amplitude fMRI data. For this reason, one discards phase fMRI data directly in most fMRI studies, and only performs amplitude fMRI data analysis. However, more and more research has shown that the phase data contains many unique and physiologically meaningful information, such as blood oxygenation levels during functional activation, the effect of large and small blood vessels, the identification of different tissue types, and so forth. Therefore, analyzing complex fMRI data helps to push brain science research into depth.

Independent Vector Analysis (IVA) is a combined Independent Component Analysis (ICA) method, and can solve the problem of order ambiguity of a single ICA method in multi-subject fMRI data analysis. Through analysis of the IVA algorithm, a plurality of (usually 20-50) brain space activation region (SM) components and corresponding time process (TC) components of each tested subject can be obtained from multi-tested fMRI data, so that difference information of each tested subject space-time response is obtained, and important basis is provided for subsequent brain function research and clinical diagnosis.

IVA achieves SM signal separation by minimizing mutual information between source vector components (SCVs, which contain the same SM components as each test, for example, in fMRI data analysis), while maximizing the correlation of the SCV internal components. The IVA algorithm typically chooses a multi-dimensional probability density function or a non-linear function to achieve an estimation of the SCV probability density distribution. Currently, IVA has real and complex algorithms. Wherein, the real IVA algorithm is used for multi-tested amplitude fMRI data analysis; the complex IVA algorithm has been used for blind separation of speech frequency domains, but has not been applied to multi-test complex fMRI data analysis, and the main problems are three:

first, complex fMRI data is noisy. Due to the lack of noise-canceling measures, the existing complex IVA algorithm cannot directly and effectively analyze the multi-test complex fMRI data.

Second, the SCV distribution of the multiple tested complex fMRI data is different from frequency domain speech. As a result, the multi-dimensional probability density function or the nonlinear function in the conventional complex IVA algorithm is only suitable for speech signals, and is not suitable for multi-test complex fMRI data.

Third, the SCV distribution of the multi-test complex fMRI data is greatly different. This is because SCV components of the complex fMRI data are various, and include dozens of kinds, such as task-related components, instantaneous task-related components, default network components, auditory components, visual components, motion components, scanner noise components, respiratory noise components, and head movement noise components. In this case, the distribution of the SCV components cannot be accurately estimated using a fixed multidimensional probability density function or a non-linear function.

Disclosure of Invention

The invention aims to provide an adaptive fixed-point IVA algorithm suitable for multi-test complex fMRI data analysis, which solves the three problems by introducing fMRI data noise elimination, the non-annular characteristic of the complex fMRI data and adaptive learning of SCV distribution, and obtains a result which is obviously better than that of the existing algorithm in the IVA analysis of the multi-test complex fMRI data.

The method adopts the technical scheme that a nonlinear function based on multi-dimensional generalized Gaussian distribution (MGGD) is adopted to estimate the SCV distribution of complex fMRI data; adaptively estimating the shape parameter of the MGGD by adopting a Maximum Likelihood Estimation (MLE) method, so as to be automatically matched with the SCV distribution which changes continuously; updating a non-linear function based on MGGD in an SCV main subspace to realize the noise elimination of complex fMRI data; and adding a pseudo-covariance matrix of input data in the algorithm updating process to directly utilize the non-annular characteristic of the complex fMRI data, and further improving the pertinence of the IVA to the analysis of the complex fMRI data. The method comprises the following concrete steps:

the first step is as follows: inputting multi-test complex fMRI dataK is 1, …, K. Wherein 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 spatial dimension.

The second step is that: for each test complex fMRI data X^(k)PCA compression and whitening were performed separately. Let the fraction of SM and TC tested be N. Complex fMRI data X of k to be tested^(k)Compressed and whitened to For compressing arrays, N<J, reducing the time dimension from J to N,to whiten the matrix, such thatThe variance of (a) is 1.

The third step: and (5) initializing. Randomly initializing a unmixing matrixK1, …, K, acting onCalculating to obtain the initial value of the nth SCV component Is W^(k)Column N, N is 1, …, N,is thatM-th column of (1), …, M, superscript T denoting transpose and superscript H denoting conjugate transpose. For simplicity, m is omitted below, and y is_n(m)、And x^(k)(m) is abbreviated to y_n、And x^(k). Calculating initial value of cost function of self-adaptive fixed point IVA algorithm by using formula (1)

Wherein E represents the mathematical expectation, p (y)_n) Is y_nA probability density function of; g (-) is a real-valued nonlinear function; | · | represents taking the magnitude of the complex number. Let λ_nIs a covariance matrixIs determined by the maximum characteristic value of the image,is λ_nCorresponding feature vector, orderThen G (-) in equation (1) takes the non-linear function based on MGGD shown in equation (2):

$G\left(q_n\right)={\left(q_n\right)}^{{\mathrm{\&beta;}}_n}={\&lsqb;{\mathrm{\&lambda;}}_n{\left({\&Sigma;}_{k=1}^K{v}_{kn}\right|y_n^{\left(k\right)}\left|\right)}^2\&rsqb;}^{{\mathrm{\&beta;}}_n}---\left(2\right)$

wherein, β_nFor the shape parameters of the MGGD, also the SCV distributions, β are initialized_n＝β₀，β₀It is preferably 0.4 to 0.5.

The fourth step: updating the demixing matrix W^(k). To W^(k)Each column ofN is 1, …, N, and is updated according to equation (3):

$\begin{array}{c}{w}_n^{\left(k\right)}\&LeftArrow;-E\left\{y_n^{\left(k\right)*}{G}^{\&prime;}\left(q_n\right)x^{\left(k\right)}\right\}\\ +E\{G^{\&prime;}\left(q_n\right)+{\left|y_n^{\left(k\right)}\right|}^2{G}^{\&prime;\&prime;}\left(q_n\right)\}{w}_n^{\left(k\right)}\\ +E\left\{x^{\left(k\right)}{\left(x^{\left(k\right)}\right)}^T\right\}E\left\{{\left(y_n^{\left(k\right)*}\right)}^2{G}^{\&prime;\&prime;}\left(q_n\right)\right\}{w}_n^{\left(k\right)*}\end{array}---\left(3\right)$

superscript denotes conjugation, and G' (. cndot.) and G "(. cndot.) are the first and second derivatives of G (. cndot.):

$G^{\&prime;}\left(q_n\right)={\mathrm{\&beta;}}_n{\left(q_n\right)}^{{\mathrm{\&beta;}}_n-1}---\left(4\right)$

$G^{\&prime;\&prime;}\left(q_n\right)={\mathrm{\&beta;}}_n({\mathrm{\&beta;}}_n-1){\left(q_n\right)}^{{\mathrm{\&beta;}}_n-2}---\left(5\right)$

W^(k)after all columns are updated, W is added^(k)The decorrelation operation shown in equation (6) is performed:

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

fifthly, estimating and updating the shape parameters β of each SCV distribution by adopting a maximum likelihood estimation method based on Newton-Raphson optimization_n：

${\mathrm{\&beta;}}_n\&LeftArrow;{\mathrm{\&beta;}}_n-\frac{logL({\mathrm{\&beta;}}_n;y_n(1),...,y_n(M\left)\right)}{\&part;logL({\mathrm{\&beta;}}_n;y_n(1),...,y_n(M\left)\right)/\&part;{\mathrm{\&beta;}}_n}---\left(7\right)$

Wherein, β_nIs estimated as

In the formula,(. is a gamma function; sigma_nIs inverse matrix of

And a sixth step: and calculating the cost function of the iteration. The iteration times are recorded as iter, and the cost function of the iteration is calculated by adopting the formula (1)

The seventh step: and judging an iteration termination condition. Calculating the cost function of the iterationDifference from last iteration cost functionWhen in useIs less than a preset threshold or reaches a maximum iteration number iter_maxAnd (4) ending the iteration by the self-adaptive fixed point IVA algorithm, and returning to the fourth step if the iteration is not ended. The value range of the preset threshold value can be set to 10^-5～10^-6Maximum number of iterations iter_maxIt is preferably 500 to 2000.

Eighth step: calculate the N SM components for each testAnd N TC componentsn＝1, …, N, K ═ 1, …, K, as follows:

$Y^{\left(k\right)}=W{\tilde{X}}^{\left(k\right)}---\left(10\right)$

$R^{\left(k\right)}={\left(U_2^{\left(k\right)}{U}_1^{\left(k\right)}\right)}^{-1}{\left(W^{\left(k\right)}\right)}^H---\left(11\right)$

wherein,

the ninth step: and carrying out phase correction and SM phase denoising. The phase correction and noise elimination method in the article "Yu, m.c., Lin, q.h., Kuang, l.d., Gong, x.f., Cong, f., Calhoun, v.d.,2015, icaoffulllcomplex-value mrIdatagushing phase information formation of spatial maps. journal of neurosciencemethod24, 75-91" is adopted, and each TC component is firstly applied to carry out phase correction, and then voxels with a phase range of-pi/4 to pi/4 in the corresponding SM component are removed to obtain the SM component with phase noise elimination.

The tenth step: and outputting N SM components after noise elimination of each tested phase and N TC components after phase correction.

Compared with the existing fast-fixed-point IVA (fixed-point IVA, FIVA) algorithm, non-ring fast-fixed-point IVA (non-ring fast-point IVA, non-FIVA) algorithm and IVA-GL algorithm (the algorithm is initialized by the IVA-G algorithm based on multi-dimensional Gaussian distribution and separated by the IVA-L algorithm based on multi-dimensional Laplace distribution), and five IVA algorithms including FIVAs and non-FIVAs, in which SCV (small-scale temporal noise reduction) is further added into the FIVA and the non-FIVA, the invention has obvious advantages under the adverse conditions of large difference between tests and low signal-to-noise ratio. In the analysis of complex fMRI data collected under the task of knocking the finger under 16 tests, the correlation coefficient between the SM component amplitude estimated by each algorithm and the prior reference signal is used as a performance index, and compared with the five IVA algorithms, the performance of the task related component estimated by the method is improved by 9-49 percent; the performance of the estimated default network DMN (default digital network) component is improved by 12-21%. In addition, the SCV internal error rate (the proportion of inconsistent components in SCV) of the estimated task-related components and the DMN components is reduced by 40-65% and 71-84% respectively compared with the five IVA algorithms. Therefore, the method can effectively analyze the multi-test complex fMRI data with high noise level but most comprehensive brain function information, and further provides better basis for brain function research and brain disease diagnosis.

Drawings

The FIGURE is a workflow diagram of the present invention for analyzing multiple test complex fMRI data.

Detailed Description

An embodiment of the present invention is described in detail below with reference to the accompanying drawings.

Existing 16 is tested for complex fMRI data acquired under the task of performing a tap on a finger, i.e., K16. Each test was performed 165J scans, each of which yielded 53 × 63 × 46 whole brain data, and the brain endosome number M was 59610. Assuming that the component number N of each SM and TC to be tested is 50, the procedure for performing multi-test complex fMRI data analysis using the present invention is shown in the drawing.

The first step is as follows: inputting multi-test complex fMRI datak＝1,…,16。

The second step is that: for each test complex fMRI data X^(k)PCA compression and whitening were performed separately. The complex fMRI data X of each tested object^(k)Compressed and whitened toCompression arrayWhitening array

The third step: and (5) initializing. Randomly initializing a unmixing matrixk 1, …,16, setting the shape parameter β of each SCV distribution_nInitial value β₀0.4, n is 1, …, 50. Calculating initial value of cost function of self-adaptive fixed point IVA algorithm by adopting formula (1)

The fourth step: updating the demixing matrix W^(k). Updating the demixing matrix W using equation (3)^(k)All columns ofn-1, …, 50; by adopting formula (6) to W^(k)A decorrelation is performed.

The fifth step is to update the shape parameter β of each SCV distribution by equation (7)_n，n＝1,…,50。

And a sixth step: calculating the cost function of the iteration by adopting the formula (1)

The seventh step: and judging an iteration termination condition. Preset threshold of 10-⁶Maximum number of iterations iter_max1000, calculating the cost function of the iterationDifference from last iteration cost functionWhen in useIs less than 10-⁶Or when the maximum iteration number is 1000, the algorithm iteration is ended, otherwise, the fourth step is returned.

Eighth step: the 50 SM components of each test piece were calculated by using equations (10) and (11)And 50 TC components

The ninth step: and carrying out phase correction and SM phase denoising. Firstly, each TC component is applied to carry out phase correction, and then voxels with the phase range outside-pi/4 in the corresponding SM component are removed to obtain the SM with the phase noise eliminated.

The tenth step: the 50 SM components after noise cancellation and the 50 TC components after phase correction for each phase to be tested are output.

Claims (1)

1. A self-adaptive fixed point IVA algorithm suitable for multi-test complex fMRI data analysis adopts a non-linear function based on multi-dimensional generalized Gaussian distribution MGGD to estimate SCV distribution of complex fMRI data; adaptively estimating the shape parameter of MGGD by adopting a Maximum Likelihood Estimation (MLE) method, and automatically matching with the SCV distribution which changes continuously; updating a non-linear function based on MGGD in an SCV main subspace to realize the noise elimination of complex fMRI data; adding a pseudo-covariance matrix of input data in the algorithm updating process, and further improving the pertinence of the IVA to the analysis of the complex fMRI data by directly utilizing the non-annular characteristic of the complex fMRI data; the method comprises the following steps:

the first step is as follows: inputting multi-test complex fMRI dataK represents the number of the tested samples; j represents the number of whole brain scans in the time dimension; m represents the number of voxels in the brain in the spatial dimension;

the second step is that: for each test complex fMRI data X^(k)Performing PCA compression and whitening respectively; let the component number of SM and TC of each test be N, and the complex fMRI data X of the test k^(k)Compressed and whitened to For compressing the matrix, N is less than J, the time dimension is reduced from J to N,to whiten the matrix, such thatThe variance of (a) is 1;

the third step: initializing; randomly initializing a unmixing matrixAct onCalculating to obtain the initial value of the nth SCV component Is W^(k)Column N, N is 1, …, N,is thatM-th column of (1), …, M, superscript T denoting transpose, superscript H denoting conjugate transpose; neglecting m, mixing y_n(m)、And x^(k)(m) is abbreviated to y_n、And x^(k)(ii) a Calculating initial value of cost function of self-adaptive fixed point IVA algorithm by using formula (1)

Wherein E represents the mathematical expectation, p (y)_n) Is y_nA probability density function of; g (-) is a real-valued nonlinear function; | · | represents taking the modulus of the complex number; let λ_nIs a covariance matrixIs determined by the maximum characteristic value of the image,is λ_nCorresponding feature vector, order In formula (1), G (-) is a non-linear function based on MGGD and shown in formula (2):

$G\left(q_n\right)={\left(q_n\right)}^{{\mathrm{\&beta;}}_n}={\&lsqb;{\mathrm{\&lambda;}}_n{\left({\mathrm{\&Sigma;}}_{k=1}^K{v}_{kn}\right|y_n^{\left(k\right)}\left|\right)}^2\&rsqb;}^{{\mathrm{\&beta;}}_n}---\left(2\right)$

in the formula, β_nFor the shape parameters of the MGGD, also the SCV distributions, β are initialized_n＝β₀；

The fourth step: updating the demixing matrix W^(k)(ii) a To W^(k)Each column ofUpdating is carried out according to the formula (3) respectively:

$\begin{array}{c}{w}_n^{\left(k\right)}\&LeftArrow;-E\left\{y_n^{\left(k\right)*}{G}^{\&prime;}\left(q_n\right)x^{\left(k\right)}\right\}\\ +E\{G^{\&prime;}\left(q_n\right)+{\left|y_n^{\left(k\right)}\right|}^2{G}^{\&prime;\&prime;}\left(q_n\right)\}{w}_n^{\left(k\right)}\\ +E\left\{x^{\left(k\right)}{\left(x^{\left(k\right)}\right)}^T\right\}E\left\{{\left(y_n^{\left(k\right)*}\right)}^2{G}^{\&prime;\&prime;}\left(q_n\right)\right\}{w}_n^{\left(k\right)*}\end{array}---\left(3\right)$

superscript denotes conjugation, and G' (. cndot.) and G "(. cndot.) are the first and second derivatives of G (. cndot.):

$G^{\&prime;}\left(q_n\right)={\mathrm{\&beta;}}_n{\left(q_n\right)}^{{\mathrm{\&beta;}}_n-1}---\left(4\right)$

$G^{\&prime;\&prime;}\left(q_n\right)={\mathrm{\&beta;}}_n({\mathrm{\&beta;}}_n-1){\left(q_n\right)}^{{\mathrm{\&beta;}}_n-2}---\left(5\right)$

W^(k)after all columns are updated, W is added^(k)The decorrelation operation shown in equation (6) is performed:

W^(k)←(W^(k)(W^(k))^H)^-1/2W^(k)(6)

fifthly, estimating and updating the shape parameters β of each SCV distribution by adopting a maximum likelihood estimation method based on Newton-Raphson optimization_n：

${\mathrm{\&beta;}}_n\&LeftArrow;{\mathrm{\&beta;}}_n-\frac{\log L({\mathrm{\&beta;}}_n;y_n(1),...,y_n(M\left)\right)}{\&part;\log L({\mathrm{\&beta;}}_n;y_n(1),...,y_n(M\left)\right)/\&part;{\mathrm{\&beta;}}_n}---\left(7\right)$

Wherein, β_nIs estimated as

In the formula,(. is a gamma function; sigma_nIs inverse matrix of

And a sixth step: calculating a cost function of the iteration; the iteration times are recorded as iter, and the cost function of the iteration is calculated by adopting the formula (1)

The seventh step: judging an iteration termination condition; calculating the cost function of the iterationDifference from last iteration cost functionWhen in useIs less than a preset threshold or reaches a maximum iteration number iter_maxIf so, the self-adaptive fixed point IVA algorithm finishes iteration, otherwise, the fourth step is returned;

eighth step: calculate the N SM components for each testAnd N TC components The following were used:

$Y^{\left(k\right)}=W{\tilde{X}}^{\left(k\right)}---\left(10\right)$

$R^{\left(k\right)}={\left(U_2^{\left(k\right)}{U}_1^{\left(k\right)}\right)}^{-1}{\left(W^{\left(k\right)}\right)}^H---\left(11\right)$

wherein,

the ninth step: carrying out phase correction and SM phase noise elimination; firstly, phase correction is carried out by using each TC component, and then voxels with the phase range outside-pi/4 in the corresponding SM component are removed to obtain the SM component with the phase noise eliminated;

the tenth step: the N SM components after noise elimination of each phase to be tested and the N TC components after phase correction are output.