CN109684970B - Window length determination method for moving principal component analysis of structural dynamic response - Google Patents

Window length determination method for moving principal component analysis of structural dynamic response Download PDF

Info

Publication number
CN109684970B
CN109684970B CN201811547587.2A CN201811547587A CN109684970B CN 109684970 B CN109684970 B CN 109684970B CN 201811547587 A CN201811547587 A CN 201811547587A CN 109684970 B CN109684970 B CN 109684970B
Authority
CN
China
Prior art keywords
principal component
matrix
window
response
component analysis
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.)
Active
Application number
CN201811547587.2A
Other languages
Chinese (zh)
Other versions
CN109684970A (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.)
Jinan University
Original Assignee
Jinan University
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 Jinan University filed Critical Jinan University
Priority to CN201811547587.2A priority Critical patent/CN109684970B/en
Publication of CN109684970A publication Critical patent/CN109684970A/en
Application granted granted Critical
Publication of CN109684970B publication Critical patent/CN109684970B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08Feature extraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/14Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Signal Processing (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a structure dynamic response shiftThe window length determining method for dynamic principal component analysis comprises the following steps: installing sensor array on the structure, measuring dynamic response signal wm(n) obtaining a response signal matrix B, carrying out principal component analysis on the matrix B to obtain an eigenvalue diagonal matrix Λ, and utilizing the eigenvalue to obtain the front k-order principal component cumulative contribution ratio CCRkDetermining a k value by accumulating the contribution rate; for the measured response wm(n) Fourier transforming to obtain a spectrogram of the measured response, thereby determining the fundamental frequency f of the structural vibration1(ii) a Calculating the number T of sampling points in the first-order vibration period of the structure according to the fundamental frequency of the structure1(ii) a From T1Obtaining the integral multiple length rT1In-window signal B ofr(ii) a To BrPerforming principal component analysis to obtain characteristic value diagonal matrix Λ in the windowrΛ is made ofrCalculating front k-order principal component cumulative contribution rate CCR in windowkThereby obtaining a front k-order principal component cumulative contribution rate convergence function CCRk(r); determining the optimal multiple r of the first-order vibration period of the structure by using the convergence spectrum of the first-k-order principal component accumulated contribution ratep(ii) a Through rpThe optimal window length l is calculated.

Description

Window length determination method for moving principal component analysis of structural dynamic response
Technical Field
The invention relates to the technical field of structural dynamics signal processing, in particular to a window length determination method for moving principal component analysis of structural dynamic response.
Background
Principal component analysis is to characterize many variables of the original data with fewer principal components and also to retain enough information of the original variables, thereby reducing the difficulty of data analysis. Principal component analysis has two main effects: firstly, data dimension can be reduced and analysis time is reduced by compressing data; the second is the interpretation of the data. By adopting the principal component analysis method, the data containing a plurality of variables can be summarized into a plurality of principal component data, redundant information is eliminated, the complex problem is simplified, the obtained data is more scientific and accurate, and the model can reflect the real situation to the maximum extent. Principal component analysis has been widely used in various industries such as economy, finance, medicine, electronic communications, the internet, and civil engineering as a typical and effective statistical method.
Principal component analysis is also widely used in the field of structure health monitoring based on dynamics, and the principal components and characteristic values are used to evaluate the health condition of the structure. However, if the global data obtained by monitoring is analyzed at one time, three disadvantages are caused: one is that a principal component matrix and a characteristic value matrix are obtained intelligently, and the health condition of the structure cannot be reflected; secondly, the requirement cannot be met on the time effect, and the real-time requirement of structural health monitoring cannot be met; and thirdly, the calculation amount is large. Therefore, in recent years, researchers have proposed a moving principal component analysis method in which measured data is captured in the form of a moving window and moved as the measurement time moves, thereby reducing the amount of calculation.
However, the selection of the length of the moving window is a key step of the moving principal component analysis method. Environmental and noise directly affect the principal component accuracy, so the size of the moving window must be large enough to minimize this effect during steady state. However, when the window length is too large, the amount of calculation increases rapidly, a delay in health diagnosis time occurs, and the health condition of the structure cannot be reflected in real time. However, in the prior art methods, the method of determining the window length is determined empirically, and there is no unified and effective theoretical support and calculation process. In order to achieve the accuracy of calculating the principal component, the existing method selects a window large enough, for example, the time window length selects one year as a period to monitor the condition of the structure, and the selection method is completely not preferable in the real-time monitoring of the structure, and the diagnosis time delay is too large.
In order to solve the above problems, a method for determining the length of the moving window by using a principal component cumulative contribution rate convergence spectrum theory is urgently needed to be proposed, so that the calculation accuracy is ensured, and meanwhile, the real-time performance of monitoring is ensured.
Disclosure of Invention
The present invention is directed to solve the above-mentioned drawbacks of the prior art, and provides a method for determining a window length for a moving principal component analysis of a structural dynamic response.
The purpose of the invention can be achieved by adopting the following technical scheme:
a method for determining the length of a moving window based on principal component cumulative contribution rate convergence spectrum theory comprises the following steps:
s1, measuring the dynamic response of the structure with the sensor array as wm(n), obtaining a response signal matrix B, wherein the response signal matrix B is as follows:
Figure BDA0001909770120000021
in the formula (1), N is 1,2, …, N, M is 1,2, …, M and N are the lengths of signal sampling points, and M is the number of the measuring points;
s2, performing Principal Component Analysis (PCA) on the response signal matrix B in the formula (1) to obtain
PCA(B)=[U,S,Λ](2)
Wherein U is principal component matrix, S is principal component score matrix corresponding to U, Λ is eigenvalue matrix corresponding to U, Λ is diagonal matrix of M × M, and eigenvalue λ isi1, M is arranged diagonally from large to small:
Figure BDA0001909770120000031
s3 passing the characteristic value lambdaiCalculating the front k-order principal component cumulative contribution rate CCRkThe number of principal components to be considered, i.e., the k value, is determined by the cumulative contribution ratio. The first k-th order principal component cumulative contribution rate is:
Figure BDA0001909770120000032
wherein CCRkThe first k principal component contribution rate. And determining the k value when the current k-order principal component cumulative contribution rate reaches 90%.
S4, Fourier transform (FFT) is carried out on the measured response signal, and the fundamental frequency f of the structure is obtained from the spectrogram1
S5 passing through fundamental frequency f1Calculate f1The number of sampling points of the first-order vibration period of the corresponding structure is as follows:
Figure BDA0001909770120000033
wherein f issIs the sampling frequency.
And S6, taking the length of the moving window as an integral multiple of the corresponding period of the fundamental frequency of the structure, and defining an independent variable r which is a multiple of the corresponding period of the fundamental frequency of the structure. Respectively have a length of rT1The measured response is intercepted, r is 1,2,3 …, and the response signal matrix B in the window is obtainedr
Figure BDA0001909770120000034
S7 response signal matrix B in the equation (6)rPerforming principal component analysis to obtain
PCA(Br)=[Ur,Srr](7)
Wherein U isrIs a principal component matrix, SrCorresponding UrΛrIs corresponding to UrΛ matrix of eigenvaluesrIs a diagonal matrix of M × M, eigenvalues
Figure BDA0001909770120000041
Arranged from large to small on the diagonal.
S8, calculating the front k-order principal component cumulative contribution rate convergence function CCRk(r) obtaining a convergence function for the principal component cumulative contribution ratio with an argument r as:
Figure BDA0001909770120000042
as r increases, the convergence function will converge to the top k principal component cumulative contribution ratio CCR determined in step S3k
S9, making a main component cumulative contribution rate convergence function curve, and determining the multiple of the optimal fundamental frequency corresponding to the period in the window. When principal component cumulative contribution rate converges to CCRk(r) when formula (9) is satisfied, rpIs the optimum value.
CCRk(rp)=CCRk(r)=CCRk,r≥rp(9)
Wherein r ispIs a multiple of the period corresponding to the optimum fundamental frequency. Satisfying the formula (9) means that when r ≧ rpThe convergence function converges to a stable value CCRk
S10, passing throughpCalculating an optimal window length of
Figure BDA0001909770120000043
Therefore, when the moving principal component analysis of the structure dynamic response is carried out, the window length is r of the corresponding period of the fundamental frequency of the structure vibrationpAnd (4) doubling.
Compared with the prior art, the invention has the following advantages and effects:
1) the invention provides a method for determining the length of a moving window by utilizing a principal component cumulative contribution rate convergence spectrum theory, which replaces the conventional method of determining the length of the moving window by taking an empirical value, ensures the calculation accuracy and ensures the real-time performance of monitoring.
2) The method provided by the invention has the advantages of simple operation and small calculation amount.
Drawings
FIG. 1 is a flow chart of a method for determining a window length of a dynamic response of a moving principal component analysis structure as disclosed in the present invention;
FIG. 2 is a schematic view of a beam bridge model according to an embodiment of the present invention;
FIG. 3 is a graph of acceleration response measured in an embodiment of the present invention;
FIG. 4 is a graph of cumulative contribution of principal component analysis in an embodiment of the present invention;
FIG. 5 is a graph of the frequency spectrum of a measured signal in an embodiment of the present invention;
FIG. 6 is a convergence spectrum of the first 3 rd order principal component cumulative contribution ratio in the example of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
Examples
As shown in fig. 1, fig. 1 is a flowchart of a method for determining a window length based on a moving principal component analysis of a structural dynamic response disclosed in the present invention. The structural model used in this example is a steel bridge model, and the schematic diagram is shown in fig. 2. The length l of the model beam is 20m, and the sampling frequency fsThe specific implementation process is as follows:
s1 measuring the dynamic response w of a structure with a sensor arraym(N), obtaining a response signal matrix B, where N is 1,2, …, N, M is 1,2, …, M is the length of signal sampling point, M is the number of measuring points, in this embodiment, N is 50000, M is 9, the measured dynamic response is as shown in fig. 3, and therefore the measured response signal matrix B is:
Figure BDA0001909770120000061
s2, performing Principal Component Analysis (PCA) on the response signal matrix B in the formula (1) to obtain
PCA(B)=[U,S,Λ](2)
Wherein U is a principal component matrix, S is a principal component score matrix corresponding to U, Λ is an eigenvalue matrix corresponding to U, Λ is a diagonal matrix of 9 × 9, and the eigenvalue λ isi1, 9 are arranged diagonally from large to small:
Figure BDA0001909770120000062
s3 passing the characteristic value lambdaiCalculating the front k-order principal component cumulative contribution rate CCRkDetermining the number of principal components to be considered by means of the cumulative contribution rate, i.e. determiningAnd (5) setting the k value. The first k-th order principal component cumulative contribution rate is:
Figure BDA0001909770120000063
wherein CCRkThe first k principal component contribution rate. When the current k-order principal component cumulative contribution rate reaches 90%, the k value can be determined. As shown in fig. 4, when k is 3, the cumulative contribution ratio CCR is 90.3%, which exceeds 90%, and thus k is determined to be 3 by the cumulative contribution ratio map.
S4, Fourier transform (FFT) is carried out on the measured response signal, and the fundamental frequency f of the structure is obtained from the spectrogram1. FIG. 5 is a Fourier spectrum of the measured dynamic response, from which it can be derived that the fundamental frequency of the beam bridge structure of the embodiment is 1.233 Hz.
S5 passing through fundamental frequency f1Calculate f1The number of sampling points in the first-order vibration period of the corresponding structure is as follows:
Figure BDA0001909770120000071
wherein f issFor sampling frequency, the number of sampling points in the first-order vibration period of the bridge structure is T1=162。
S6, taking the length of the moving window as the integral multiple of the corresponding period of the structure fundamental frequency, defining an independent variable r, wherein r is the multiple of the corresponding period of the structure fundamental frequency. Respectively have a length of rT1The measured response is intercepted, r is 1,2,3 …, and the response signal matrix B in the window is obtainedr
Figure BDA0001909770120000072
S7 response signal matrix B in the equation (6)rPerforming principal component analysis to obtain
PCA(Br)=[Ur,Srr](7)
Wherein U isrIs a principal component matrix, SrCorresponding UrΛrIs corresponding to UrThe eigenvalue matrix of 9 × 9 is a diagonal matrix, and the eigenvalues are arranged from large to small on the diagonal.
S8, calculating the first 3-order principal component cumulative contribution rate convergence function CCR3(r) obtaining a convergence function for the principal component cumulative contribution ratio with an argument r as:
Figure BDA0001909770120000073
as r increases, the convergence function converges to the first 3 < rd > principal component cumulative contribution ratio CCR3=90.3%。
S9, making a main component cumulative contribution rate convergence function curve, and determining the optimal multiple r of the fundamental frequency corresponding to the period in the windowp. When principal component cumulative contribution rate converges to CCR3=90.3%, that is, when the formula (9) is satisfied, rpIs the optimum value.
CCRk(rp)=CCRk(r)=CCRk,r≥rp(9)
Wherein r ispIs a multiple of the period corresponding to the optimum fundamental frequency. Satisfying the formula (9) means that when r ≧ rpThe convergence function converges to a stable value CCRk90.3%. FIG. 6 is a calculated convergence function curve of the first three principal component cumulative contribution ratios, which can be used to determine rp=58。
S10, passing throughpThe optimal window length is calculated as:
Figure BDA0001909770120000081
therefore, when the moving principal component analysis of the structure dynamic response is carried out, the window length is r of the corresponding period of the fundamental frequency of the structure vibrationpDoubling the window length l to 9396.
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. A window length determination method for moving principal component analysis of structural dynamic response is characterized by comprising the following steps:
s1, mounting a sensor array on the structure, and measuring the dynamic response signal wm(N), N is 1,2, …, N, M is 1,2, …, M, N is the length of signal sampling points, M is the number of measuring points, a response signal matrix B is obtained,
Figure FDA0002499508850000011
s2, performing principal component analysis on the response signal matrix B to obtain a characteristic value diagonal matrix Λ;
s3, calculating front k-order principal component cumulative contribution rate CCR by using eigenvaluekDetermining a k value by accumulating the contribution rate;
s4 response w to measuredm(n) Fourier transforming to obtain a spectrogram of the measured response, thereby determining the fundamental frequency f of the structural vibration1
S5, calculating the number T of sampling points in the first-order vibration period of the structure according to the fundamental frequency of the structure1
S6, from T1Obtaining the integral multiple length rT1In-window response signal matrix BrWherein r is a self-defined independent variable, r is a multiple of a corresponding period of the fundamental frequency of the structure, and r is 1,2,3 …;
s7 matrix B of response signalsrPerforming principal component analysis to obtain characteristic value diagonal matrix Λ in the windowr
S8, use ΛrCalculating front k-order principal component cumulative contribution rate CCR in windowk(r) to obtain a front k-th order principal component cumulative contribution rate convergence function CCRk(r);
S9, determining the optimal multiple r of the first-order vibration period of the structure by using the front k-order principal component accumulated contribution rate convergence spectrump
S10, passing throughpCalculating an optimal windowThe length l.
2. The method for determining the window length of the moving principal component analysis of structural dynamical response of claim 1, wherein the step S2 is performed as follows:
performing principal component analysis on the response signal matrix B to obtain
PCA(B)=[U,S,Λ](2)
Wherein U is principal component matrix, S is principal component score matrix corresponding to U, Λ is eigenvalue matrix corresponding to U, Λ is diagonal matrix of M × M, and eigenvalue λ isi1, M is arranged diagonally from large to small:
Figure FDA0002499508850000021
3. the method as claimed in claim 1, wherein the first k-th principal component cumulative contribution rate CCR is determined by a window length of a moving principal component analysis of a structural dynamical responsekComprises the following steps:
Figure FDA0002499508850000022
wherein λjAnd λiIs the eigenvalue.
4. The method for determining window length of moving Principal Component Analysis (PCA) of structural dynamical response (SRR) of claim 1, wherein the number of sampling points T in the first-order vibration period of the structure is calculated from the fundamental frequency of the structure in step S51The process is as follows:
Figure FDA0002499508850000023
wherein f issIs the sampling frequency.
5. A constructional movement as claimed in claim 1The method for determining the window length of the moving principal component analysis of force response is characterized in that in step S6, the length is rT1The measured response is intercepted, r is 1,2,3 …, and the response signal matrix B in the window is obtainedr
Figure FDA0002499508850000031
6. The method for determining the window length of the moving principal component analysis of structural dynamical response of claim 1, wherein the step S7 is performed as follows:
for response signal matrix BrPerforming principal component analysis to obtain
PCA(Br)=[Ur,Srr](7)
Wherein U isrIs a principal component matrix, SrCorresponding UrΛrIs corresponding to UrΛ matrix of eigenvaluesrIs a diagonal matrix of M × M, eigenvalues
Figure FDA0002499508850000032
Arranged from large to small on the diagonal.
7. The method of claim 1, wherein Λ is employed in step S8 to determine the window length of the moving principal component analysis of the structural dynamic responserCalculating front k-order principal component cumulative contribution rate convergence function CCR in windowk(r) obtaining a convergence function for the principal component cumulative contribution ratio with an argument r as:
Figure FDA0002499508850000033
wherein
Figure FDA0002499508850000034
And
Figure FDA0002499508850000035
for the eigenvalues, the convergence function CCR increases with rk(r) cumulative contribution ratio CCR of principal components converging to the first k-th orderk
8. The method for determining window length for moving principal component analysis of structural dynamical response of claim 1, wherein the step S10 is performed by rpThe optimal window length l is calculated as follows:
Figure FDA0002499508850000036
wherein f issIs the sampling frequency.
CN201811547587.2A 2018-12-18 2018-12-18 Window length determination method for moving principal component analysis of structural dynamic response Active CN109684970B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811547587.2A CN109684970B (en) 2018-12-18 2018-12-18 Window length determination method for moving principal component analysis of structural dynamic response

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811547587.2A CN109684970B (en) 2018-12-18 2018-12-18 Window length determination method for moving principal component analysis of structural dynamic response

