CN107358156B - Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform - Google Patents

Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform Download PDF

Info

Publication number
CN107358156B
CN107358156B CN201710417075.3A CN201710417075A CN107358156B CN 107358156 B CN107358156 B CN 107358156B CN 201710417075 A CN201710417075 A CN 201710417075A CN 107358156 B CN107358156 B CN 107358156B
Authority
CN
China
Prior art keywords
imf
hilbert
frequency
spectrum
feature
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201710417075.3A
Other languages
Chinese (zh)
Other versions
CN107358156A (en
Inventor
林春漪
曹文雄
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN201710417075.3A priority Critical patent/CN107358156B/en
Publication of CN107358156A publication Critical patent/CN107358156A/en
Application granted granted Critical
Publication of CN107358156B publication Critical patent/CN107358156B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/253—Fusion techniques of extracted features

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform, which comprises the following steps of: acquiring multiple frames of ultrasonic echo RF signals of a tissue sample, selecting a region of interest, and constructing an ultrasonic RF time sequence; performing feature extraction on the ultrasonic RF time sequence by using a HHT algorithm to obtain a sample feature vector based on the relation among time, frequency and energy; and (4) performing feature screening and feature fusion by using a rankfeatures algorithm, calculating a sample feature fusion index, and taking the feature fusion index as the feature of the characterization organization. The characteristic fusion index provided by the invention can effectively reflect the difference between different types of tissue samples, and is suitable for the field of characteristic extraction of ultrasonic tissue characterization.

Description

Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform
Technical Field
The invention relates to the technical field of ultrasonic tissue characterization, in particular to a feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform.
Background
Because the interaction mechanism of the ultrasonic wave after entering the biological tissue is not completely understood, people only can indirectly achieve the purpose of identifying the structure, the components and the state of the tissue by extracting the ultrasonic echo information and making an explanation, thereby promoting people to research the ultrasonic tissue characterization feature extraction.
The existing ultrasonic tissue characterization feature extraction methods mainly comprise three types: based on ultrasound B-mode maps, based on backscattered acoustic beam RF signals and based on ultrasound RF time series analysis. In the research based on the ultrasonic B-mode image, the texture characteristic of the gray level distribution of the B-mode image is mainly utilized, and the method uses the gray level of the ultrasonic image, so the method is easily influenced by imaging parameters of an ultrasonic diagnostic apparatus, TGC adjustment and the like; the method based on the back scattering sound beam RF signal needs depth attenuation compensation, the attenuation capacities of different tissues are different, and the individual difference can cause the sound wave propagation path to be different, so that the effect of the method is seriously influenced. And based on the method of the ultrasonic RF time sequence, the RF time sequence is derived from the RF data of different frames at the same position and the same depth, thereby reducing the influence of individual difference to a certain extent.
Conventional spectral analysis methods are mostly based on Fast Fourier Transform (FFT) and Wavelet Transform (WT). Due to the existence of frequency spectrum leakage and the barrier effect, the FFT analysis of the non-integer subharmonics has larger errors, and although the frequency spectrum leakage and the barrier effect can be better eliminated through the windowing interpolation algorithm, the algorithms usually reduce the frequency resolution. In the WT, since high-pass and low-pass decomposition filters are not perfectly ideal filters, there is a severe crossover phenomenon between respective frequency bands, and thus it is difficult to achieve strict division of respective frequency bands. The essence of the above methods is a theory based on the expansion of basis functions, and the signal analysis result largely depends on the experience of the designer.
Compared with the traditional time-frequency analysis method, Hilbert-Huang Transform (HHT for short) has obvious advantages in processing nonlinear non-stationary signals as a relatively emerging time-frequency analysis method, the HHT can adaptively decompose a series of components with different frequencies, namely Inherent Mode Functions (IMF), according to the characteristics of the signals, interference of other non-self factors such as selection of basis functions or window functions and the like is not considered, and the obtained result can better reflect the characteristics of the signals.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provide a method for extracting characteristics of ultrasonic tissue characterization based on Hilbert-Huang transform, the provided characteristic fusion index can effectively reflect the difference between different types of tissue samples, and the method is suitable for the field of characteristic extraction of ultrasonic tissue characterization.
The purpose of the invention is realized by the following technical scheme:
a feature extraction method based on Hilbert-Huang transform and applicable to ultrasonic tissue characterization comprises the following steps:
s1: acquiring ultrasonic echo RF signal multiframes of a tissue sample, selecting an ROI, and constructing an ultrasonic RF time sequence;
s2: performing feature extraction on the ultrasonic RF time sequence by using a HHT algorithm to obtain a sample feature vector based on the relation among time, frequency and energy;
s3: and (4) performing feature screening and feature fusion by using a rankfeatures algorithm, calculating a sample feature fusion index, and taking the feature fusion index as the feature of the characterization organization.
Preferably, in step S1, the specific process of constructing the ultrasonic RF time series includes:
s1-1: scanning tissues to obtain multiple frames of ultrasonic echo RF signals;
s1-2: demodulating a certain frame of ultrasonic echo RF signal and displaying a B-type image;
s1-3: selecting an ROI with the size of M multiplied by N on the B-type picture;
s1-4: taking the first L frames of RF signals for each point in the ROI constructs M N ultrasound RF time series with length L.
Preferably, in step S2, the method for calculating the sample feature vector by using the HHT algorithm includes:
s2-1: taking a time sequence x (t) in a sample ROI, and carrying out ensemble empirical mode decomposition to obtain n IMF components;
s2-2: extracting time domain features based on IMF components:
the number of IMF zero-crossing points is 5: IMF1-ZCs, IMF2-ZCs, IMF3-ZCs, IMF4-ZCs, IMF 5-ZCs;
IMF variance 5: IMF1-Var, IMF2-Var, IMF3-Var, IMF4-Var, IMF 5-Var;
IMF variance contribution rate of 5 IMF1-VarR, IMF2-VarR, IMF3-VarR, IMF4-VarR, IMF 5-VarR;
s2-3: respectively performing Hilbert spectrum analysis on the n-order IMFs of the time series x (t);
specifically, step S2-3 includes:
step a: hilbert transform of an IMF component c (t)
Figure BDA0001313958120000031
In which P is Cauchy main value, and then obtaining the analytic signal of c (t)
Figure BDA0001313958120000032
Is shown as polar sittingSubject form
Figure BDA0001313958120000033
Wherein the amplitude function
Figure BDA0001313958120000034
Function of energy
Figure BDA0001313958120000035
Phase function
Figure BDA0001313958120000036
The instantaneous frequency of c (t) can be calculated from the phase function
Figure BDA0001313958120000037
Step b: respectively carrying out Hilbert transformation on n-order IMFs of x (t) according to the step a, constructing an analytic signal, expressing the analytic signal in a polar coordinate form, and finally reconstructing the x (t) to obtain
Figure BDA0001313958120000038
The residual component r (t) is omitted here, Re representing the real part;
step c: balance
Figure BDA0001313958120000039
The Hilbert amplitude spectrum is called as a Hilbert spectrum for short, and the Hilbert marginal spectrum is obtained by integrating the Hilbert spectrum with time
Figure BDA00013139581200000310
In the formula
Figure BDA00013139581200000311
Represents the k-th order IMF component ck(t) Hilbert margin spectrum;
s2-4: extracting time domain-energy characteristics based on the amplitude function and the energy function:
5 IMF mean intensities: IMF1-AvgA, IMF2-AvgA, IMF3-AvgA, IMF4-AvgA, IMF5-AvgA,
5 IMF energies: IMF1-Egy, IMF2-Egy, IMF3-Egy, IMF4-Egy, IMF 5-Egy;
s2-5: extracting time domain-frequency features based on instantaneous frequency:
5 IMF highest frequencies: IMF1-MaxF, IMF2-MaxF, IMF3-MaxF, IMF4-MaxF, IMF 5-MaxF;
s2-6: extracting energy-frequency characteristics based on Hilbert marginal spectrum:
5 IMF mean center frequencies: IMF1-MCF, IMF2-MCF, IMF3-MCF, IMF4-MCF, IMF 5-MCF; 1 x (t) mean center frequency Orig-MCF;
s2-7: extracting frequency domain-energy characteristics based on Hilbert marginal spectrum:
hilbert marginal spectral entropy origin-EgyS of x (t);
the normalized Hilbert marginal spectrum of x (t) has low-band energy, medium-low band energy, medium-high band energy, high-band energy origin-MargS 1, origin-MargS 2, origin-MargS 3 and origin-MargS 4;
s2-8: extracting frequency domain-energy curve fitting characteristics based on Hilbert marginal spectrum:
normalized Hilbert marginal spectrum for x (t) straight-line fit slope, intercept: O-MLFSlope, O-MLFinterp;
slope and intercept of curve fitting of a first-order polynomial exponential function: O-MEFSlope, O-MEFInterp;
7 sixth order curve fitting coefficients: O-MSOFa0, O-MSOFa1, O-MSOFa2, O-MSOFa3, O-MSOFa4, O-MSOFa5, O-MSOFa 6;
s2-9: extracting time domain-frequency domain-energy characteristics based on the Hilbert spectrum:
statistical characteristics of Hilbert spectral gray statistical histogram of x (t) -mean TFImgmean, variance TFImgSD, skewness TFImgSkaew, kurtosis TFImgKurto;
s2-10: taking M × N time sequences x in sample ROI respectivelyi(t), where i is 0,1, 2., mxn-1, repeating steps S2-1 to S2-9 to obtain mxn eigenvectors corresponding to the mxn time series, and finally averaging the mxn eigenvectors to obtain one eigenvector, i.e., the sample eigenvector.
Preferably, n IMF components are obtained in step S2-1, where n is greater than or equal to 5.
Specifically, the specific process of step S2-1 includes:
step a: adding a gaussian white noise sequence to X (t) to obtain X (t) ═ X (t) + ω (t);
step b: determining all local maximum values and local minimum values of X (t), and obtaining an upper envelope line b and a lower envelope line b by adopting cubic spline interpolationmax(t) and bmin(t) calculating the average value
Figure BDA0001313958120000051
Find h1(t)＝X(t)-m1(t)；
Step c: if h is1(t) satisfies the condition of eigenmode function IMF, then h1(t) is the first IMF component of X (t); otherwise, repeating step b and h1(t) as a new X (t), calculate m11(t), judging again h11(t)＝h1(t)-m11(t) whether IMF conditions are met, if not, cycling k times until h1k(t)＝h1(k-1)(t)-m1k(t) satisfies the IMF condition. Note c1(t)＝h1k(t) is the first order IMF component of X (t), in which case the remaining component r1(t)＝X(t)-c1(t)；
Step d: will r is1(t) repeating steps a to c as new X (t) to obtain a second IMF-conditioned component c of X (t)2(t) cyclically executing the steps a to c until rn(t) becoming a monotonous function, and finishing the cycle, thereby obtaining a set of n (n is more than or equal to 5) order IMF components of the signal X (t);
step e: the steps a to d are repeatedly executed for N times, and different Gaussian white noise sequences omega are added each timei(t), decomposing to obtain a set of N groups of IMF components, and averaging the components of corresponding orders in the N groups of IMF components
Figure BDA0001313958120000052
As the final IMF component, ckAnd (t) represents a k-order IMF component obtained by EEMD decomposition of x (t).
Specifically, step S2-2 includes:
step a: for IMF component c of length Li(t) The method for counting the number of zero-crossing points is as follows: when j is more than or equal to 1 and less than or equal to L, if c (j) × c (j-1)<0, then it is marked as a zero crossing, and traverses [1, L ]]Finding IMF component ciAll zero crossings of (t) are obtainedi(t) zero-crossing points, so that zero-crossing points IMF1-ZCs, IMF2-ZCs, IMF3-ZCs, IMF4-ZCs and IMF5-ZCs of the first 5-order IMF can be obtained;
step b: according to the formula
Figure BDA0001313958120000053
Calculating the variances IMF1-Var, IMF2-Var, IMF3-Var, IMF4-Var and IMF5-Var of the first 5-order IMF, and then calculating the variances according to the formula
Figure BDA0001313958120000054
Variance contribution ratios of the first 5 th order IMF, IMF1-VarR, IMF2-VarR, IMF3-VarR, IMF4-VarR, and IMF5-VarR, were calculated.
Preferably, in step S2-4, after the step A of step S2-3 is performed on the time series x (t) of IMFs of order n, n amplitude functions a with length L can be obtainedi(t), energy function ei(t), (i ═ 0,1, 2.., n-1), according to the formula
Figure BDA0001313958120000061
Calculating the average intensity of the IMF of the first 5 th order, and then according to the formula
Figure BDA0001313958120000062
The energy of the first 5 th order IMF is calculated.
Preferably, in step S2-5, after the n-th order IMF of the time series x (t) is processed by step a of step S2-3, n instantaneous frequency functions f with length L can be obtainedi(t), (i ═ 0,1, 2.., n-1), according to the formula
Figure BDA0001313958120000063
The highest frequency of the first 5 IMF orders is calculated.
Preferably, in step S2-6, after the n-th order IMF of the time series x (t) is processed by steps a-c of step S2-3, n Hilbert marginal spectrums h containing m frequency point numbers can be obtainedi(f)，And x (t) Hilbert marginal spectrum h (f) containing m frequency points respectively according to a formula
Figure BDA0001313958120000064
And
Figure BDA0001313958120000065
wherein h isi(fj) Magnitude, f, of the j-th frequency point of the Hilbert marginal spectrum representing the i-th IMF componentjRepresenting the j-th frequency point, the mean center frequency of the first 5 IMF orders and the mean center frequency of x (t) are calculated.
Preferably, step S2-7 includes the steps of:
step a: according to Shannon entropy formula
Figure BDA0001313958120000066
Computing Hilbert marginal spectral entropy of x (t)
Figure BDA0001313958120000067
Wherein
Figure BDA0001313958120000068
Step b: normalizing the marginal spectrum h (f) of x (t) to obtain
Figure BDA0001313958120000069
Step c: h is to beo(f) Divided into 4 frequency bins: low frequency band
Figure BDA00013139581200000610
Middle and low frequency band
Figure BDA00013139581200000611
Middle and high frequency band
Figure BDA00013139581200000612
And high frequency band
Figure BDA00013139581200000613
According to the formula respectively
Figure BDA00013139581200000614
Figure BDA00013139581200000615
And
Figure BDA00013139581200000616
and calculating the energy of the normalized Hilbert marginal spectrum low frequency band, medium and high frequency band of x (t).
Preferably, the normalized Hilbert marginal spectrum h of x (t) in the step S2-8 is obtained by using a least square methodo(f) Performing straight line fitting to obtain O-MLFSlope and O-MLFinter, performing first order polynomial exponential function curve fitting to obtain O-MEFSlope and O-MEFinter, and performing sixth order curve fitting to obtain O-MSOFa0, O-MSOFa1, O-MSOFa2, O-MSOFa3, O-MSOFa4, O-MSOFa5 and O-OFMSa 6.
Preferably, step S2-9 specifically includes:
step a: after the n-order IMF of the time sequence x (t) is processed in the steps a to c of the step S2-3, a Hilbert amplitude spectrum H (f, t) of the time sequence x (t) can be obtained, the time t is taken as a horizontal axis, the frequency f is taken as a vertical axis, and the gray scale is controlled by using the amplitude to obtain a Hilbert spectrum time-frequency gray scale image of the time sequence x (t);
step b: then, a Hilbert spectrum time-frequency gray level image histogram h of the time sequence x (t) is obtained by taking the gray level of the image as a horizontal axis and the occurrence frequency of each gray level on the image as a vertical axisiWherein i is 0,1, 2. h isiExpressing the occurrence frequency of the ith gray level pixel of the gray histogram, and l expressing the gray level number of the time-frequency gray image;
step c: respectively using formulas
Figure BDA0001313958120000071
Figure BDA0001313958120000072
And
Figure BDA0001313958120000073
and calculating the mean, variance, skewness and kurtosis of the Hilbert spectrum gray statistical histogram of x (t).
Preferably, in step S3, the rankfeatures algorithm uses a toolbox function of MATLAB software, and the algorithm sorts the key features according to the category separation criterion, specifically using the function
[IDX,Z]＝rankfeatures(X,Group)，
Sorting the features using a binary-sorted independence evaluation criterion, the input X being a matrix, each column representing a vector of observations and the number of rows corresponding to the original number of features; the input Group comprises class labels, the output IDX is a list formed by row subscripts of X, the subscript list is arranged from high to low according to the significance of the features in the matrix X, and the output Z is the weight corresponding to each feature after the use criterion.
Preferably, in step S3, the sample feature fusion index is calculated by using a rankfeatures algorithm, and the specific steps are as follows:
s3-1: assigning a feature data matrix formed by all sample feature vectors to X, assigning vectors formed by all sample labels to Group, inputting the vectors into an [ IDX, Z ] ═ rankfeatures (X, Group) function for processing, and obtaining a feature weight vector Z containing 56 components;
s3-2: and selecting m characteristics with larger weights to carry out weighted summation to obtain a sample characteristic fusion index.
Compared with the prior art, the invention has the following advantages and beneficial effects:
1. the invention adopts Hilbert-Huang transform to process the ultrasonic RF time sequence, and extracts characteristics from the angle of the relation among time, frequency and energy.
2. The method utilizes a rankfeatures algorithm to carry out feature screening, calculates the weight of each feature, selects the feature with larger weight to carry out weighted summation, and constructs a feature fusion index.
3. The invention is based on the ultrasonic RF time sequence, because the signals of the ultrasonic RF time sequence are taken from the same position, the error caused by the difference of the ultrasonic wave propagation paths is reduced, and the deep attenuation compensation is not needed, thereby improving the characterization precision; has no special requirement on whether the tissue is uniform or not, and has wider application range.
4. The ultrasonic RF time sequence based on the invention can be obtained in the routine B ultrasonic examination without generating the additional hardware expense of the equipment.
Drawings
FIG. 1 is a flow chart of the method of example 2;
FIG. 2 is a pathological map of the breast tumor area of the control group and the treated group in example 2: FIG. 2(a) control group; FIG. 2(b) treatment group;
FIG. 3 is a type B plot of breast tumor areas of the control and treated groups of example 2: FIG. 3(a) control group; FIG. 3(b) treatment group;
FIG. 4 is the RF time series distribution diagram of the mammary gland of the rat of example 2.
Detailed Description
The present invention will be described in further detail with reference to examples and drawings, but the present invention is not limited thereto.
Example 1
A feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform comprises the following steps:
s1: collecting multiple frames of ultrasonic echo RF signals of a tissue sample, selecting a Region of interest (ROI for short), and constructing an ultrasonic RF time sequence;
s2: performing feature extraction on the ultrasonic RF time sequence by using a HHT algorithm to obtain a sample feature vector based on the relation among time, frequency and energy;
s3: and (4) performing feature screening and feature fusion by using a rankfeatures algorithm, calculating a sample feature fusion index, and taking the feature fusion index as the feature of the characterization organization.
In step S1, the specific process of constructing the ultrasonic RF time sequence is:
s1-1: scanning tissues to obtain multiple frames of ultrasonic echo RF signals;
s1-2: demodulating a certain frame of ultrasonic echo RF signal and displaying a B-type image;
s1-3: selecting an ROI with the size of M multiplied by N on the B-type picture;
s1-4: taking the first L frames of RF signals for each point in the ROI constructs M N ultrasound RF time series with length L.
In step S2, the method for calculating the sample feature vector using the HHT algorithm includes:
s2-1: taking a time sequence x (t) in a sample ROI, and carrying out ensemble empirical Mode Decomposition (EEMD for short) to obtain n (n is more than or equal to 5) IMF components, wherein the specific process is as follows:
step a: adding a gaussian white noise sequence to X (t) to obtain X (t) ═ X (t) + ω (t);
step b: determining all local maximum values and local minimum values of X (t), and obtaining an upper envelope line b and a lower envelope line b by adopting cubic spline interpolationmax(t) and bmin(t) calculating the average value
Figure BDA0001313958120000091
Find h1(t)＝X(t)-m1(t)；
Step c: if h is1(t) satisfies the condition of Intrinsic Mode Function (IMF), then h1(t) is the first IMF component of X (t); otherwise, repeating step b and h1(t) as a new X (t), calculate m11(t), judging again h11(t)＝h1(t)-m11(t) whether IMF conditions are met, if not, cycling k times until h1k(t)＝h1(k-1)(t)-m1k(t) satisfies the IMF condition. Note c1(t)＝h1k(t) is the first order IMF component of X (t), in which case the remaining component r1(t)＝X(t)-c1(t)；
Step d: will r is1(t) repeating steps a to c as new X (t) to obtain a second IMF-conditioned component c of X (t)2(t) cyclically executing the steps a to c until rn(t) becomes a monotonic function, and the cycle ends, fromObtaining a set of n (n is more than or equal to 5) order IMF components of the signal X (t);
step e: the steps a to d are repeatedly executed for N times, and different Gaussian white noise sequences omega are added each timei(t), decomposing to obtain a set of N groups of IMF components, and averaging the components of corresponding orders in the N groups of IMF components
Figure BDA0001313958120000101
As the final IMF component, ckAnd (t) represents a k-order IMF component obtained by EEMD decomposition of x (t).
S2-2: extracting temporal features (based on IMF components): IMF1-ZCs, IMF2-ZCs, IMF3-ZCs, IMF4-ZCs, IMF5-ZCs, IMF variance (5) IMF1-Var, IMF2-Var, IMF3-Var, IMF4-Var, IMF5-Var, IMF variance contribution rate (5) IMF1-VarR, IMF2-VarR, IMF3-VarR, IMF4-VarR and IMF5-VarR, wherein the specific process is as follows:
step a: for IMF component c of length Li(t), the method for counting the number of zero-crossing points is as follows: when j is more than or equal to 1 and less than or equal to L, if c (j) × c (j-1)<0, then it is marked as a zero crossing, and traverses [1, L ]]Finding IMF component ciAll zero crossings of (t) are obtainedi(t) zero-crossing points, so that zero-crossing points IMF1-ZCs, IMF2-ZCs, IMF3-ZCs, IMF4-ZCs and IMF5-ZCs of the first 5-order IMF can be obtained;
step b: according to the formula
Figure BDA0001313958120000102
Calculating the variances IMF1-Var, IMF2-Var, IMF3-Var, IMF4-Var and IMF5-Var of the first 5-order IMF, and then calculating the variances according to the formula
Figure BDA0001313958120000103
Calculating the variance contribution rates IMF1-VarR, IMF2-VarR, IMF3-VarR, IMF4-VarR and IMF5-VarR of the first 5-order IMF;
s2-3: hilbert spectrum analysis (HSA for short) is respectively carried out on the n-order IMFs of the time sequence x (t), and the specific process is as follows:
step a: hilbert transform of an IMF component c (t)
Figure BDA0001313958120000111
In which P is Cauchy main value, and then obtaining the analytic signal of c (t)
Figure BDA0001313958120000112
A (t) e in polar coordinate formiθ(t)Wherein the amplitude function
Figure BDA0001313958120000113
Function of energy
Figure BDA0001313958120000114
Phase function
Figure BDA0001313958120000115
The instantaneous frequency of c (t) can be calculated from the phase function
Figure BDA0001313958120000116
Step b: respectively carrying out Hilbert transformation on n-order IMFs of x (t) according to the step a, constructing an analytic signal, expressing the analytic signal in a polar coordinate form, and finally reconstructing the x (t) to obtain
Figure BDA0001313958120000117
The residual component r (t) is omitted here, Re representing the real part;
step c: balance
Figure BDA0001313958120000118
The Hilbert amplitude spectrum is called as a Hilbert spectrum for short, and the Hilbert marginal spectrum is obtained by integrating the Hilbert spectrum with time
Figure BDA0001313958120000119
In the formula
Figure BDA00013139581200001110
Represents the k-th order IMF component ck(t) Hilbert margin spectrum.
S2-4: extracting time-domain-energy features (based on magnitude and energy functions): IMF1-AvgA, IMF2-AvgA, IMF3-AvgA, IMF4-AvgA, IMF5-AvgA, IMF energy (5) IMF1-Egy, IMF2-Egy, IMF3-Egy, IMF4-Egy and IMF5-Egy, wherein the specific method is as follows:
after the n-th order IMF of the time series x (t) is processed by the step a of the step S2-3, n amplitude functions a with the length L can be obtainedi(t), energy function ei(t), (i ═ 0,1, 2.., n-1), according to the formula
Figure BDA00013139581200001111
Calculating the average intensity of the IMF of the first 5 th order, and then according to the formula
Figure BDA00013139581200001112
The energy of the first 5 th order IMF is calculated.
S2-5: time-frequency feature extraction (based on instantaneous frequency): the IMF with the highest IMF frequency (5) 1-MaxF, IMF2-MaxF, IMF3-MaxF, IMF4-MaxF and IMF5-MaxF specifically comprises the following steps:
after the n-th order IMF of the time series x (t) is processed by the step a of the step S2-3, n instantaneous frequency functions f with the length L can be obtainedi(t), (i ═ 0,1, 2.., n-1), according to the formula
Figure BDA0001313958120000121
The highest frequency of the first 5 IMF orders is calculated.
S2-6: extraction of energy-frequency features (based on Hilbert marginal spectrum): the method comprises the following steps of carrying out IMF1-MCF, IMF2-MCF, IMF3-MCF, IMF4-MCF and IMF5-MCF on the average central frequency (5), and carrying out x (t) origin-MCF on the average central frequency (1), wherein the specific method comprises the following steps:
after the n-order IMF of the time series x (t) is processed by the steps a to c of the step S2-3, n Hilbert marginal spectrums h containing m frequency point numbers can be obtainedi(f) And x (t) Hilbert marginal spectrum h (f) containing m frequency points according to a formula
Figure BDA0001313958120000122
And
Figure BDA0001313958120000123
wherein h isi(fj) Magnitude, f, of the j-th frequency point of the Hilbert marginal spectrum representing the i-th IMF componentjRepresenting the j-th frequency point, the mean center frequency of the first 5 IMF orders and the mean center frequency of x (t) are calculated.
S2-7: extracting frequency domain-energy features (based on Hilbert marginal spectrum): the method comprises the following steps of x (t) Hilbert marginal spectrum entropy (1) origin-EgyS, x (t) normalized Hilbert marginal spectrum low-frequency band energy, medium-low frequency band energy, medium-high frequency band energy, high-frequency band energy (4) origin-MargS 1, origin-MargS 2, origin-MargS 3 and origin-MargS 4, and specifically comprises the following steps:
step a: according to Shannon entropy formula
Figure BDA0001313958120000124
Computing Hilbert marginal spectral entropy of x (t)
Figure BDA0001313958120000125
Wherein
Figure BDA0001313958120000126
Step b: normalizing the marginal spectrum h (f) of x (t) to obtain
Figure BDA0001313958120000127
Step c: h is to beo(f) Divided into 4 frequency bins: low frequency band
Figure BDA0001313958120000128
Middle and low frequency band
Figure BDA0001313958120000129
Middle and high frequency band
Figure BDA0001313958120000131
And high frequency band
Figure BDA0001313958120000132
According to the formula respectively
Figure BDA0001313958120000133
Figure BDA0001313958120000134
And
Figure BDA0001313958120000135
and calculating the energy of the normalized Hilbert marginal spectrum low frequency band, medium and high frequency band of x (t).
S2-8: extracting frequency domain-energy curve fitting characteristics (based on Hilbert marginal spectrum): x (t) normalized Hilbert marginal spectrum straight-line fitting slope, intercept (2) O-MLFSlope, O-MLFinter, first order polynomial exponential function curve fitting slope, intercept (2) O-MEFSlope, O-MEFInterp, sixth order curve fitting coefficient (7) O-MSOFa0, O-MSOFa1, O-MSOFa2, O-MSOFa3, O-MSOFa4, O-MSOFa5 and O-MSOFa6, and the concrete method is as follows:
normalized Hilbert marginal spectrum h of x (t) by using least square methodo(f) Performing straight line fitting to obtain O-MLFSlope and O-MLFinter, performing first order polynomial exponential function curve fitting to obtain O-MEFSlope and O-MEFinter, and performing sixth order curve fitting to obtain O-MSOFa0, O-MSOFa1, O-MSOFa2, O-MSOFa3, O-MSOFa4, O-MSOFa5 and O-OFMSa 6.
S2-9: extracting time-frequency domain-energy features (based on Hilbert spectra): statistical characteristics of Hilbert spectrum gray statistical histogram of x (t), namely mean, variance, skewness, kurtosis (4), TFImgMean, TFImgSD, TFImgSkaew and TFImgKurto, comprise the following steps:
step a: after the n-order IMF of the time sequence x (t) is processed in the steps a to c of the step S2-3, a Hilbert amplitude spectrum H (f, t) of the time sequence x (t) can be obtained, the time t is taken as a horizontal axis, the frequency f is taken as a vertical axis, and the gray scale is controlled by using the amplitude to obtain a Hilbert spectrum time-frequency gray scale image of the time sequence x (t);
step b: then, a Hilbert spectrum time-frequency gray level image histogram h of the time sequence x (t) is obtained by taking the gray level of the image as a horizontal axis and the occurrence frequency of each gray level on the image as a vertical axisiWherein i is 0,1, 2. h isiRepresenting a grey scaleThe number of times of the ith gray level pixel appears, and l represents the gray level number of the time-frequency gray level image;
step c: respectively using formulas
Figure BDA0001313958120000136
Figure BDA0001313958120000141
And
Figure BDA0001313958120000142
and calculating the mean, variance, skewness and kurtosis of the Hilbert spectrum gray statistical histogram of x (t).
S2-10: taking M × N time sequences x in sample ROI respectivelyi(t), where i is 0,1, 2., mxn-1, repeating steps S2-1 to S2-9 to obtain mxn eigenvectors corresponding to the mxn time series, and finally averaging the mxn eigenvectors to obtain one eigenvector, i.e., the sample eigenvector.
In step S3, the rankfeatures algorithm uses a toolbox function of MATLAB software, and the algorithm sorts the key features according to the category separation criterion, specifically using the function
[IDX,Z]＝rankfeatures(X,Group)，
Sorting the features using a binary-sorted independence evaluation criterion, the input X being a matrix, each column representing a vector of observations and the number of rows corresponding to the original number of features; the input Group comprises class labels, the output IDX is a list formed by row subscripts of X, the subscript list is arranged from high to low according to the significance of the features in the matrix X, and the output Z is the weight corresponding to each feature after the use criterion. The two-sample t-test criterion is estimated by default using the joint variance.
Calculating a sample characteristic fusion index by using a rankfeatures algorithm, wherein the method comprises the following specific steps:
s3-1: assigning a feature data matrix formed by all sample feature vectors to X, assigning vectors formed by all sample labels to Group, inputting the vectors into an [ IDX, Z ] ═ rankfeatures (X, Group) function for processing, and obtaining a feature weight vector Z containing 56 components;
s3-2: m (m <56) features with larger weights are selected for weighted summation to obtain a sample feature fusion index.
Example 2
The present embodiment has the same structure as embodiment 1 except for the following features:
as shown in fig. 1, the present embodiment is a method for extracting features of ultrasonic tissue characterization based on hilbert-yellow transform, which is applied to detecting microstructure changes of tissues before and after chemotherapy of breast cancer tissues, and includes the following steps:
the method comprises the following steps: constructing an ultrasonic RF time series
The chest cavity subcutaneous tissue of a female BALB/C nude mouse is scanned by using Sonix TOUCH produced by Ultrasonix Canada and a broadband linear array ultrasonic probe with the center frequency of 6.6MHz, and a multi-frame ultrasonic echo RF signal is recorded. 27 female BALB/C nude mice aged 4 to 6 weeks were used as subjects in this experiment. Firstly culturing human breast cancer cell MCF-7, adding 10% fetal bovine serum and 100 unit/ml streptomycin, placing in an environment with temperature of 37 deg.C and CO2 concentration of 5%, and then aseptically adding about 4 x 108One MCF-7 cell was inoculated subcutaneously into the thoracic cavity of a nude mouse. From a medical perspective, these nude mouse models can be studied as undifferentiated samples. When the tumor grows to about 5mm in diameter, the breast cancer animal model is divided into 14 treatment groups and 13 control groups, respectively, the treatment group nude mice are treated with 3mg/kg of paclitaxel per day by chemotherapy, and the control group is given physiological saline only. FIG. 2 is a pathological diagram of a breast tumor region in a control group and a treatment group, and changes of tissue microstructure such as significant reduction of tumor cell density, increase of gaps among cells and the like in the treatment group can be observed;
s101: scanning the breast tumor tissue of the experimental nude mouse to obtain an ultrasonic echo RF signal;
s102: selecting a 100 th frame of ultrasonic echo RF signal to demodulate and display a B-type image, wherein FIG. 3 is a B-type image of a breast tumor area of a control group and a treatment group, and comparing the two B-type images, the difference of the images of the tumor area is not obvious, and the change of the microstructure of the tissue cannot be identified;
s103: selecting a 20 x 70 point size ROI on the B-mode map;
s104: taking the echo signals of the previous 256 frames of each point in the ROI to obtain an ultrasonic RF time sequence with the length of 256 of 20 multiplied by 70 points, wherein the figure 4 is an ultrasonic RF time sequence distribution diagram;
step two: calculating the characteristic vector of the experimental nude mouse sample by using HHT algorithm
According to a method for calculating sample feature vectors, performing feature extraction on the time sequence of each experimental nude mouse sample ROI to obtain sample feature vectors containing 56 components;
step three: calculating the fusion index of the characteristics of the nude mouse sample by using a rankfeatures algorithm
S301: respectively calculating sample characteristic vectors for 27 nude mice, assigning a transpose matrix of a characteristic data matrix of 27 multiplied by 56 to X, assigning a vector formed by 27 sample labels to Group, and calculating to obtain weight vectors of 56 characteristics by using a formula [ IDX, Z ] ═ rankfeatures (X, Group);
s302: the first 16 features with larger weights are selected for weighted summation, and through calculation, in this embodiment, the feature fusion indexes of 14 samples in the control group are respectively: 0.8334, 1.2672, 0.9598, 16.4901, 9.1725, 1.0033, 3.0090, 1.2001, 0.8029, 1.6267, 1.2665, 2.0368, 1.3963, 1.3334; the characteristic fusion indexes of the 13 samples in the treatment group are respectively as follows: 1.6160,1.5381,0.7148,0.7618,0.7112,1.0661,0.8636,1.5183,0.9405,1.0340,0.8058,0.8384,0.6947.
The mean value of the characteristic fusion index of the control group was 3.0284, and the variance was 19.6270; the mean of the characteristic fusion index for the treatment group was 1.0079 and the variance was 0.1119.
It can be seen that the mean value of the characteristic fusion index of the control group is greater than that of the characteristic fusion index of the treatment group, and the variance of the characteristic fusion index of the control group is greater than that of the characteristic fusion index of the treatment group, and the result shows that the characteristic fusion index of the treatment group sample subjected to chemotherapy is smaller than that of the control group sample and has smaller fluctuation. It is demonstrated that the feature fusion index can be used to distinguish tissues with altered microstructure and thus can be used as a feature for tissue characterization.
The above embodiments are preferred embodiments of the present invention, but the present invention is not limited to the above embodiments, and any other changes, modifications, substitutions, combinations, and simplifications which do not depart from the spirit and principle of the present invention should be construed as equivalents thereof, and all such changes, modifications, substitutions, combinations, and simplifications are intended to be included in the scope of the present invention.

