A method and device for high-precision estimation of amplitude and phase of rapidly dynamic changing signals

Abstract

The application discloses a kind of amplitude and phase high-precision estimation method and device suitable for fast dynamic change signal, the method of the application includes the discrete signal sequence of continuous signal sampling, constructs Hankel matrix and decomposes left singular matrix U, diagonal matrix D, right singular matrix V;The estimated value of the number of frequency components p is estimated by D, combined with extracting the subspace of right singular matrix V and using least square method to obtain diagonal matrix Q;The damping coefficient of frequency component is calculated by Q and constructs system matrix A;The dynamic phasor coefficient vector S is calculated by A and discrete signal sequence using least square method;The amplitude and phase of continuous signal are calculated by S and damping coefficient.The application aims at solving the problem that the main lobe overlap of dense harmonic / interharmonic with small frequency interval in power system causes spectral leakage interference, extremely serious even can submerge part of main spectrum line, and then lead to the accuracy of spectral leakage affecting the amplitude and phase estimation of dynamic change signal.

Description

Technical Field

[0001] This invention relates to the field of dynamic signal testing and analysis technology, specifically to a method and apparatus for high-precision estimation of amplitude and phase of rapidly changing signals. Background Technology

[0002] In power systems, we often need to consider the amplitude and phase of voltage and current, as they constitute a complete description of the power signal. However, in practical power systems, signals typically need to be sampled, truncated, and discretized for analysis, which can lead to spectral leakage. Spectral leakage in power systems can manifest as interference between different frequency components, potentially affecting power quality analysis, fault detection, and system control. Ideally, the spectrum of a power waveform should be clear, with each frequency component corresponding to a distinct spectral line. However, since the sampling length is not always an integer multiple of the signal period, practical power signals may contain additional components in the spectrum. These additional components are caused by truncation and sampling, resulting in spectral leakage across the entire frequency domain. Spectral leakage in power signals can affect the spectral analysis and harmonic analysis of power systems, adversely impacting the operational stability of the power system and the normal operation of equipment. Therefore, high-precision estimation of the amplitude and phase of rapidly changing signals is a critical challenge that urgently needs to be addressed.

[0003] High-precision estimation of signal amplitude and phase is a crucial issue in many engineering and technological fields. Currently, commonly used methods include quasi-synchronous algorithms, adaptive sampling frequency adjustment methods, windowed interpolation algorithms, corrected sampling sequences, and spectral leakage cancellation algorithms. However, these methods struggle to achieve high-precision estimation of amplitude and phase in rapidly changing signals. The main reasons are as follows: 1) Existing methods for estimating harmonic and interharmonic parameters mostly require prior conditions such as the visibility of the main spectrum of the frequency component to be estimated and the known number or approximate frequency of the frequency component. 2) Related improved methods are prone to serious errors due to requirements for synchronized sampling, difficulty in selecting the order of the computational model, or the need to ignore some spectral leakage interference. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a method and apparatus for high-precision estimation of amplitude and phase of rapidly changing signals, addressing the aforementioned problems of the prior art. This invention aims to solve the problem of spectral leakage interference caused by the overlap of main lobes of dense harmonics / interharmonics with small frequency intervals in power systems, which is extremely serious and may even submerge part of the main spectral lines, thus affecting the accuracy of amplitude and phase estimation of dynamically changing signals.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0006] A high-precision estimation method for amplitude and phase of signals suitable for rapidly changing dynamic signals includes the following steps:

[0007] S101, sampling continuous signals at a frequency F s Sampling yields a discrete signal sequence y(n) of length J;

[0008] S102, construct the Hankel matrix Y using the discrete signal sequence y(n), and decompose the singular values ​​of the Hankel matrix Y into a left singular matrix U, a diagonal matrix D, and a right singular matrix V;

[0009] S103, estimate the number of frequency components p using the diagonal matrix D.

[0010] S104, based on the estimated value of the number of frequency components p Given a right singular matrix V, extract the subspace of the right singular matrix V, and use the least squares method to obtain the diagonal matrix Q based on the subspace of the right singular matrix V;

[0011] S105, calculate the damping coefficient β of the k-th frequency component from the imaginary and real parts of the diagonal elements of the diagonal matrix Q. k The system matrix A is constructed using the eigenvalues ​​of the diagonal matrix Q;

[0012] S106, the dynamic phasor coefficient vector S is calculated using the least squares method from the system matrix A and the discrete signal sequence y(n);

[0013] S107, consisting of the dynamic phasor coefficient vector S and the damping coefficient β k Calculate the amplitude X and phase φ of a continuous signal.

[0014] Optionally, step S103 includes:

[0015] S201, extract singular value pairs from the diagonal matrix D shown in the following equation:

[0016] D = diag{c 11 ,c 12 ,c 21 ,c 22 ,…,c i1 ,c i2 ,…,c m1 ,c m2},

[0017] In the above formula, c 11 ,c 12 For the first pair of singular values, c 21 ,c 22 For the second pair of singular values, c i1 ,c i2 For the i-th singular value, cm1 ,c m2 Let y(n) be the m-th singular value, and i = 1, ..., m, m = (J+1) / 4, where J is the length of the discrete signal sequence y(n);