Publications (2)

Publication Number Publication Date
CN109684970A CN109684970A (en) 2019-04-26
CN109684970B true CN109684970B (en) 2020-08-07

Family

ID=66186265

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811547587.2A Active CN109684970B (en) 2018-12-18 2018-12-18 Window length determination method for moving principal component analysis of structural dynamic response

Country Status (1)

Country Link
CN (1) CN109684970B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110657882B (en) * 2019-09-23 2021-07-27 暨南大学 Bridge real-time safety state monitoring method utilizing single-measuring-point response

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101523337B1 (en) * 2012-09-20 2015-05-28 주식회사 바이오센 Methods for Assessing Virus Infections against a Plant using Raman Spectroscopy and Chemometrics
CN104331595B (en) * 2014-09-04 2018-02-13 天津大学 The mobile principal component correlation coefficient process of bridge damnification early warning
CN105911153B (en) * 2016-04-08 2018-07-13 暨南大学 A kind of Signal separator and denoising method and device based on mobile window function
CN106091972B (en) * 2016-06-30 2018-09-21 河海大学 A kind of building change detecting method projecting dot density based on moving window
CN107918053B (en) * 2017-10-21 2020-03-27 广东科瑞德电气科技有限公司 Fast slip calculation method based on window movement
CN108573224B (en) * 2018-04-04 2021-07-23 暨南大学 Bridge structure damage positioning method for mobile reconstruction of principal components by using single sensor information