Claims (8)

1. The method for extracting the characteristics of the ultrasonic tissue characterization based on the Hilbert-Huang transform is characterized by comprising the following steps of:
s1: acquiring ultrasonic echo RF signal multiframes of a tissue sample, selecting an ROI, and constructing an ultrasonic RF time sequence; the specific process for constructing the ultrasonic RF time sequence comprises the following steps:
s1-1: scanning tissues to obtain multiple frames of ultrasonic echo RF signals;
s1-2: demodulating a certain frame of ultrasonic echo RF signal and displaying a B-type image;
s1-3: selecting an ROI with the size of M multiplied by N on the B-type picture;
s1-4: taking previous L frames of RF signals of each point in the ROI to construct M multiplied by N ultrasonic RF time sequences with the length of L;
s2: performing feature extraction on the ultrasonic RF time sequence by using a HHT algorithm to obtain a sample feature vector based on the relation among time, frequency and energy; the specific process of calculating the sample feature vector comprises the following steps:
s2-1: taking a time sequence x (t) in a sample ROI, and carrying out ensemble empirical mode decomposition to obtain n IMF components;
s2-2: extracting time domain features based on IMF components:
the number of IMF zero-crossing points is 5: IMF1-ZCs, IMF2-ZCs, IMF3-ZCs, IMF4-ZCs, IMF 5-ZCs;
IMF variance 5: IMF1-Var, IMF2-Var, IMF3-Var, IMF4-Var, IMF 5-Var;
IMF variance contribution rate of 5 IMF1-VarR, IMF2-VarR, IMF3-VarR, IMF4-VarR, IMF 5-VarR;
s2-3: respectively performing Hilbert spectrum analysis on the n-order IMFs of the time series x (t);
specifically, step S2-3 includes:
step a: hilbert transform of an IMF component c (t)
Figure FDA0002287102550000011
In which P is Cauchy main value, and then obtaining the analytic signal of c (t)
Figure FDA0002287102550000012
A (t) e in polar coordinate formiθ(t)Wherein the amplitude function
Figure FDA0002287102550000013
Function of energy
Figure FDA0002287102550000014
Phase function
Figure FDA0002287102550000015
The instantaneous frequency of c (t) can be calculated from the phase function
Figure FDA0002287102550000016
Step b: respectively carrying out Hilbert transformation on n-order IMFs of x (t) according to the step a, constructing an analytic signal, expressing the analytic signal in a polar coordinate form, and finally reconstructing the x (t) to obtain
Figure FDA0002287102550000021
The residual component r (t) is omitted here, Re representing the real part;
step c: balance
Figure FDA0002287102550000022
The Hilbert amplitude spectrum is called as a Hilbert spectrum for short, and the Hilbert marginal spectrum is obtained by integrating the Hilbert spectrum with time
Figure FDA0002287102550000023
In the formula
Figure FDA0002287102550000024
Represents the k-th order IMF component ck(t) Hilbert margin spectrum;
s2-4: extracting time domain-energy characteristics based on the amplitude function and the energy function:
5 IMF mean intensities: IMF1-AvgA, IMF2-AvgA, IMF3-AvgA, IMF4-AvgA, IMF 5-AvgA;
5 IMF energies: IMF1-Egy, IMF2-Egy, IMF3-Egy, IMF4-Egy, IMF 5-Egy;
s2-5: extracting time domain-frequency features based on instantaneous frequency:
5 IMF highest frequencies: IMF1-MaxF, IMF2-MaxF, IMF3-MaxF, IMF4-MaxF, IMF 5-MaxF;
s2-6: extracting energy-frequency characteristics based on Hilbert marginal spectrum:
5 IMF mean center frequencies: IMF1-MCF, IMF2-MCF, IMF3-MCF, IMF4-MCF, IMF 5-MCF; 1 x (t) mean center frequency Orig-MCF;
s2-7: extracting frequency domain-energy characteristics based on Hilbert marginal spectrum:
hilbert marginal spectral entropy origin-EgyS of x (t);
the normalized Hilbert marginal spectrum of x (t) has low-band energy, medium-low band energy, medium-high band energy, high-band energy origin-MargS 1, origin-MargS 2, origin-MargS 3 and origin-MargS 4;
s2-8: extracting frequency domain-energy curve fitting characteristics based on Hilbert marginal spectrum:
normalized Hilbert marginal spectrum for x (t) straight-line fit slope, intercept: O-MLFSlope, O-MLFinterp;
slope and intercept of curve fitting of a first-order polynomial exponential function: O-MEFSlope, O-MEFInterp;
7 sixth order curve fitting coefficients: O-MSOFa0, O-MSOFa1, O-MSOFa2, O-MSOFa3, O-MSOFa4, O-MSOFa5, O-MSOFa 6;
s2-9: extracting time domain-frequency domain-energy characteristics based on the Hilbert spectrum:
statistical characteristics of Hilbert spectral gray statistical histogram of x (t) -mean TFImgmean, variance TFImgSD, skewness TFImgSkaew, kurtosis TFImgKurto;
s2-10: taking M × N time sequences x in sample ROI respectivelyi(t), wherein i is 0,1, 2., mxn-1, repeating steps S2-1 to S2-9 to obtain mxn eigenvectors corresponding to the mxn time series, and finally averaging the mxn eigenvectors to obtain one eigenvector, that is, a sample eigenvector;
s3: and (4) performing feature screening and feature fusion by using a rankfeatures algorithm, calculating a sample feature fusion index, and taking the feature fusion index as the feature of the characterization organization.
2. The feature extraction method according to claim 1, wherein n IMF components are obtained in step S2-1, n ≧ 5.
3. The feature extraction method according to claim 1, wherein the specific process of step S2-1 includes:
step a: adding a gaussian white noise sequence to X (t) to obtain X (t) ═ X (t) + ω (t);
step b: determining all local maximum values and local minimum values of X (t), and obtaining an upper envelope line b and a lower envelope line b by adopting cubic spline interpolationmax(t) and bmin(t) calculating the average value
Figure FDA0002287102550000031
Find h1(t)＝X(t)-m1(t)；
Step c: if h is1(t) satisfies the condition of eigenmode function IMF, then h1(t) is the first IMF component of X (t); otherwise, repeating step b and h1(t) as a new X (t), calculate m11(t), judging again h11(t)＝h1(t)-m11(t) whether IMF conditions are met, if not, cycling k times until h1k(t)＝h1(k-1)(t)-m1k(t) IMF conditions are satisfied; note c1(t)＝h1k(t) is the first order IMF component of X (t), in which case the remaining component r1(t)＝X(t)-c1(t)；
Step d: will r is1(t) repeating steps a to c as new X (t) to obtainTo the second component c of X (t) satisfying the IMF condition2(t) cyclically executing the steps a to c until rn(t) becoming a monotonous function, and finishing the cycle, thereby obtaining a set of n (n is more than or equal to 5) order IMF components of the signal X (t);
step e: the steps a to d are repeatedly executed for N times, and different Gaussian white noise sequences omega are added each timei(t), decomposing to obtain a set of N groups of IMF components, and averaging the components of corresponding orders in the N groups of IMF components
Figure FDA0002287102550000041
As the final IMF component, ckAnd (t) represents a k-order IMF component obtained by EEMD decomposition of x (t).
4. The feature extraction method according to claim 1, wherein step S2-2 includes:
step a: for IMF component c of length Li(t), the method for counting the number of zero-crossing points is as follows: when j is more than or equal to 1 and less than or equal to L, if c (j) × c (j-1)<0, then it is marked as a zero crossing, and traverses [1, L ]]Finding IMF component ciAll zero crossings of (t) are obtainedi(t) zero-crossing points, so that zero-crossing points IMF1-ZCs, IMF2-ZCs, IMF3-ZCs, IMF4-ZCs and IMF5-ZCs of the first 5-order IMF can be obtained;
step b: according to the formula
Figure FDA0002287102550000042
Calculating the variances IMF1-Var, IMF2-Var, IMF3-Var, IMF4-Var and IMF5-Var of the first 5-order IMF, and then calculating the variances according to the formula
Figure FDA0002287102550000043
Variance contribution ratios of the first 5 th order IMF, IMF1-VarR, IMF2-VarR, IMF3-VarR, IMF4-VarR, and IMF5-VarR, were calculated.
5. The feature extraction method of claim 1, wherein the steps S2-4 to S2-8 specifically include:
in step S2-4, with respect to timeAfter the n-th order IMF of the inter-sequence x (t) is processed by the step a of the step S2-3, n amplitude functions a with the length L can be obtainedi(t), energy function ei(t), (i ═ 0,1, 2.., n-1), according to the formula
Figure FDA0002287102550000044
Calculating the average intensity of the IMF of the first 5 th order, and then according to the formula
Figure FDA0002287102550000045
Calculating the energy of the IMF of the first 5 orders;
in step S2-5, after the n-th order IMF of the time series x (t) is processed in step A of step S2-3, n instantaneous frequency functions f with length L can be obtainedi(t), (i ═ 0,1, 2.., n-1), according to the formula
Figure FDA0002287102550000046
Calculating the highest frequency of the IMF of the first 5 orders;
in step S2-6, after the n-order IMF of the time series x (t) is processed by the steps a-c of the step S2-3, n Hilbert marginal spectrums h containing m frequency points can be obtainedi(f) And x (t) Hilbert marginal spectrum h (f) containing m frequency points according to a formula
Figure FDA0002287102550000051
And
Figure FDA0002287102550000052
wherein h isi(fj) Magnitude, f, of the j-th frequency point of the Hilbert marginal spectrum representing the i-th IMF componentjRepresenting the j-th frequency point, calculating the average center frequency of the IMF of the first 5 orders and the average center frequency of x (t);
the step S2-7 includes the steps of:
step a: according to Shannon entropy formula
Figure FDA0002287102550000053
Computing Hilbert marginal spectral entropy of x (t)
Figure FDA0002287102550000054
Wherein
Figure FDA0002287102550000055
Step b: normalizing the marginal spectrum h (f) of x (t) to obtain
Figure FDA0002287102550000056
Step c: h is to beo(f) Divided into 4 frequency bins: low frequency band
Figure FDA0002287102550000057
Middle and low frequency band
Figure FDA0002287102550000058
Middle and high frequency band
Figure FDA0002287102550000059
And high frequency band
Figure FDA00022871025500000510
According to the formula respectively
Figure FDA00022871025500000511
Figure FDA00022871025500000512
And
Figure FDA00022871025500000513
calculating the energy of a low frequency band, a medium-high frequency band and a high frequency band of the normalized Hilbert marginal spectrum of x (t);
step S2-8 is implemented by using a least square method to normalize the Hilbert marginal spectrum h of x (t)o(f) Performing straight line fitting to obtain O-MLFSlope and O-MLFInterp, performing first-order polynomial exponential function curve fitting to obtain O-MEFSlope and O-MEFInterp, and performing sixth-order polynomial exponential function curve fitting to obtain O-MEFSlope and O-MEFInterpThe curve fitting of the order yields O-MSOFa0, O-MSOFa1, O-MSOFa2, O-MSOFa3, O-MSOFa4, O-MSOFa5, O-MSOFa 6.
6. The feature extraction method according to claim 1, wherein the step S2-9 specifically includes:
step a: after the n-order IMF of the time sequence x (t) is processed in the steps a to c of the step S2-3, a Hilbert amplitude spectrum H (f, t) of the time sequence x (t) can be obtained, the time t is taken as a horizontal axis, the frequency f is taken as a vertical axis, and the gray scale is controlled by using the amplitude to obtain a Hilbert spectrum time-frequency gray scale image of the time sequence x (t);
step b: then, a Hilbert spectrum time-frequency gray level image histogram h of the time sequence x (t) is obtained by taking the gray level of the image as a horizontal axis and the occurrence frequency of each gray level on the image as a vertical axisiWherein i ═ 0,1,2iExpressing the occurrence frequency of the ith gray level pixel of the gray histogram, and l expressing the gray level number of the time-frequency gray image;
step c: respectively using formulas
Figure FDA0002287102550000061
Figure FDA0002287102550000062
And
Figure FDA0002287102550000063
and calculating the mean, variance, skewness and kurtosis of the Hilbert spectrum gray statistical histogram of x (t).
7. The feature extraction method of claim 1, wherein in step S3, the rankfeatures algorithm uses a toolbox function of MATLAB software, which sorts key features according to category separation criteria, in particular uses a function
[IDX,Z]＝rankfeatures(X,Group)，
Sorting the features using a binary-sorted independence evaluation criterion, the input X being a matrix, each column representing a vector of observations and the number of rows corresponding to the original number of features; the input Group comprises class labels, the output IDX is a list formed by row subscripts of X, the subscript list is arranged from high to low according to the significance of the features in the matrix X, and the output Z is the weight corresponding to each feature after the use criterion.
8. The feature extraction method according to claim 1, wherein the step S3 of calculating the sample feature fusion index by using a rankfeatures algorithm comprises the following specific steps:
s3-1: assigning a feature data matrix formed by all sample feature vectors to X, assigning vectors formed by all sample labels to Group, and inputting the vectors into an [ IDX, Z ] ═ rankfeatures (X, Group) function for processing to obtain a feature weight vector Z containing a plurality of components;
s3-2: and selecting the first m characteristics with large weight from the obtained characteristic weight vector Z containing a plurality of components to carry out weighted summation to obtain a sample characteristic fusion index.
CN201710417075.3A 2017-06-06 2017-06-06 Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform Expired - Fee Related CN107358156B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710417075.3A CN107358156B (en) 2017-06-06 2017-06-06 Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710417075.3A CN107358156B (en) 2017-06-06 2017-06-06 Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform

Publications (2)

Publication Number Publication Date
CN107358156A CN107358156A (en) 2017-11-17
CN107358156B true CN107358156B (en) 2020-05-19

Family

ID=60271035

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710417075.3A Expired - Fee Related CN107358156B (en) 2017-06-06 2017-06-06 Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform

Country Status (1)

Country Link
CN (1) CN107358156B (en)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108614926A (en) * 2018-04-12 2018-10-02 西安交通大学 A kind of modal parameters discrimination method being combined with Hilbert-Huang transform based on manifold learning
CN109214092A (en) * 2018-09-11 2019-01-15 吉林大学 The four component borehole strain data exception extraction methods based on Hilbert-Huang transform
CN110032585B (en) * 2019-04-02 2021-11-30 北京科技大学 Time sequence double-layer symbolization method and device
CN113205000A (en) * 2021-04-07 2021-08-03 江苏师范大学 Modulation mode identification method based on Cauchy Score polar coordinate diagram
CN113191232B (en) * 2021-04-21 2023-04-18 西安交通大学 Electro-hydrostatic actuator fault identification method based on multi-mode homologous features and XGboost model
CN114509604B (en) * 2022-04-18 2022-09-02 国网江西省电力有限公司电力科学研究院 GIS shell transient state ground potential rise waveform spectrum analysis method and system
CN116028882B (en) * 2023-03-29 2023-06-02 深圳市傲天科技股份有限公司 User labeling and classifying method, device, equipment and storage medium

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102855623A (en) * 2012-07-19 2013-01-02 哈尔滨工业大学 Method for measuring myocardium ultrasonic angiography image physiological parameters based on empirical mode decomposition (EMD)
CN105030279A (en) * 2015-06-24 2015-11-11 华南理工大学 Ultrasonic RF (radio frequency) time sequence-based tissue characterization method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102855623A (en) * 2012-07-19 2013-01-02 哈尔滨工业大学 Method for measuring myocardium ultrasonic angiography image physiological parameters based on empirical mode decomposition (EMD)
CN105030279A (en) * 2015-06-24 2015-11-11 华南理工大学 Ultrasonic RF (radio frequency) time sequence-based tissue characterization method

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Bouaguel Waad 等.A three-stage feature selection using quadratic programming for credit scoring.《Applied Artificial Intelligence》.2013,第27卷(第8期), *
Nishant Uniyal 等.Ultrasound RF Time Series for Classification of Breast Lesions.《IEEE TRANSACTIONS ON MEDICAL IMAGING》.2015,第34卷(第2期), *
R. Jurkonis 等.Algorithms and Results of Eye Tissues Differentiation Based on RF Ultrasound.《The Scientific World Journal》.2012,第2012卷(第2期), *