[0018] S202, summing each pair of singular values ​​according to the following formula:

[0019] g i =c i1 +c i2 ,

[0020] In the above formula, g i The sum of the i-th pair of singular values;

[0021] S203, arrange the sums of m pairs of singular values ​​in descending order to form a descending sequence:

[0022] [g1,g2,…,g i ,…,g (J+1) / 4 ],

[0023] In the above formula, g1~g (J+1) / 4 These are the elements from 1 to m in the descending sequence;

[0024] S204, Calculate the relative ratio sequence based on the results after sorting in descending order:

[0025]

[0026] In the above formula, σ1~σ (J+1) / 4 Let be the logarithms of the sums of the first to m pairs of singular values, and we have:

[0027] σ i =log g i / log g1,

[0028] In the above formula, σ i Let g be the logarithm of the sum of the i-th pair of singular values. i Let g1 be the i-th element in the descending sequence, and g1 be the 1-th element in the descending sequence.

[0029] S205, take the largest index value from the relative ratio sequence as the estimate of the number of frequency components p.

[0030] Optionally, step S104 includes:

[0031] S301, based on the estimated value of the number of frequency components p Given a right singular matrix V, extract two subspaces Vi from V according to the following formula. 11 and V 12 :

[0032]

[0033]

[0034] In the above formula, v(1) is the first column of the right singular matrix V, and v(2) is the second column of the right singular matrix V. The right singular matrix V is the first List, The right singular matrix V is the first The column extracts two subspaces V of the right singular matrix V according to the following formula. 21 and V 22 :

[0035]

[0036]

[0037] In the above formula, v(2) is the second column of the right singular matrix V, and v(3) is the third column of the right singular matrix V. The right singular matrix V is the first List, The right singular matrix V is the first List;

[0038] S302, combined with subspace V 11 and V 12 Matrix Q1 is obtained using the least squares method based on the following formula:

[0039]

[0040] In the above formula, eig is the function for solving eigenvalues ​​and eigenvectors. For subspace V 11 The pseudo-inverse matrix, For subspace V 11 The conjugate transpose of . These are the eigenvalues ​​and eigenvectors corresponding to the first frequency component. The eigenvalues ​​and eigenvectors corresponding to the second frequency component are: Let i be the eigenvalue and eigenvector corresponding to the i-th frequency component. For the first The eigenvalues ​​and eigenvectors corresponding to each frequency component; combined with the subspace V 21 and V 22 Matrix Q2 is obtained using the least squares method based on the following formula:

[0041]

[0042] In the above formula, For subspace V 21 The pseudo-inverse matrix, For subspace V 21 The conjugate transpose of; These are the eigenvalues ​​and eigenvectors corresponding to the first frequency component. The eigenvalues ​​and eigenvectors corresponding to the second frequency component are: Let i be the eigenvalue and eigenvector corresponding to the i-th frequency component. For the first Eigenvalues ​​and eigenvectors corresponding to each frequency component;

[0043] S303, combine matrix Q1 and matrix Q2 to form a diagonal matrix Q, and the eigenvalues ​​corresponding to each frequency component in the diagonal matrix Q are the average of the corresponding eigenvalues ​​in matrix Q1 and matrix Q2.

[0044] Optionally, in step S105, the damping coefficient β of the k-th frequency component is calculated. k The function expression is:

[0045]

[0046] In the above formula, β k1 and β k2 Let be the damping coefficients of the k-th frequency components in matrices Q1 and Q2, respectively, and we have:

[0047]

[0048]

[0049] In the above formula, Re and Im represent the real and imaginary parts of the complex number, respectively, Q1(k,k) is the element in the kth row and kth column of matrix Q1, and Q2(k,k) is the element in the kth row and kth column of matrix Q2.

[0050] Optionally, in step S105, the function expression for constructing the system matrix A using the eigenvalues ​​of the diagonal matrix Q is:

[0051] A = [A1 A2 A3],

[0052]

[0053]

[0054]

[0055] In the above formula, J is the length of the discrete signal sequence y(n), λ1+jω1 represents the first eigenvalue of the diagonal matrix Q, and λ2+jω2 represents the second eigenvalue of the diagonal matrix Q. Describe the first diagonal matrix Q. There are eigenvalues, where λ1 is the attenuation factor of the fundamental phasor, j is the imaginary unit, and ω1 is the angular frequency of the fundamental phasor; T s For the sampling period, and F s =1 / T s λ2 is the attenuation factor corresponding to the second frequency component, and ω2 is the angular frequency corresponding to the second frequency component; For the first The attenuation factor corresponding to each frequency component For the first The angular frequencies corresponding to each frequency component, and λ k =β k T s ω k =2πf k T s , λ k Let ω be the attenuation factor for the k-th frequency component. k Let T be the angular frequency of the k-th frequency component. s The sampling period is related to the sampling frequency F. s Satisfy F s =1 / T s f k This is the kth frequency component.

[0056] Optionally, the function expression for calculating the dynamic phasor coefficient vector S from the system matrix A and the discrete signal sequence y(n) using the least squares method in step S106 is as follows:

[0057] S=(A H A) -1 A H y,