Also Published As

Publication number Publication date
CN109684970A (en) 2019-04-26

Similar Documents

Publication Publication Date Title
KR102323514B1 (en) Method for measuring the axial force of bolts
CN116186634B (en) Intelligent management system for construction data of building engineering
CN110455490B (en) Method and device for calculating supersonic velocity and hypersonic velocity wind tunnel flow field turbulence
CN110160524B (en) Sensor data acquisition method and device of inertial navigation system
CN113901379B (en) Real-time data dynamic online quick processing method for edge terminal
CN112834193B (en) Operation bridge vibration and health state abnormity early warning method based on three-dimensional graph
CN112307619B (en) Construction method of early warning model, equipment fault early warning method and device
JP2023075189A (en) Condition monitoring system, method, and program
CN108898050A (en) A kind of flexible material process equipment roll shaft performance index calculation method
CN109684970B (en) Window length determination method for moving principal component analysis of structural dynamic response
CN116861313B (en) Kalman filtering working condition identification method and system based on vibration energy trend
CN117150217A (en) Data compression processing method based on big data
CN117109487B (en) Automatic nondestructive measurement method for metal thickness
CN107368457B (en) The instantaneous Frequency Estimation method examined based on LoG operator and Grubbs
CN117387884A (en) Bridge deflection measurement method based on multi-sensor data fusion
CN107014482B (en) Online monitoring device and method for vibration state
CN113933055B (en) Rolling bearing raceway defect size quantification method, device and system
CN114993671A (en) Vibration fault diagnosis method and system based on Q factor wavelet transform
JP2015096831A (en) Information processing device, information processing method, and program
CN109241147B (en) Method for evaluating variability of statistical value
JP2023508075A (en) Integrating physical sensors within a data assimilation framework
CN112016235A (en) Impact load identification method and system of flexible antenna structure
CN117889945B (en) Highway bridge construction vibration testing method
CN111735528B (en) Method for accurately positioning disturbance point based on KNN voting
WO2017130417A1 (en) Biological sound analysis device, biological sound analysis method, computer program, and recording medium

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