Also Published As

Publication number Publication date
CN107358156A (en) 2017-11-17

Similar Documents

Publication Publication Date Title
CN107358156B (en) Feature extraction method for ultrasonic tissue characterization based on Hilbert-Huang transform
EP3419516B1 (en) Ultrasound blood flow imaging
Zhou et al. High spatial–temporal resolution reconstruction of plane-wave ultrasound images with a multichannel multiscale convolutional neural network
JP5433097B2 (en) Ultrasonic observation apparatus, operation method of ultrasonic observation apparatus, and operation program of ultrasonic observation apparatus
CN107346541B (en) Tissue characterization method based on ultrasonic radio frequency time series wavelet analysis
EP2904975B1 (en) Ultrasound observation device, operation method for ultrasound observation device, and operation program for ultrasound observation device
Besson et al. A sparse reconstruction framework for Fourier-based plane-wave imaging
CN104411250B (en) The method of operating of ultrasound observation apparatus, ultrasound observation apparatus
KR101610874B1 (en) Module for Processing Ultrasonic Signal Based on Spatial Coherence and Method for Processing Ultrasonic Signal
JP5881905B2 (en) Ultrasonic observation apparatus, operation method of ultrasonic observation apparatus, and operation program of ultrasonic observation apparatus
CN102905624A (en) Ultrasound observation device, operation method of ultrasound observation device, and operation program of ultrasound observation device
CN103381096A (en) Blood perfusion separation detecting and imaging method for bone surface capillary
Pham et al. Joint blind deconvolution and robust principal component analysis for blood flow estimation in medical ultrasound imaging
CN111248858B (en) Photoacoustic tomography reconstruction method based on frequency domain wave number domain
JPWO2012063930A1 (en) Ultrasonic diagnostic apparatus, method for operating ultrasonic diagnostic apparatus, and operation program for ultrasonic diagnostic apparatus
CN110840484B (en) Ultrasonic imaging method and device for adaptively matching optimal sound velocity and ultrasonic equipment
US9386965B2 (en) Ultrasound machine providing composite image data
Chahuara et al. Regularized framework for simultaneous estimation of ultrasonic attenuation and backscatter coefficients
Sultan et al. Estimation of Cortical Bone Strength Using CNN-based Regression Model
Wang et al. An easily-achieved time-domain beamformer for ultrafast ultrasound imaging based on compressive sensing
Hansen et al. In-vivo validation of fast spectral velocity estimation techniques
CN106037799B (en) Elastic parameter imaging method based on ultrasonic RF backscatter signal time frequency analysis
Wang et al. Gaussian wavelet based dynamic filtering (GWDF) method for medical ultrasound systems
Qu et al. The effect of sound-speed-image resolution on phase aberration correction for ultrasound computed tomography
CN109886907B (en) Baseband harmonic B-ultrasonic fusion method based on weight matrix algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200519