[0058] In the above formula, A H Let y be the conjugate transpose of the system matrix A, and y be the discrete signal sequence y(n); and the resulting computational dynamic phasor coefficient vector S has the following form:

[0059]

[0060] In the above formula, Let be the complex constant obtained by least squares fitting, and T denote the transpose.

[0061] Optionally, in step S107, the dynamic phasor coefficient vector S and the damping coefficient β are used. k The functional expressions for calculating the amplitude X and phase φ of a continuous signal are as follows:

[0062]

[0063]

[0064]

[0065]

[0066]

[0067]

[0068] In the above formula, X (0) φ is the amplitude of the fundamental phasor of the dynamic signal. (0) X represents the phase of the fundamental phasor of the dynamic signal. (1) φ is the rate of change of the fundamental phasor amplitude of the dynamic signal. (1) X is the frequency of the fundamental phasor of the dynamic signal. (2) φ is the rate of change of the amplitude of the fundamental phasor of the dynamic signal. (2) Let β1 be the frequency change rate of the fundamental phasor of the dynamic signal, Re and Im represent the real and imaginary parts of the complex number, respectively, and β1 be the damping coefficient of the fundamental phasor. and is a complex constant obtained by least squares fitting in the dynamic phasor coefficient vector S.

[0069] Optionally, the functional expression of the Hankel matrix Y constructed in step S102 is:

[0070]

[0071] In the above formula, y(1) is the first discrete signal, y(2) is the second discrete signal, y((J+1) / 2) is the (J+1) / 2th discrete signal, y((J+3) / 2) is the (J+3) / 2th discrete signal, y(J) is the Jth discrete signal, and J is the length of the discrete signal sequence y(n). The function expression for decomposing the singular values ​​of the Hankel matrix Y into the left singular matrix U, the diagonal matrix D, and the right singular matrix V in step S102 is as follows:

[0072] UDV = svd(Y),

[0073] In the above formula, U is the left singular matrix, D is the diagonal matrix, V is the right singular matrix, svd is the singular value decomposition, Y is the Hankel matrix, and the size of the left singular matrix U, the diagonal matrix D, and the right singular matrix V are all ((J+1) / 2)×((J+1) / 2), where J is the length of the discrete signal sequence y(n).

[0074] Furthermore, the present invention also provides a high-precision amplitude and phase estimation device suitable for rapidly changing signals, comprising a microprocessor and a memory interconnected thereto, wherein the microprocessor is programmed or configured to execute the high-precision amplitude and phase estimation method suitable for rapidly changing signals.

[0075] Furthermore, the present invention also provides a computer-readable storage medium storing a computer program that is programmed or configured by a microprocessor to execute the high-precision amplitude and phase estimation method suitable for rapidly changing signals.

[0076] Compared with existing technologies, this invention has the following advantages: Addressing the problem of spectral leakage interference caused by the overlap of main lobes of dense harmonics / interharmonics with small frequency intervals in power systems, which can be extremely severe and even submerge part of the main spectral lines, thus affecting the accuracy of amplitude and phase estimation of dynamically changing signals, this invention's method includes sampling a continuous signal to obtain a discrete signal sequence, constructing a Hankel matrix and decomposing it into a left singular matrix U, a diagonal matrix D, and a right singular matrix V; and estimating the number of frequency components p through D. Combination The subspace of the right singular matrix V is extracted and the diagonal matrix Q is obtained using the least squares method. The damping coefficients of the frequency components are calculated from Q, and the system matrix A is constructed. The dynamic phasor coefficient vector S is calculated from A and the discrete signal sequence using the least squares method. The amplitude and phase of the continuous signal are calculated from S and the damping coefficients. That is, by selecting an appropriate rectangular window length, singular value decomposition is performed on the sampled discrete signal. The relative ratio is calculated based on the singular value decomposition results to determine the model order. The system matrix is ​​established using the least squares method, and the amplitude and phase are estimated with high precision to suppress spectral leakage, thereby improving the accuracy of amplitude and phase estimation of dynamically changing signals. Attached Figure Description

[0077] Figure 1 This is a schematic diagram of the basic process of the method in an embodiment of the present invention. Detailed Implementation

[0078] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0079] like Figure 1 The high-precision estimation method for amplitude and phase of rapidly changing signals in this embodiment includes the following steps:

[0080] S101, sampling continuous signals at a frequency F s Sampling yields a discrete signal sequence y(n) of length J;

[0081] S102, construct the Hankel matrix Y using the discrete signal sequence y(n), and decompose the singular values ​​of the Hankel matrix Y into a left singular matrix U, a diagonal matrix D, and a right singular matrix V;

[0082] S103, estimate the number of frequency components p using the diagonal matrix D.

[0083] S104, based on the estimated value of the number of frequency components p Given a right singular matrix V, extract the subspace of the right singular matrix V, and use the least squares method to obtain the diagonal matrix Q based on the subspace of the right singular matrix V;

[0084] S105, calculate the damping coefficient β of the k-th frequency component from the imaginary and real parts of the diagonal elements of the diagonal matrix Q. k The system matrix A is constructed using the eigenvalues ​​of the diagonal matrix Q;

[0085] S106, the dynamic phasor coefficient vector S is calculated using the least squares method from the system matrix A and the discrete signal sequence y(n);

[0086] S107, consisting of the dynamic phasor coefficient vector S and the damping coefficient β k Calculate the amplitude X and phase φ of a continuous signal.

[0087] A continuous signal can be represented as consisting of an effective signal and noise, and its functional expression is:

[0088] y(n) = s(n) + z(n),

[0089]

[0090] In the above formula, z(n) is unknown Gaussian white noise, and s(n) is the effective signal. For an unknown parameter with a normalized damping coefficient, λ k =β k T s and ω k =2πf k T s All parameters are unknown (k = 2, 3, ..., p), where p is the number of frequency components. To verify that this embodiment is suitable for high-precision amplitude and phase estimation of rapidly changing signals, this embodiment uses a sample of a continuous signal output from a Fluke 6100A three-phase standard power source. Its signal model expression is:

[0091]

[0092] In the above formula, u(t) is the voltage of the continuous signal at time t, X(t) is the amplitude of the continuous signal at time t, and f b Let X be the fundamental frequency of the power system (50Hz or 60Hz), φ(t) be the phase of the continuous signal at time t, and X be the frequency of the fundamental frequency of the power system (50Hz or 60Hz). i f is the amplitude of the i-th frequency component (i-th harmonic). i φ represents the frequency magnitude of the i-th frequency component. iThe i-th frequency component phase, d is the number of harmonics. Step 101 in this embodiment includes adjusting the continuous signal using F... s =1 / T s The sampling frequency is used to sample the data. The sequence y(n) = {y(1), y(2), ..., y(J)} of length J is taken as the discrete signal sequence y(n) of length J to be analyzed, where J = 8b-1 (b is a positive integer). In this embodiment, the sequence y(n) = {y(1), y(2), ..., y(199)} of length J is taken as the discrete signal sequence y(n) to be analyzed.

[0093] The functional expression of the Hankel matrix Y constructed in step S102 of this embodiment is:

[0094]

[0095] In the above formula, y(1) is the first discrete signal, y(2) is the second discrete signal, y((J+1) / 2) is the (J+1) / 2th discrete signal, y((J+3) / 2) is the (J+3) / 2th discrete signal, y(J) is the Jth discrete signal, and J is the length of the discrete signal sequence y(n). When J = 199, then:

[0096]

[0097] In step S102, the functional expression for decomposing the singular values ​​of the Hankel matrix Y into a left singular matrix U, a diagonal matrix D, and a right singular matrix V is as follows:

[0098] UDV = svd(Y),

[0099] In the above formula, U is the left singular matrix, D is the diagonal matrix, V is the right singular matrix, svd is the singular value decomposition, Y is the Hankel matrix, and the size of the left singular matrix U, the diagonal matrix D, and the right singular matrix V are all ((J+1) / 2)×((J+1) / 2), where J is the length of the discrete signal sequence y(n). In this embodiment, the specific sizes of the left singular matrix U, the diagonal matrix D, and the right singular matrix V are all 100×100. The singular value decomposition svd is a well-known method, so its implementation details will not be described here.

[0100] In step S103, the diagonal matrix D in the singular value decomposition result is used to estimate the number p of the signal frequency components. In traditional methods, the estimation results may be unstable due to the potential for noise interference affecting threshold selection. This patent uses a logarithmic method to calculate the relative ratio to determine the number of signal frequency components. It adaptively estimates the number of frequency components based on the relative magnitude of singular values, eliminating biases caused by numerical variations when using only absolute values. This method offers higher robustness against noise and interference and stronger adaptability to different signals. The elements of the diagonal matrix D are:

[0101] D = diag{c 11 ,c 12 ,c 21 ,c 22 ,…,c i1 ,c i2 ,…,c m1 ,c m2},

[0102] In the above formula, c i1 c i2 There are a pair of singular values ​​(i = 1, ..., m, and m = (J+1) / 4), and a total of (J+1) / 2 singular values. For example, when J = 199, m = 50, resulting in a total of 100 singular values. Therefore, step S103 in this embodiment includes:

[0103] S201, extract singular value pairs from the diagonal matrix D shown in the following equation:

[0104] D = diag{c 11 ,c 12 ,c 21 ,c 22 ,…,c i1 ,c i2 ,…,c m1 ,c m2},

[0105] In the above formula, c 11 ,c 12 For the first pair of singular values, c 21 ,c 22 For the second pair of singular values, c i1 ,c i2 For the i-th singular value, c m1 ,c m2 Let y(n) be the m-th singular value, and i = 1, ..., m, m = (J+1) / 4, where J is the length of the discrete signal sequence y(n);

[0106] S202, summing each pair of singular values ​​according to the following formula:

[0107] g i =c i1 +c i2 ,

[0108] In the above formula, g i The sum of the i-th pair of singular values;

[0109] S203, arrange the sums of m pairs of singular values ​​in descending order to form a descending sequence:

[0110] [g1,g2,…,g i ,…,g (J+1) / 4 ],

[0111] In the above formula, g1~g (J+1) / 4 These are the elements from 1 to m in the descending sequence;

[0112] S204, Calculate the relative ratio sequence based on the results after sorting in descending order:

[0113]

[0114] In the above formula, σ1~σ (J+1) / 4 Let be the logarithms of the sums of the first to m pairs of singular values, and we have:

[0115] σ i =log g i / log g1,

[0116] In the above formula, σ i Let g be the logarithm of the sum of the i-th pair of singular values. i Let g1 be the i-th element in the descending sequence, and g1 be the 1-th element in the descending sequence.

[0117] S205, take the largest index value from the relative ratio sequence as the estimate of the number of frequency components p. As an optional implementation, if multiple identical maximum values ​​appear, this embodiment preferentially selects the one with the smallest index value among the identical values ​​as the number of signal frequency components to ensure model stability. In this embodiment, the index value of the maximum value in the relative ratio sequence is 4, so the estimated value of the number of frequency components p is...

[0118] It should be noted that when extracting the subspace of the right singular matrix V in step S104 and obtaining the diagonal matrix Q using the least squares method based on the subspace of the right singular matrix V, the number of subspaces extracted each time and the number of times the subspaces are extracted can be set as needed. Each extracted subspace can be obtained by using the least squares method to obtain a submatrix, and then the eigenvalues ​​in the submatrix are summed to obtain the total matrix, i.e., the diagonal matrix Q.

[0119] For example, as an optional implementation, in this embodiment, the number of subspaces extracted each time is 2, and the number of subspace extractions is 2. Each extracted subspace can be obtained by using the least squares method to obtain a submatrix, namely: matrix Q1 and matrix Q2. Matrix Q1 and matrix Q2 are combined to form a diagonal matrix Q, and the eigenvalues ​​corresponding to each frequency component in the diagonal matrix Q are the average of the corresponding eigenvalues ​​in matrix Q1 and matrix Q2. Specifically, step S104 in this embodiment includes:

[0120] S301, based on the estimated value of the number of frequency components p Given a right singular matrix V, extract two subspaces Vi from V according to the following formula. 11 and V 12 :

[0121]

[0122]

[0123] In the above formula, v(1) is the first column of the right singular matrix V, and v(2) is the second column of the right singular matrix V. The right singular matrix V is the first List, The right singular matrix V is the first To improve the robustness of coarse frequency estimation performance under colored noise interference, the estimated value is based on the number of signal frequency components. Given a right singular matrix V, extract the two subspaces V of the right singular matrix V by two shifts. 11 and V 12 Since the first and last lines typically contain more noise information against a background of colored noise, V will be... 11 and V 12 The last and first rows are deleted as shown in the above formula; similarly, following the above method, the two subspaces of the right singular matrix V are extracted again by shifting, and the last and first rows are deleted respectively, resulting in two subspaces V. 21 and V 22 Specifically, the two subspaces V of the right singular matrix V are extracted according to the following formula. 21 and V 22 :

[0124]

[0125]

[0126] In the above formula, v(2) is the second column of the right singular matrix V, and v(3) is the third column of the right singular matrix V. The right singular matrix V is the first List, The right singular matrix V is the first List;

[0127] The estimated value of the number p of frequency components in this embodiment Therefore, we have:

[0128] V 11 =[v(1),v(2)],

[0129] V 12 =[v(2),v(3)],

[0130] V 21 =[v(2),v(3)],

[0131] V 22 =[v(3),v(4)],

[0132] S302, combined with subspace V 11 and V 12 Matrix Q1 is obtained using the least squares method based on the following formula:

[0133]

[0134] In the above formula, eig is the function for solving eigenvalues ​​and eigenvectors. For subspace V 11 The pseudo-inverse matrix, For subspace V 11 The conjugate transpose of . These are the eigenvalues ​​and eigenvectors corresponding to the first frequency component. The eigenvalues ​​and eigenvectors corresponding to the second frequency component are: Let i be the eigenvalue and eigenvector corresponding to the i-th frequency component. For the first The eigenvalues ​​and eigenvectors corresponding to each frequency component; combined with the subspace V 21 and V 22 Matrix Q2 is obtained using the least squares method based on the following formula:

[0135]

[0136] In the above formula, For subspace V 21 The pseudo-inverse matrix, For subspace V 21 The conjugate transpose of; These are the eigenvalues ​​and eigenvectors corresponding to the first frequency component. The eigenvalues ​​and eigenvectors corresponding to the second frequency component are: Let i be the eigenvalue and eigenvector corresponding to the i-th frequency component. For the first Eigenvalues ​​and eigenvectors corresponding to each frequency component;

[0137] S303, combine matrix Q1 and matrix Q2 to form a diagonal matrix Q, and the eigenvalues ​​corresponding to each frequency component in the diagonal matrix Q are the average of the corresponding eigenvalues ​​in matrix Q1 and matrix Q2.

[0138] The frequency f of the kth harmonic can be calculated from the imaginary and real parts of the diagonal elements of diagonal matrices Q1 and Q2. k (Hz) and damping coefficient β k ( / s), The corresponding results are then averaged to obtain a more accurate value. Specifically, in step S105, the damping coefficient β of the k-th frequency component is calculated. k The function expression is:

[0139]

[0140] In the above formula, β k1 and β k2 Let be the damping coefficients of the k-th frequency components in matrices Q1 and Q2, respectively, and we have:

[0141]

[0142]

[0143] In the above formula, Re and Im represent the real and imaginary parts of the complex number, respectively; Q1(k,k) is the element in the k-th row and k-th column of matrix Q1; and Q2(k,k) is the element in the k-th row and k-th column of matrix Q2. The expression for the calculation function of the k-th frequency component is:

[0144]

[0145]

[0146]

[0147] In the above formula, f k1 and f k2 These are the frequency magnitudes of the k-th frequency components in matrices Q1 and Q2, respectively.

[0148] Specifically, in this embodiment:

[0149] f 21 =59.5Hz, f 22 =53.4Hz,

[0150] β 21 =0.11s, β 22 =0.10s,

[0151] f 31 =99.6Hz, f 32 =99.9Hz,

[0152] β 31 =0.36s, β 32 =0.35s,

[0153] f 41 =150.1Hz, f 42 =150.1Hz,

[0154] β 41 =0.15s, β 42 =0.16s,

[0155] In step S105, the system matrix A is constructed using the eigenvalues ​​of the diagonal matrix Q; the system matrix A is a second-order Taylor series extended Vandermonde matrix of order J×(p+4), and λ is calculated. k =β k T s and ω k =2πf k T s The estimated value, The system matrix A is obtained. Specifically, the functional expression for constructing the system matrix A using the eigenvalues ​​of the diagonal matrix Q in step S105 is as follows:

[0156] A = [A1 A2 A3],

[0157]

[0158]

[0159]

[0160] In the above formula, J is the length of the discrete signal sequence y(n), λ1+jω1 represents the first eigenvalue of the diagonal matrix Q, and λ2+jω2 represents the second eigenvalue of the diagonal matrix Q. Describe the first diagonal matrix Q. There are eigenvalues, where λ1 is the attenuation factor of the fundamental phasor, j is the imaginary unit, and ω1 is the angular frequency of the fundamental phasor; T s For the sampling period, and F s =1 / T s λ2 is the attenuation factor corresponding to the second frequency component, and ω2 is the angular frequency corresponding to the second frequency component; For the first The attenuation factor corresponding to each frequency component For the first The angular frequencies corresponding to each frequency component, and λ k =β k T s ω k =2πf k T s , λ k Let ω be the attenuation factor for the k-th frequency component. k Let T be the angular frequency of the k-th frequency component. s The sampling period is related to the sampling frequency F. s Satisfy F s =1 / T s f k Let be the frequency magnitude of the k-th frequency component.

[0161] In step S106 of this embodiment, the function expression for calculating the dynamic phasor coefficient vector S using the least squares method from the system matrix A and the discrete signal sequence y(n) is as follows:

[0162] S=(A H A) -1 A H y,

[0163] In the above formula, A H Let y be the conjugate transpose of the system matrix A, and y be the discrete signal sequence y(n); and the resulting computational dynamic phasor coefficient vector S has the following form:

[0164]

[0165] In the above formula, Let be the complex constant obtained by least squares fitting, and T denote the transpose.

[0166] In step S107 of this embodiment, the dynamic phasor coefficient vector S and the damping coefficient β are used. k The functional expressions for calculating the amplitude X and phase φ of a continuous signal are as follows:

[0167]

[0168]

[0169]

[0170]

[0171]

[0172]

[0173] In the above formula, X (0) φ is the amplitude of the fundamental phasor of the dynamic signal. (0) X represents the phase of the fundamental phasor of the dynamic signal. (1) φ is the rate of change of the fundamental phasor amplitude of the dynamic signal. (1) X is the frequency of the fundamental phasor of the dynamic signal. (2) φ is the rate of change of the amplitude of the fundamental phasor of the dynamic signal. (2) Let β1 be the frequency change rate of the fundamental phasor of the dynamic signal, Re and Im represent the real and imaginary parts of the complex number, respectively, and β1 be the damping coefficient of the fundamental phasor. and is a complex constant obtained by least squares fitting in the dynamic phasor coefficient vector S.

[0174] In this embodiment, three evaluation parameters—total vector error (TVE), frequency error (FE), and permissible rate of change of frequency error (RFEs)—are used to evaluate the phasor, frequency, and rate of change of frequency estimation errors of the method in this embodiment. The results are compared with those of other existing methods for high-precision estimation and extraction of amplitude and phase of rapidly changing signals. Other existing methods include DFT-DPM (see: D. Petri, D. Fontanelli and D. Macii, "A Frequency-Domain Algorithm for Dynamic Synchrophasor and Frequency Estimation," in IEEE Transactions on Instrumentation and Measurement, vol. 63, no. 10, pp. 2330-2340, Oct. 2014, doi: 10.1109 / TIM.2014.2308996.) and FAST-PMU (see: P. Castello, J. Liu, C. Muskas, P.A. Pegoraro, F. Ponci and A. Monti, "AFast and Accurate PMU Algorithm for P+M)." Class Measurement of Synchrophasor and Frequency," in IEEETransactions on Instrumentation and Measurement, vol.63, no.12, pp.2837-2845, Dec.2014, doi:10.1109 / TIM.2014.2323137.), Hybrid P / M-class (see literature: AJRoscoe, "Exploring the Relative Performance of Frequency-Tracking andFixed-Filter Phasor Measurement Unit Algorithms Under C37.118 TestProcedures, the Effects of Interharmonics, and Initial Attempts at Merging P-Class Response With M-Class Filtering," in IEEE Transactions on Instrumentation and Measurement, vol.62, no.8, pp.(2140-2153, Aug.2013, doi:10.1109 / TIM.2013.2265431.), and the final results are shown in Table 1.

[0175] Table 1: Comparison of amplitude and phase estimation results between the method in this embodiment and existing methods.

[0176]

[0177]

[0178] As can be seen from Table 1, compared with other existing methods, the amplitude and phase estimation results obtained based on the method of this embodiment effectively suppress spectral leakage and obtain more accurate amplitude and phase estimation results.

[0179] Furthermore, this embodiment also provides a high-precision amplitude and phase estimation device suitable for rapidly changing signals, including a microprocessor and a memory interconnected thereto, wherein the microprocessor is programmed or configured to execute the high-precision amplitude and phase estimation method suitable for rapidly changing signals.

[0180] Furthermore, this embodiment also provides a computer-readable storage medium storing a computer program that is programmed or configured by a microprocessor to execute the high-precision amplitude and phase estimation method suitable for rapidly changing signals.

[0181] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-readable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create a machine for implementing the process. Figure 1 One or more processes and / or boxes Figure 1The computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to operate in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The functions specified in one or more boxes. These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable apparatus for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0182] The above description is merely a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should also be considered within the scope of protection of the present invention.

Claims

1. A high-precision estimation method for amplitude and phase of rapidly changing signals, characterized in that... Includes the following steps: S101, sampling frequency for continuous signals Sampling yields a length of J The discrete signal sequence y(n); S102, Construct the Hankel matrix using the discrete signal sequence y(n). Y , the Hankel matrix Y Singular value decomposition consists of a left singular matrix U, a diagonal matrix D, and a right singular matrix. ; S103, Estimating the number of frequency components using the diagonal matrix D. p The estimated value ; S104, based on the number of frequency components p The estimated value And right singular matrix Extracting the right singular matrix The subspace of, and according to the right singular matrix The diagonal matrix is ​​obtained by using the least squares method in the subspace. The extraction of the right singular matrix When extracting subspaces, the number of subspaces extracted each time is 2, and the number of extractions is 2. Each extracted subspace can be obtained as a submatrix using the least squares method: matrix Q1 and matrix Q2. Matrix Q1 and matrix Q2 are then combined to form a diagonal matrix. And diagonal matrix The eigenvalues ​​corresponding to each frequency component are the average of the corresponding eigenvalues ​​in matrices Q1 and Q2. S105, composed of a diagonal matrix The damping coefficient of the k-th frequency component is calculated using the imaginary and real parts of the diagonal elements. Using a diagonal matrix Eigenvalues ​​construct system matrix A Damping coefficient of the k-th frequency component The expression for the computation function is: ， In the above formula, and Let be the damping coefficients of the k-th frequency components in matrices Q1 and Q2, respectively, and we have: ， ， In the above formula, and These represent the real and imaginary parts of a complex number, respectively. Let Q1 be the element in the kth row and kth column. Let Q2 be the element in the kth row and kth column. The sampling period; S106, from the system matrix A The dynamic phasor coefficient vector S is calculated using the least squares method with respect to the discrete signal sequence y(n); S107, consisting of the dynamic phasor coefficient vector S and the damping coefficient Calculate the amplitude X and phase of a continuous signal. .

2. The high-precision amplitude and phase estimation method suitable for rapidly dynamically changing signals according to claim 1, characterized in that, Step S103 includes: S201, extract singular value pairs from the diagonal matrix D shown in the following equation: ， In the above formula, For the first pair of singular values, For the second pair of singular values, For the i-th pair of singular values, Let m be the m-th singular value, and have i =1,…,m, , Let y(n) be the length of the discrete signal sequence; S202, summing each pair of singular values ​​according to the following formula: ， In the above formula, The sum of the i-th pair of singular values; S203, arrange the sums of m pairs of singular values ​​in descending order to form a descending sequence: ， In the above formula, ~ These are the elements from 1 to m in the descending sequence; S204, Calculate the relative ratio sequence based on the results after sorting in descending order: ， In the above formula, ~ Let be the logarithms of the sums of the first to m pairs of singular values, and we have: ， In the above formula, Let be the logarithm of the sum of the i-th pair of singular values. Let i be the i-th element in the descending sequence. It is the first element in the descending sequence; S205, Take the largest index value from the relative ratio sequence as the number of frequency components obtained. p The estimated value .

3. The high-precision amplitude and phase estimation method suitable for rapidly dynamically changing signals according to claim 1, characterized in that, Step S104 includes: S301, based on the number of frequency components p The estimated value And right singular matrix Extract the right singular matrix according to the following formula. Two subspaces : ， ， In the above formula, A right singular matrix Column 1 A right singular matrix Column 2 A right singular matrix The List, A right singular matrix The Column; extract the right singular matrix according to the following formula. Two subspaces : ， ， In the above formula, A right singular matrix Column 2 A right singular matrix Column 3 A right singular matrix The List, A right singular matrix The List; S302, combined with subspace Matrix Q1 is obtained using the least squares method based on the following formula: ， In the above formula, To solve for the eigenvalues ​​and eigenvectors, For subspace The pseudo-inverse matrix, For subspace The conjugate transpose of . These are the eigenvalues ​​and eigenvectors corresponding to the first frequency component. The eigenvalues ​​and eigenvectors corresponding to the second frequency component are: Let i be the eigenvalue and eigenvector corresponding to the i-th frequency component. For the first The eigenvalues ​​and eigenvectors corresponding to each frequency component; combined with the subspace Matrix Q2 is obtained using the least squares method based on the following formula: ， In the above formula, For subspace The pseudo-inverse matrix, For subspace The conjugate transpose of; These are the eigenvalues ​​and eigenvectors corresponding to the first frequency component. The eigenvalues ​​and eigenvectors corresponding to the second frequency component are: Let i be the eigenvalue and eigenvector corresponding to the i-th frequency component. For the first Eigenvalues ​​and eigenvectors corresponding to each frequency component; S303, combine matrix Q1 and matrix Q2 to form a diagonal matrix And diagonal matrix The eigenvalues ​​corresponding to each frequency component are the average of the corresponding eigenvalues ​​in matrices Q1 and Q2.

4. The high-precision amplitude and phase estimation method suitable for rapidly dynamically changing signals according to claim 1, characterized in that, Step S105 uses a diagonal matrix Eigenvalues ​​construct system matrix A The function expression is: ， ， ， ， In the above formula, Let y(n) be the length of the discrete signal sequence. Represents a diagonal matrix The first eigenvalue, Represents a diagonal matrix The second eigenvalue, Represents a diagonal matrix The 1 eigenvalue, The attenuation factor of the fundamental phasor. The imaginary unit, The angular frequency of the fundamental phasor; This is the attenuation factor corresponding to the second frequency component. This is the angular frequency corresponding to the second frequency component; For the first The attenuation factor corresponding to each frequency component For the first The angular frequencies corresponding to each frequency component, and have , , Let be the attenuation factor for the k-th frequency component. Let be the angular frequency of the k-th frequency component. The sampling period is related to the sampling frequency. satisfy , Let be the frequency magnitude of the k-th frequency component.

5. The high-precision amplitude and phase estimation method suitable for rapidly dynamically changing signals according to claim 1, characterized in that, In step S106, the system matrix is ​​used. A The functional expression for calculating the dynamic phasor coefficient vector S using the least squares method for the discrete signal sequence y(n) is as follows: ， In the above formula, For the system matrix A The conjugate transpose of . The discrete signal sequence is y(n); and the resulting computational dynamic phasor coefficient vector S has the following form: ， In the above formula, These are the complex constants obtained by least squares fitting. This indicates transpose.

6. The high-precision amplitude and phase estimation method suitable for rapidly dynamically changing signals according to claim 5, characterized in that, In step S107, the dynamic phasor coefficient vector S and the damping coefficient are used. Calculate the amplitude X and phase of a continuous signal. The function expression is: ， ， ， ， ， ， In the above formula, The amplitude of the fundamental phasor of the dynamic signal. The phase of the fundamental phasor of the dynamic signal. The rate of change of the fundamental phasor amplitude of the dynamic signal. The frequency of the fundamental phasor of the dynamic signal. The rate of change of the amplitude of the fundamental phasor of the dynamic signal. The rate of change of the frequency of the fundamental phasor of the dynamic signal. and These represent the real and imaginary parts of a complex number, respectively. The damping coefficient of the fundamental phasor; , and is a complex constant obtained by least squares fitting in the dynamic phasor coefficient vector S.

7. The high-precision amplitude and phase estimation method suitable for rapidly dynamically changing signals according to claim 1, characterized in that, The Hankel matrix constructed in step S102 Y The function expression is: ， In the above formula, For the first discrete signal, The second discrete signal, For the first A discrete signal For the first A discrete signal For the first A discrete signal Let y(n) be the length of the discrete signal sequence. In step S102, the Hankel matrix is... Y Singular value decomposition consists of a left singular matrix U, a diagonal matrix D, and a right singular matrix. The function expression is: ， In the above formula, It is a left singular matrix. It is a diagonal matrix. It is a right singular matrix. For singular value decomposition, It is a Hankel matrix, and a left singular matrix. diagonal matrix Right singular matrix The size is , Let y(n) be the length of the discrete signal sequence.

8. A high-precision amplitude and phase estimation device suitable for rapidly changing signals, comprising a microprocessor and a memory interconnected, characterized in that, The microprocessor is programmed or configured to execute the high-precision amplitude and phase estimation method for rapidly changing signals as described in any one of claims 1 to 7.

9. A computer-readable storage medium storing a computer program, characterized in that, The computer program is used to be programmed or configured by a microprocessor to execute the high-precision amplitude and phase estimation method for rapidly changing signals as described in any one of claims 1 to 7.

Images