Harmonic fast detection system based on improved PASTD-ESPRIT

Abstract

The application provides a harmonic rapid detection system based on improved PASTD-ESPRIT, which comprises a preprocessing module based on QR rank revealing, performs column-pivot QR decomposition on a constructed Hankel data matrix, extracts a diagonal element sequence of an upper triangular matrix R, determines an effective order r of a signal subspace by using a second-order difference criterion of the diagonal elements, and transfers the obtained order to a PASTD-ESPRIT main detection module in real time; the PASTD-ESPRIT main detection module initializes a weight matrix according to the effective order r, updates a characteristic vector of the signal subspace in a way of a projection approximation subspace tracking technology PASTD recursion, solves frequencies of each harmonic by using an ESPRIT (Estimation of Signal Parameters via Rotational Invariance Techniques), completes joint estimation of harmonic amplitudes and phases by using a least square method, and finally outputs online closed-loop detection results of three parameters of frequencies, amplitudes and phases.

Description

Technical Field

[0001] This invention relates to the field of AC power generation, transmission, distribution and consumption technology, and specifically to a rapid harmonic detection system based on an improved PASTD-ESPRIT. Background Technology

[0002] The widespread integration of nonlinear loads and distributed generation into power systems has led to increasingly prominent problems of current and voltage waveform distortion in the power grid, with harmonic pollution becoming a major factor affecting power quality. Accurate and rapid detection of the amplitude, frequency, and phase of harmonic components is a crucial prerequisite for achieving harmonic mitigation, reactive power compensation, and system protection. Especially under complex operating conditions such as power grid frequency fluctuations, background noise interference, and time-varying harmonic components (e.g., interharmonics and transient harmonics), stringent requirements are placed on the speed of detection algorithms and the accuracy of parameter identification.

[0003] Currently, commonly used harmonic detection techniques include Fast Fourier Transform (FFT), Prony analysis, and modern spectral estimation methods based on subspace decomposition (such as the rotation-invariant subspace method, ESPRIT). While FFT is computationally simple, it suffers from spectral leakage and the picket-fence effect, limiting accuracy under asynchronous sampling. Furthermore, its time and frequency resolutions are mutually restrictive, making it difficult to meet the requirements of dynamic harmonic tracking. Prony analysis can directly identify harmonic parameters, but it is extremely sensitive to noise, and the choice of model order has a significant impact on the results, lacking robustness. The ESPRIT algorithm, with its advantages of super-resolution and no need for peak search, has been introduced into harmonic detection. However, the traditional ESPRIT algorithm requires singular value decomposition (SVD) of the Hankel matrix or eigenvalue decomposition of the covariance matrix, resulting in a computational complexity as high as O(N³), making real-time updates between samples difficult to achieve on embedded systems or DSPs. Simultaneously, its batch processing method cannot recursively utilize historical information, leading to significant detection delays or even tracking failures when harmonic parameters change abruptly or over time.

[0004] To overcome the problems of poor speed and high computational cost of batch processing ESPRIT, existing improved schemes introduce Projective Approximate Subspace Tracking (PAST) and its variant PASTD algorithm. PASTD estimates the signal subspace recursively, avoiding direct eigenvalue decomposition and significantly reducing the computational burden. However, the standard PASTD-ESPRIT scheme still has two inherent defects: First, PASTD uses a fixed forgetting factor, which leads to a contradiction between the noise suppression capability of steady-state harmonic detection and the response speed of dynamic harmonic tracking—a small forgetting factor, while fast, easily amplifies noise disturbances, increasing the estimation variance; a large forgetting factor, while smoothing noise, lags in updating when harmonic components change rapidly, resulting in significant tracking errors. Second, and more critically, this scheme lacks an online estimation and update mechanism for the effective rank of the signal subspace (i.e., the order of the harmonic model). In actual power systems, the number of harmonics (fundamental, integer harmonics, and interharmonics) often changes dynamically with load variations, while existing PASTD-ESPRIT schemes typically preset a fixed order or rely on offline coarse estimation. When noise levels are high or the signal-to-noise ratio fluctuates, PASTD's estimation of eigenvalues ​​is prone to "order drift," leading to angle mismatches when ESPRIT solves the rotationally invariant equations, resulting in frequency estimation errors or spurious harmonic components. This problem severely limits the algorithm's robustness and practicality in complex time-varying environments.

[0005] In summary, existing harmonic detection schemes based on PASTD-ESPRIT have significant shortcomings in terms of subspace tracking adaptability, real-time model order update capability, and robustness in high-noise environments, making it difficult to simultaneously achieve both speed and accuracy in detection.

[0006] Literature: Feng Liu, Ligong Sun, Zhihao Cheng. The Improved PASTD andESPRIT in Power System Harmonic Estimation Application[C]. 2nd International Conference on Electrical, Control and Automation Engineering (ECAE 2017), 2017: 1-5. The proposed technical solution still has the following shortcomings: First, its model order estimation relies entirely on the eigenvalue sequence updated during the PASTD iteration process and the MDL criterion for determination. When the signal-to-noise ratio is low or there are interharmonics with similar frequencies in the signal, the eigenvalue decay law is not obvious, which can easily lead to order estimation deviation. The order deviation will directly affect the correctness of the subspace dimension, thereby reducing the frequency estimation accuracy of ESPRIT. Second, order estimation and subspace tracking are in the same iterative loop and are tightly coupled. When the harmonic components in the signal change dynamically (such as a sudden increase or decrease in harmonic components), the MDL criterion needs a certain amount of time to accumulate before it can respond to the order change. During this transient period, the subspace dimension and the actual signal dimension are mismatched, which leads to an increase in transient estimation error and limits the fast tracking performance. In addition, the order estimation in this paper relies entirely on the eigenvalue information in the iteration process and lacks a preprocessing mechanism to use the structural features of the observed data for order prediction. This makes it take a long time for the algorithm to converge to the correct order in the initial stage, affecting the speed of fast detection.

[0007] Patent application CN106405487A discloses a general spatial spectrum estimation method based on extended ESPRIT technology. This method introduces a virtual reference array, whose position is determined by the element positions of the actual array. A phase compensation matrix is ​​defined by the positional difference between corresponding elements in the actual and virtual arrays and a phase compensation angle. The phase compensation matrix is ​​a diagonal unitary matrix. This unitary matrix is ​​used to perform phase compensation on the original signal feature space. Then, the classic ESPRIT algorithm is applied to the compensated signal feature space to obtain multiple angle estimates, but only one angle estimate is output. This angle is closest to the input phase compensation angle, and the eigenvalues ​​of the fitted matrix corresponding to this angle estimate are also output. The newly defined general ESPRIT spatial spectrum can be calculated using the output angle estimate and the corresponding eigenvalues. This general ESPRIT spatial spectrum has a spectral peak at the actual signal incident angle. Therefore, searching for the spectral peak of the general ESPRIT spatial spectrum in the parameter space yields the direction-of-arrival (DOA) estimate of the spatial signal. However, this patent cannot completely solve the existing technical problems and does not meet the requirements of this invention. Summary of the Invention

[0008] In view of the deficiencies in the prior art, the purpose of this invention is to provide a rapid harmonic detection system based on the improved PASTD-ESPRIT.

[0009] The fast harmonic detection system based on the improved PASTD-ESPRIT provided by the present invention includes a preprocessing module based on QR rank and a PASTD-ESPRIT main detection module. The input of the QR-based preprocessing module receives the three-phase voltage or current sampling signal after coordinate transformation, and its output is connected to the order input of the PASTD-ESPRIT main detection module. The QR-based preprocessing module performs column pivoting QR decomposition on the constructed Hankel data matrix, extracts the diagonal element sequence of the upper triangular matrix R, and uses the second-order difference criterion of the diagonal elements to determine the effective order r of the signal subspace. The obtained order is then transmitted to the PASTD-ESPRIT main detection module in real time. The PASTD-ESPRIT main detection module initializes the weight matrix based on the effective order r, updates the feature vector of the signal subspace using the PASTD recursive method of the compact projection approximate subspace tracking technique, then uses the rotation-invariant technique ESPRIT to solve for the frequency of each harmonic, and completes the joint estimation of harmonic amplitude and phase through the least squares method, finally outputting the online closed-loop detection results of the three parameters of frequency, amplitude and phase.

[0010] Preferably, the Hankel data matrix is ​​constructed as follows: A discrete-time series is obtained by sampling a continuous signal:

[0011] In the formula, It is the signal value of the nth discrete sampling point. N This represents the total number of sampling points. T s The sampling interval is... r The number of harmonics contained in the signal. e ( n () has a mean of 0 and a variance of σ 2 Gaussian white noise, A i 、f i φ i The first i The amplitude, frequency, and initial phase of the subharmonic; ω i = 2πf i , is the angular frequency of the i-th harmonic; j is the imaginary unit; Take the length as m Using a sliding window, construct a data vector:

[0012] in, m > r superscript T Indicates the conjugate transpose, the vector is decomposed into signal components. X (n ) and noise components E ( n The expression is: Y ( n )= BX ( n )+ E ( n ),in B The Vandermonde matrix is ​​composed of harmonic frequencies. X ( n It contains the complex envelope of each harmonic; the corresponding autocorrelation matrix R is:

[0013] in, D = E [ X ( n ) X H ( n )], superscript H denoted as conjugate transpose, where I is the identity matrix.

[0014] Preferably, the column-major QR decomposition satisfies: for any m × n real matrix of order A and m ≥ n There exists a column permutation matrix Π such that the permuted matrix is... A Π is obtained through standard QR decomposition. A Π= QR Where Q is an orthogonal matrix, the upper triangular matrix R The blocks are as follows:

[0015] in, R 11 ∈C r×r The rank information corresponding to the signal subspace, and its minimum singular value reflects the lower bound of the energy of the principal components of the signal; R 22 ∈C (n-r)×(n-r) Its spectral norm reflects the upper bound of the energy of the noise subspace; R 11 ∈C r×(n-r) The top right corner sub-block, 0 represents (n The zero matrix of r)×r.

[0016] Preferably, the method for extracting the diagonal element sequence of the upper triangular matrix R is as follows: Perform QR decomposition of the Hankel data matrix using column-selective pivoting, extract the main diagonal elements of matrix R, and arrange them in descending order to obtain the pivot sequence of the QR decomposition. d i =| r ii |, i =1,2…, n -1

[0017] in, r ii This represents the i-th main diagonal element of the upper triangular matrix R obtained after performing column-selective pivoting QR decomposition on the Hankel data matrix; d i It is the absolute value of the diagonal element, forming the principal component sequence; The initial order of the model; The amplitudes of the first few diagonal elements are larger values ​​within a preset range, corresponding to the main components in the signal. As the index increases, the amplitude changes abruptly at a certain position r, dropping to a smaller value within the preset range. This point of change is the critical boundary between the signal and noise.

[0018] Preferably, the second-order difference criterion for diagonal elements is: through second-order difference... To depict its changing trend, among which d ( k ) is the first k The size of the main diagonal elements is expressed as:

[0019] along with k The value asymptotically approximates the initial order of the model. r c The absolute value of the second difference | | will continue to decrease, and when this value converges to zero, then:

[0020] Corresponding k Located at the point of abrupt change in the curve, this position corresponds to the critical boundary between signal and noise, and is also the effective order. r .

[0021] Preferably, the process of initializing the weight matrix is ​​as follows: Set the number of iterations k =0, maximum number of iterations L = N m +1, the weight matrix is ​​initialized as follows:

[0022] in, The initial values ​​for the weight matrix are m×r; Simultaneously initialize eigenvalues ​​and eigenvectors: λ i (0)=1 indicates that the initial value of the i-th eigenvalue is set to 1. w i (0)= W i (0) Initial weighting matrix of the i-th eigenvector The i-th column, i ∈[1, r ]; Iteration from k Starting with =1, first assign the current sampled data vector the value: Y(k) is a data vector constructed from the sampled signal.

[0023] Preferably, the PASTD recursive method is for i =1 to i = r Perform the following calculations in sequence: Find the projection of the current data onto the existing feature vectors:

[0024] Among them, scalar projection This indicates the magnitude of the component of the current data in the direction of the i-th feature vector; This represents the current data vector; Using the forgetting factor β Update eigenvalues, 0 < β <1. Balancing the impact of historical information and new data:

[0025] in, Represents the feature value at the current moment; Update the feature vector using residuals and projections:

[0026]

[0027] Among them, residual This represents the portion of the current data that is not represented by the i-th feature vector; This represents the updated feature vector; Update the data vector to remove the influence of the current feature vector and prepare for the next feature vector estimation:

[0028] After the iteration is complete, the updated feature vectors are arranged into a weight matrix column by column. W =[ w 1,..., wr This matrix is ​​an approximate estimate of the signal subspace.

[0029] Preferably, the process of solving for the frequencies of each harmonic using rotational invariance techniques is as follows: weight matrix W Perform block processing and remove W The last row yields the submatrix. W e ∈C (m 1)×r Remove W The first row yields the upper submatrix. W f ∈C (m 1)×r Constructing rotation operators S ∈C r×r It satisfies the rotation invariance relation:

[0030] Solving for the rotation operator using the least squares method:

[0031] Eigenvalue decomposition is performed on the rotation operator S to obtain the eigenvalues. z i , i =1,2,..., r Its mapping relationship with harmonic frequencies is as follows:

[0032] The frequencies of each harmonic can be obtained by reverse calculation:

[0033] In the formula, arg() represents the argument of a complex number, and its value ranges from -1 to 1. π , π The value is adjusted based on the harmonic frequency range of the power system.

[0034] Preferably, the process of jointly estimating the harmonic amplitude and phase using the least squares method is as follows: Obtaining harmonic frequencies f i Then, construct the observation matrix Φ∈C N×r :

[0035] Organize the sampling sequence into an observation vector y =[ y (0), y (1),...,y ( N 1)] T Construct a least squares optimization problem:

[0036] Where A=[ A 1 , A 2 ,..., A r ] T Given a complex amplitude vector, solve using canonical equations:

[0037] For complex amplitude C i =A i Decompose the data to obtain the amplitude and phase:

[0038] in, For the first r The initial phase of the subharmonic, Φ is the observation matrix, and Re() and Im() represent taking the real and imaginary parts, respectively.

[0039] Preferably, during the PASTD iteration process, the order adjustment is triggered based on the eigenvalue change pattern, and an online adaptive update mechanism for the subspace dimension is established: when the harmonic components change dynamically, the column pivot QR decomposition and diagonal second-order difference criterion are re-executed through the QR rank preprocessing module to obtain the updated effective order, and the updated order is transmitted to the PASTD-ESPRIT main detection module in real time to reinitialize the dimension of the weight matrix and the number of eigenvectors in the PASTD recursion process.

[0040] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention directly obtains the effective dimension of the signal subspace through the rank-revelation process of QR decomposition, providing a rigorous theoretical basis and numerically stable order input for the PASTD-ESPRIT algorithm, avoiding the computational redundancy and misjudgment risk caused by the traditional information theory criteria (such as AIC, MDL) requiring multiple feature decompositions or subjective threshold settings. (2) The preprocessing structure based on QR rank-revelation given in this invention is independent of the subsequent subspace tracking calculation. It only triggers parameter adjustment when the order changes significantly. This ensures the high dynamic response capability of harmonic detection to load switching and sudden changes in system operating conditions. It also significantly reduces the noise subspace pollution caused by overestimation of the order and the frequency component loss caused by underestimation of the order. It is far superior to the traditional harmonic detection strategy with fixed order or a posteriori adjustment in terms of detection accuracy and fast tracking speed. (3) This invention adjusts only the model order without changing the core recursive structure of PASTD-ESPRIT, which has good algorithm compatibility and lightweight hardware implementation advantages. It is particularly suitable for embedded application scenarios that are sensitive to computing resources and require fast response, such as power quality monitoring terminals and active power filter instruction current extraction. Attached Figure Description

[0041] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 A schematic diagram of the PASTD-ESPRIT algorithm; Figure 2 The charts show a comparison of the order determination results based on SVD and QR, respectively. Figure 3 To improve the PASTD-ESPRIT algorithm flowchart; Figure 4a This is a graph showing the frequency estimation results. Figure 4b This is a diagram showing the amplitude estimation result. Figure 4c This is a diagram showing the phase estimation results; Figure 5 The graph shows the effect of transient response frequency estimation. Detailed Implementation

[0042] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that those skilled in the art can make several changes and improvements without departing from the concept of the present invention. These all fall within the scope of protection of the present invention.

[0043] Example To achieve rapid and high-precision detection of harmonics in power systems, existing solutions incorporate the PASTD-ESPRIT algorithm to reduce the computational burden of the traditional ESPRIT. However, this approach suffers from a critical flaw: the model order (i.e., the effective rank of the signal subspace, corresponding to the number of harmonic components) is typically preset to a fixed value or relies solely on offline coarse estimation, failing to track the dynamic changes of actual harmonic components online. The number of harmonics changes in real time when load switching, inverter operating conditions change, or interharmonics occur in the power system. A fixed order leads to subspace dimension mismatch—too small an order will miss true harmonic components (false detection), while too large an order will introduce noise subspace components, generating spurious harmonics (false detection), and simultaneously cause angular mismatch in the ESPRIT algorithm when solving rotationally invariant equations, resulting in a significant deviation of frequency estimation from the true value. Furthermore, the existing PASTD-ESPRIT lacks a preprocessing mechanism for the input data matrix. When the signal-to-noise ratio is low or the matrix is ​​ill-conditioned, fixed-order subspace tracking is prone to "order drift," further degrading detection performance. To address the shortcomings of the prior art, this invention proposes a rapid harmonic detection system based on an improved PASTD-ESPRIT.

[0044] To address the shortcomings of existing technologies, this invention proposes a fast harmonic detection system based on an improved PASTD-ESPRIT algorithm. The core improvement lies in the introduction of a rank reveal preprocessing module based on QR decomposition. This module performs orthogonal triangulation (QR decomposition) on the rapidly sampled signal data matrix. On one hand, this improves the matrix condition number and filters out some noise disturbances through orthogonal transformation, thereby enhancing numerical stability. On the other hand, it utilizes a rank reveal mechanism (such as column-dominant QR decomposition) to quickly estimate the effective rank of the data matrix, i.e., dynamically determining the number of harmonic components (model order) at the current moment. This dynamic order serves as a core parameter, directly guiding the subspace dimension setting in the PASTD-ESPRIT algorithm, ensuring accurate matching between the subspace tracking and the dimensions of the actual harmonic components.

[0045] Compared with the existing PASTD-ESPRIT scheme, the improved strategy proposed in this invention suffers from several drawbacks: the existing scheme, due to its fixed order, suffers from missed detections or spurious harmonics when the number of harmonics changes; while this invention, through a QR rank reveal preprocessing module, can adaptively update the model order online, providing accurate order guidance for PASTD-ESPRIT and fundamentally solving the "order mismatch" problem. This invention effectively overcomes the technical bottleneck of the inability to dynamically update the model order in existing technologies, significantly improving the detection accuracy and robustness of the algorithm in time-varying harmonic environments. This scheme has the following advantages: The orthogonal triangulation process of the QR preprocessing module improves the ill-conditioning of the input data matrix and suppresses noise disturbances, providing more robust initial conditions for subsequent subspace tracking, reducing fluctuations in order estimation under low signal-to-noise ratio environments, and enhancing the numerical stability and noise robustness of the overall algorithm. By utilizing the rank revelation mechanism of QR decomposition, the number of harmonic components (including fundamental, integer harmonics, and interharmonics) can be quickly estimated without prior information or offline training. The dynamically updated accurate order is then passed to the PASTD-ESPRIT algorithm to avoid missed detections due to too small an order or false harmonics caused by too large an order, thus ensuring the accuracy of frequency, amplitude, and phase estimation.

[0046] This invention provides a fast harmonic detection system based on an improved PASTD-ESPRIT algorithm. The overall structure of the improved PASTD-ESPRIT harmonic detection system is shown in Figure 4, which includes a QR-rank-based preprocessing module and a PASTD-ESPRIT main detection module. The QR-rank-based preprocessing module receives the three-phase voltage / current sampling signals after coordinate transformation. It performs column-pivoting QR decomposition on the constructed Hankel data matrix to extract the diagonal element sequence of the upper triangular matrix R. It then automatically determines the effective order r of the signal subspace using the second-order difference criterion for the diagonal elements and transmits the obtained order to the PASTD-ESPRIT algorithm in real time. The PASTD-ESPRIT main detection module initializes the weight matrix W based on this order r, updates the eigenvectors of the signal subspace using the PASTD (Plain Projection Approximate Subspace Tracking) recursive method, and then uses the rotation-invariant technique (ESPRIT) to solve for the frequency of each harmonic. Finally, it uses the least squares method to jointly estimate the harmonic amplitude and phase, achieving online closed-loop detection of the three parameters: frequency, amplitude, and phase.

[0047] Compared with existing rapid harmonic detection methods, the rapid harmonic detection system based on the improved PASTD-ESPRIT proposed in this patent can achieve adaptive estimation and dynamic updating of the harmonic model order without relying on preset order prior information and complex eigenvalue decomposition, fundamentally improving the accuracy and robustness of harmonic detection under time-varying conditions. The traditional PASTD-ESPRIT algorithm requires manually setting a fixed model order. When harmonic sources are switched or their components change in the actual power grid, order mismatch will directly lead to incorrect subspace partitioning and serious deviations in parameter estimation. Although some improved methods introduce information theory criteria or singular value decomposition for order estimation, they are computationally burdensome, have poor speed, and the order criterion is unstable in high-noise environments.

[0048] This invention proposes a fast harmonic detection system based on QR rank revealing preprocessing and adaptive order PASTD-ESPRIT, focusing on the online identification of subspace dimension and the linkage of main algorithm parameters. First, this invention utilizes column-pivoting QR decomposition to fully map the singular value attenuation characteristics of the received data matrix to the diagonal element sequence of the upper triangular matrix, achieving rank reveal capability equivalent to singular value decomposition with extremely low computational overhead, providing an efficient numerical foundation for order estimation. Then, based on the second-order difference convergence characteristics of the diagonal element sequence, an automatic order criterion is designed, accurately identifying the critical boundary between the signal subspace and the noise subspace under noise interference without manual threshold setting, providing a reliable and physically consistent initial order value for the PASTD recursive process. Furthermore, an online adaptive subspace dimension update mechanism is established, triggering order adjustment based on eigenvalue changes during PASTD iteration, enabling the detection strategy to respond quickly to dynamic changes in harmonic components. Specific implementation methods are as follows: A. Principles of the PASTD-ESPRIT Algorithm In the scenario of harmonic detection in power systems, a discrete time series is obtained after sampling a continuous signal: (1) In the formula, It is the signal value of the nth discrete sampling point. N This represents the total number of sampling points. T s The sampling interval is... r The number of harmonics contained in the signal. e ( n () has a mean of 0 and a variance of σ 2 Gaussian white noise, A i 、f i φ i The first i The amplitude, frequency, and initial phase of the subharmonic. ω i = 2πf i , is the angular frequency of the i-th harmonic; n ∈[0, N 1], j is the imaginary unit.

[0049] To adapt to the subspace tracking algorithm, the length is taken as... m ( m > r Using a sliding window, construct a data vector: (2) superscript TIndicates the conjugate transpose, the vector is decomposed into signal components. X ( n ) and noise components E ( n The expression is: Y ( n )= BX ( n )+ E ( n ),in B The Vandermonde matrix is ​​composed of harmonic frequencies. X ( n It contains the complex envelope of each harmonic. E ( n () represents the noise vector. The corresponding autocorrelation matrix is: (3) in, D = E [ X ( n ) X H ( n E[ ] represents the expected value, and the superscript indicates the expected value. H This represents the conjugate transpose, where I is the identity matrix. σ 2 It is Gaussian white noise e ( n The variance of the matrix is ​​calculated. The principal eigenvectors of this matrix span the corresponding signal components in the subspace, which is the core foundation for subsequent subspace tracking.

[0050] The Projection Approximation Subspace Tracking (PAST) algorithm, proposed by Yang in 1995, departs from the traditional eigenvalue decomposition framework, transforming subspace decomposition into an unconstrained optimization problem. It approximates the signal subspace by minimizing the objective function, which takes the form: (4) in, tr ( ) represents the matrix trace operation, and E is the expected value. W ∈ C m×r Let be the weight matrix to be updated. As can be seen from the expansion, this function essentially measures the original data vector. Y Rather than WW H The expected value of the squared error of the projection onto the Zhang Cheng space provides convenience for gradient iteration.

[0051] The PAST algorithm has limitations when tracking multiple feature vectors. The compressed projective approximate subspace tracking algorithm (PASTD, PAST Deflated) improves upon this, offering convergence performance and numerical stability better suited for practical engineering applications. Its core idea is to first iteratively obtain the principal feature vector using the PAST framework, then remove the projection component of this vector onto the data. This process is repeated in the remaining data until all features are tracked. r The estimation of eigenvectors, the algorithm flow is as follows: Figure 1 As shown.

[0052] Set the number of iterations k =0, maximum number of iterations L = N m +1, the weight matrix is ​​initialized as follows: (5) in, The initial values ​​for the weight matrix are m×r; Simultaneously initialize eigenvalues ​​and eigenvectors: λ i (0)=1 indicates that the initial value of the i-th eigenvalue is set to 1. w i (0)= W i (0) Initial weighting matrix of the i-th eigenvector The i-th column, i ∈[1, r Iteration from k Starting with =1, first assign the current sampled data vector the value: (6) against i =1 to i = r Perform the following calculations in sequence: Find the projection of the current data onto the existing feature vectors: (7) Among them, scalar projection This indicates the magnitude of the component of the current data in the direction of the i-th feature vector; This represents the current data vector; Using the forgetting factor β (0< β <1) Update feature values, taking into account both historical information and the impact of new data: (8) in, Represents the feature value at the current moment; Updating feature vectors using residuals and projections ε (9) (10) Among them, residual This represents the portion of the current data that is not represented by the i-th feature vector; This represents the updated feature vector; Update the data vector to remove the influence of the current feature vector and prepare for the next feature vector estimation: (11) After the iteration is complete, the updated feature vectors are arranged into a weight matrix column by column. W =[ w 1,..., w r This matrix is ​​an approximate estimate of the signal subspace. The ESPRIT algorithm utilizes the rotation invariance of the signal subspace to achieve high-precision frequency estimation without peak searching. Combined with the least squares method, it can further solve for the amplitude and phase.

[0053] First, the weight matrix... W Perform block processing: remove W The last row yields the submatrix. W e ∈C (m 1)×r Remove W The first row yields the upper submatrix. W f ∈C (m 1)×r Constructing rotation operators S ∈C r×r It satisfies the rotation invariance relation: (12) Solving for the rotation operator using the least squares method: (13) Eigenvalue decomposition is performed on the rotation operator S to obtain the eigenvalues. z i ( i =1,2,..., r ), which satisfies a mapping relationship with harmonic frequencies: (14) The frequencies of each harmonic can be obtained by reverse calculation: (15) In the formula, arg() represents the argument of a complex number, and its value ranges from -1 to 1. π , π The values ​​need to be corrected in conjunction with the harmonic frequency range of the power system (e.g., the fundamental frequency is 50Hz, and the harmonics are integer multiples of the fundamental frequency) to ensure that the frequency values ​​conform to the actual physical meaning.

[0054] Obtaining harmonic frequencies f i Then, construct the observation matrix Φ∈C N×r : (16) Organize the sampling sequence into an observation vector y =[ y (0), y (1),..., y ( N 1)] T Construct a least squares optimization problem: (17) Where A=[ A 1 , A 2 ,..., A r ] T It is a complex amplitude vector. Solved using canonical equations: (18) For complex amplitude C i =A i Decompose the data to obtain the amplitude and phase: (19) in, For the first r The initial phase of the subharmonic, Φ is the observation matrix, and Re() and Im() represent taking the real and imaginary parts, respectively.

[0055] This method combines the subspace tracking advantage of the PASTD algorithm with the parameter estimation advantage of the ESPRIT algorithm, ensuring both speed and high-precision synchronous estimation of harmonic frequency, amplitude, and initial phase.

[0056] B. Improved PASTD-ESPRIT Algorithm—An Automatic Order Determination Method Based on QR Decomposition of Diagonal Element Sequences Sub-model order rThe accurate selection of the order is crucial for the application of the PASTD algorithm, directly determining the modeling accuracy, computational efficiency, and reliability of the analysis results. A reasonable order selection allows the algorithm to accurately capture the core features of the data, avoiding underfitting or overfitting.

[0057] As shown in Table 1, when detecting 67Hz, 167Hz, and 267Hz multi-frequency signals under 20dB white noise interference, and with other parameters remaining constant, different orders... r The order of a frequency estimation method has a significant impact on the results: an excessively high order can introduce redundant components and increase computational costs; an excessively low order will fail to accurately estimate the frequency, leading to biased results. Therefore, determining the optimal order is crucial. r Balancing modeling accuracy and computational efficiency is an important prerequisite for fully leveraging the performance of the PASTD algorithm.

[0058] Table 1 Comparison of frequency detection performance and computation time of the PASTD-ESPRIT algorithm at different model orders.

[0059] In view of this, this application introduces a preprocessing mechanism based on rank-revealing QR decomposition. By extracting structured rank information from the sampling covariance matrix, it provides adaptive order guidance parameters for the PASTD-ESPRIT algorithm, thereby avoiding the performance limitations caused by the fixed order assumption. The core idea of ​​this rank-revealing QR decomposition is to completely map the singular value decay characteristics of the original matrix A to the block structure of the upper triangular matrix R through a specific column permutation matrix. The specific principle is as follows: The core idea of ​​rank-based QR decomposition lies in transforming the original matrix through a specific column permutation matrix. A The singular value decay feature is fully mapped to the upper triangular matrix. R In the block structure. For any m × n real matrix of order A ( m ≥ n Given a column permutation matrix Π, there exists a matrix such that the resulting matrix is... A Π is obtained through standard QR decomposition. A Π= QR , upper triangular matrix R The blocks are as follows: (20) in, R 11 ∈C r×r The rank information corresponding to the signal subspace, and its minimum singular value reflects the lower bound of the energy of the principal components of the signal; R 22 ∈C (n-r)×(n-r)Its spectral norm reflects the upper bound of the energy of the noise subspace; R 11 ∈C r×(n-r) The top right corner sub-block, 0 represents (n A zero matrix of r × r. Mathematically, it has been proven that if the optimal column permutation is chosen, the sub-block... R 11 and R 22 Singular values ​​and the original matrix A The singular values ​​satisfy strict quantization bound relations: (twenty one) (twenty two) The above relationship reveals a key conclusion: when the original moments Formation A In the r There is a significant order of magnitude gap between the first and (r+1)th singular values ​​(i.e. σ r (A) σ r+1 When (A) is reached, the gap will be reflected without loss in the result of QR decomposition—that is... R 11 The minimum singular value will be much greater than R 22 spectral norm ( σ min ( R 11 ) || R 22 ||2). This is the mathematical basis for how rank reveals QR decomposition can replace SVD in determining the numerical rank.

[0060] In practical calculations, direct solution is used. σ min ( R 11 ) and || R 22 ||2 still involves high computational overhead. However, thanks to the introduction of column-major QR decomposition, the decision process can be greatly simplified.

[0061] The column-pivoted QR decomposition employs a "greedy strategy" during orthogonalization—always selecting the column with the largest norm among the remaining columns as the current pivot for elimination. This operation forcibly guarantees the generation of an upper triangular matrix. R It has the property that the absolute value of its diagonal elements is non-increasing, that is: (twenty three) More importantly, for an upper triangular matrix, its diagonal elements rii It is a good approximation of the singular values ​​of the matrix, especially the minimal diagonal elements at the ends of the matrix, whose order of magnitude directly determines the value. R 22 The magnitude of the spectral norm.

[0062] Based on the rank information representation characteristics inherent in the aforementioned matrix block structure, the effective rank of the original matrix can be determined by quantitatively analyzing the diagonal elements of the upper triangular matrix R. More importantly, for an upper triangular matrix, its diagonal elements... r ii It is a good approximation of the singular values ​​of the matrix, especially the minimal diagonal elements at the ends of the matrix, whose order of magnitude directly determines the value. R 22 The magnitude of the spectral norm. Applying this property to the Hankel matrix constructed from the observed signals, a QR decomposition with column-pivoted elements is performed on the Hankel matrix to extract... R Arranging the main diagonal elements of the matrix in descending order yields the principal component sequence of the QR decomposition. The specific principle is as follows: Perform column-pivoting QR decomposition on the Hankel matrix to extract... R Arrange the main diagonal elements of the matrix in descending order to obtain the principal component sequence of the QR decomposition: (twenty four) in, d i =| r ii |, i =1,2…, n -1, r ii This represents the i-th main diagonal element of the upper triangular matrix R obtained after performing column-selective pivoting QR decomposition on the Hankel data matrix; d i It is the absolute value of the diagonal element, forming the principal component sequence; This represents the initial order of the model.

[0063] The main diagonal elements of the upper triangular matrix obtained by QR decomposition | r ii There are significant differences in numerical values: the first few diagonal elements have larger amplitudes, corresponding to the main components of the signal; as the index increases, the amplitude will decrease at a certain position. r A sudden change occurs at a certain point, causing the signal to drop rapidly from a large value to a very small value. This point of change is the critical boundary between signal and noise.

[0064] After this critical point, the amplitude changes of the main diagonal elements tend to level off, and the overall value is extremely small, reflecting noise interference. Since the main diagonal element sequence can be regarded as a curve composed of discrete points, it can be obtained through second-order difference. To depict its changing trend, among which d ( k ) is the first k The size of the main diagonal elements, i.e.: (25) along with k The value asymptotically approximates the initial order of the model. r c The absolute value of the second difference | | will continue to decrease, until the value converges to zero, that is: (26) Corresponding k It is located precisely at the point of abrupt change in the curve, which corresponds to the critical boundary between signal and noise, i.e., the effective order. r .

[0065] To visually verify this method, test signals containing three signal components were sampled and observed by adding different white noise intensities (20dB / 40dB / 60dB) at a sampling frequency of 2000Hz. The relationship between the sequence number and amplitude obtained from the sample matrix is ​​as follows: Figure 2 As shown, when the sequence number is 6, the amplitude change gradually stabilizes, and this position can be roughly considered as the effective dividing point between the main signal and the interference signal.

[0066] An automatic order determination algorithm based on QR-SDM is used to perform column-pivoting QR decomposition on the sample matrix, utilizing... R This method automatically identifies the effective rank boundary positions of the second derivative mutation points in the diagonal principal component sequence of a matrix, without requiring manual threshold setting or reliance on computationally intensive singular value decomposition. While maintaining the computational efficiency advantage of QR decomposition, this method achieves objective determination of the model order, providing accurate and adaptive order preprocessing and online updates for the PASTD-ESPRIT algorithm. Figure 3 As shown, applying QR-SDM to the PASTD-ESPRIT algorithm enables more accurate prediction of amplitude and frequency information when facing unknown signals, effectively improving the algorithm's identification accuracy and engineering adaptability.

[0067] To verify the fast performance and algorithmic advantages of the improved PASTD-ESPRIT detection strategy described in this application, the system sampling rate was set to 2000Hz. For both the FFT and PRONY algorithms, the selection of the sliding window length needs to consider both the frequency physical resolution constraint of the FFT algorithm and the noise sensitivity suppression requirements of the PRONY algorithm. Therefore, the sliding window length parameter was set to 80 sampling points (corresponding to two fundamental cycles). The improved PASTD-ESPRIT detection algorithm described in this application, based on subspace tracking and QR rank-revealing preprocessing mechanisms, does not have a strict dependence on the sliding window length for parameter estimation accuracy. Even with the window length parameter reduced to 20 sampling points, it still maintains effective and accurate harmonic parameter identification capabilities. Specific algorithm parameters are shown in Table 2.

[0068] Table 2 Algorithm Parameter Settings

[0069] Simulation verifications were conducted under both steady-state and transient disturbance conditions. The steady-state disturbance experiment aimed to quantitatively evaluate the parameter estimation accuracy and noise robustness of the FFT algorithm, PRONY algorithm, and the improved PASTD-ESPRIT algorithm described in this application under different noise intensities. The transient disturbance experiment was used to examine the dynamic response delay and tracking convergence speed of each algorithm to changes in signal parameters after abrupt changes in harmonic components.

[0070] a) Steady-state perturbation verification To verify the effectiveness of the improved PASTD-ESPRIT detection algorithm described in this application, Gaussian white noise with intensities of 60dB, 40dB, and 20dB was injected into the experimental environment. The FFT detection algorithm, the PRONY detection algorithm, and the improved PASTD-ESPRIT detection algorithm described in this application were used to estimate the parameters of an interharmonic signal with an amplitude of 5A and a frequency of 67Hz. The estimated parameters included frequency, amplitude, and phase. The frequency, amplitude, and phase estimation curves obtained from the experiment are attached. Figure 4a , Figure 4b and Figure 4c As shown.

[0071] like Figure 4aAs shown, it should be noted that the frequency parameter estimation performance of the three algorithms differs significantly under different noise intensity conditions. Considering both the need for fast response and computational cost, this experiment uniformly uses a sliding data window length of 80 sampling points. For the sliding window FFT algorithm, its frequency physical resolution is limited by the window length parameter, corresponding to a resolution of 25Hz. Therefore, it can only output a coarse estimate of about 75Hz for the 67Hz interharmonic. In comparison, under different noise conditions, the improved PASTD-ESPRIT algorithm described in this application has better frequency parameter estimation accuracy than the PRONY algorithm. Its estimation results show lower sensitivity to changes in noise intensity, and the accuracy of parameter identification is higher.

[0072] As attached Figure 4b As shown, noise interference of different intensities has varying effects on the amplitude parameter estimation performance of the three algorithms. Due to its inherent noise sensitivity, the PRONY algorithm shows an increasing trend of deviation from the true value as the white noise intensity gradually increases. The FFT algorithm, limited by the spectral leakage effect, fails to effectively concentrate energy at interharmonic frequencies, resulting in a persistent and non-negligible systematic bias in amplitude estimation. In contrast, the improved PASTD-ESPRIT algorithm described in this application significantly outperforms the aforementioned two algorithms in amplitude parameter identification, maintaining a stable estimation result of 5A and high accuracy and stability even under varying noise intensity conditions.

[0073] As attached Figure 4c As shown, noise intensity variations have a relatively limited impact on algorithm performance in phase parameter estimation. However, the results also reveal a significant asynchrony between the FFT and PRONY algorithms in phase estimation, leading to a considerably larger phase error. Compared to the aforementioned two algorithms, the improved PASTD-ESPRIT algorithm described in this application achieves more accurate dynamic tracking of phase parameters, and its estimation results show better consistency with the true phase values.

[0074] The aforementioned performance advantages stem from the subspace projection mechanism employed by the PASTD-ESPRIT algorithm. This algorithm effectively suppresses interference energy components distributed within the noise subspace by projecting the observed signal onto the signal subspace, thereby achieving efficient separation of signal and noise components. This inherent mechanism endows the improved PASTD-ESPRIT algorithm described in this application with superior noise immunity, and its parameter estimation accuracy exhibits significantly low sensitivity to changes in the sliding data window length.

[0075] b) Transient perturbation verification Transient disturbance experiments were conducted to compare the transient response characteristics of three algorithms under a sudden increase in harmonics. It should be noted that under ideal, noise-free conditions, the PRONY algorithm exhibits a relatively fast response speed; however, when noise interference of a certain intensity exists in the system, the PRONY algorithm requires accumulating longer data samples to obtain accurate estimation results, resulting in a significant lag in its dynamic response speed. To more realistically simulate the actual operating environment under complex conditions, this experiment, under a 20dB white noise background, introduced a sudden increase of 5A and a frequency of 67Hz interharmonic components at time t1 to examine the transient tracking performance of each algorithm. The experimental results are attached. Figure 5 As shown.

[0076] From the appendix Figure 5 The experimental results show that the FFT algorithm outputs the estimation result at time t2, exhibiting the best response speed. However, the frequency value estimated by this algorithm fluctuates between 75Hz and 77Hz, showing a significant deviation from the actual interharmonic frequency of 67Hz. Considering both parameter estimation accuracy and dynamic response speed, the improved PASTD-ESPRIT algorithm described in this application demonstrates a significantly better overall transient response performance than the FFT and PRONY algorithms, effectively tracking abrupt harmonic components while ensuring fast response.

[0077] In summary, through dual experimental verification using steady-state noise disturbances and transient sudden harmonics, the improved PASTD-ESPRIT algorithm demonstrates significant comprehensive advantages in rapid harmonic detection scenarios. Compared to the FFT and PRONY algorithms, which rely on long data windows to balance resolution and noise resistance, the PASTD-ESPRIT algorithm used in this scheme requires only a short sliding window of 20 sampling points to achieve high-precision parameter estimation. Its subspace projection mechanism effectively suppresses energy interference in the noise subspace, thus maintaining extremely high accuracy and stability in tracking frequency, amplitude, and phase under different intensities of white noise. Regarding transient response, although the FFT algorithm has certain advantages in response speed, its frequency estimation results are limited by the low-resolution window function, resulting in a large inherent bias. In contrast, the improved PASTD-ESPRIT algorithm, under complex conditions with noise interference, can balance response speed and estimation accuracy with smaller data delays, achieving unbiased and rapid tracking of harmonic parameters. Therefore, this strategy effectively solves the contradiction between window length, noise immunity and fast response speed that is difficult to reconcile in traditional frequency domain or parameterization methods, and provides a more superior fast detection method for high dynamic power quality monitoring and active filter control.

[0078] Those skilled in the art will understand that, in addition to implementing the system, apparatus, and their modules provided by this invention in purely computer-readable program code, the same program can be implemented in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps. Therefore, the system, apparatus, and their modules provided by this invention can be considered a hardware component, and the modules included therein for implementing various programs can also be considered structures within the hardware component; alternatively, modules for implementing various functions can be considered both software programs implementing the method and structures within the hardware component.

[0079] Specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. Unless otherwise specified, the embodiments and features described in this application can be arbitrarily combined with each other.

Claims

1. A rapid harmonic detection system based on an improved PASTD-ESPRIT, characterized in that, It includes a QR-based rank-based preprocessing module and a PASTD-ESPRIT main detection module; The input of the QR-based preprocessing module receives the three-phase voltage or current sampling signal after coordinate transformation, and its output is connected to the order input of the PASTD-ESPRIT main detection module. The QR-based preprocessing module performs column pivoting QR decomposition on the constructed Hankel data matrix, extracts the diagonal element sequence of the upper triangular matrix R, and uses the second-order difference criterion of the diagonal elements to determine the effective order r of the signal subspace. The obtained order is then transmitted to the PASTD-ESPRIT main detection module in real time. The PASTD-ESPRIT main detection module initializes the weight matrix based on the effective order r, updates the feature vector of the signal subspace using the PASTD recursive method of the compact projection approximate subspace tracking technique, then uses the rotation-invariant technique ESPRIT to solve for the frequency of each harmonic, and completes the joint estimation of harmonic amplitude and phase through the least squares method, finally outputting the online closed-loop detection results of the three parameters of frequency, amplitude and phase.

2. The harmonic rapid detection system based on the improved PASTD-ESPRIT as described in claim 1, characterized in that, The Hankel data matrix is ​​constructed as follows: A discrete-time series is obtained by sampling a continuous signal: In the formula, It is the signal value of the nth discrete sampling point. N This represents the total number of sampling points. T s The sampling interval is... r The number of harmonics contained in the signal. e ( n () has a mean of 0 and a variance of σ 2 Gaussian white noise, A i 、f i φ i The first i The amplitude, frequency, and initial phase of the subharmonic; ω i = 2πf i , is the angular frequency of the i-th harmonic; j is the imaginary unit; Take the length as m Using a sliding window, construct a data vector: in, m > r superscript T Indicates the conjugate transpose, the vector is decomposed into signal components. X ( n ) and noise components E ( n The expression is: Y ( n )= BX ( n )+ E ( n ),in B The Vandermonde matrix is ​​composed of harmonic frequencies. X ( n It contains the complex envelope of each harmonic; the corresponding autocorrelation matrix R is: in, D = E [ X ( n ) X H ( n )], superscript H denoted as conjugate transpose, where I is the identity matrix.

3. The rapid harmonic detection system based on the improved PASTD-ESPRIT as described in claim 2, characterized in that, Column-major QR decomposition satisfies: for any m × n real matrix of order A and m ≥ n There exists a column permutation matrix Π such that the permuted matrix is... A Π is obtained through standard QR decomposition. A Π= QR Where Q is an orthogonal matrix, the upper triangular matrix R The blocks are as follows: in, R 11 ∈C r×r The rank information corresponding to the signal subspace, and its minimum singular value reflects the lower bound of the energy of the principal components of the signal; R 22 ∈C (n-r)×(n-r) Its spectral norm reflects the upper bound of the energy of the noise subspace; R 11 ∈C r×(n-r) The top right corner sub-block, 0 represents (n The zero matrix of r)×r.

4. The harmonic rapid detection system based on the improved PASTD-ESPRIT according to claim 3, characterized in that, The method for extracting the diagonal element sequence of the upper triangular matrix R is as follows: Perform QR decomposition of the Hankel data matrix with column-selected pivots, extract the main diagonal elements of matrix R, and arrange them in descending order to obtain the pivot sequence of the QR decomposition. d i =| r ii |， i =1,2…, n -1 in, r ii This represents the i-th main diagonal element of the upper triangular matrix R obtained after performing column-selective pivoting QR decomposition on the Hankel data matrix; d i It is the absolute value of the diagonal element, forming the principal component sequence; The initial order of the model; The amplitudes of the first few diagonal elements are larger values ​​within a preset range, corresponding to the main components in the signal. As the index increases, the amplitude changes abruptly at a certain position r, dropping to a smaller value within the preset range. This point of change is the critical boundary between the signal and noise.

5. The harmonic rapid detection system based on the improved PASTD-ESPRIT according to claim 4, characterized in that, diagonally... The second-order difference criterion is: through second-order difference To depict its changing trend, among which d ( k ) is the first k The size of the main diagonal elements is expressed as: along with k The value asymptotically approximates the initial order of the model. r c The absolute value of the second difference | | will continue to decrease, and when this value converges to zero, then: Corresponding k Located at the point of abrupt change in the curve, this position corresponds to the critical boundary between signal and noise, and is also the effective order. r .

6. The rapid harmonic detection system based on the improved PASTD-ESPRIT as described in claim 5, characterized in that, The process of initializing the weight matrix is ​​as follows: Set the number of iterations k =0, maximum number of iterations L = N m +1, the weight matrix is ​​initialized as follows: in, The initial values ​​for the weight matrix are m×r; Simultaneously initialize eigenvalues ​​and eigenvectors: λ i (0)=1 indicates that the initial value of the i-th eigenvalue is set to 1. w i (0)= W i (0) Initial weighting matrix of the i-th eigenvector The i-th column, i ∈[1, r ]; Iteration from k Starting with =1, first assign the current sampled data vector the value: Y(k) is a data vector constructed from the sampled signal.

7. The rapid harmonic detection system based on the improved PASTD-ESPRIT as described in claim 6, characterized in that, PASTD recursion method is for i =1 to i = r Perform the following calculations in sequence: Find the projection of the current data onto the existing feature vectors: Among them, scalar projection This indicates the magnitude of the component of the current data in the direction of the i-th feature vector; This represents the current data vector; Using the forgetting factor β Update eigenvalues, 0 < β <1. Balancing the impact of historical information and new data: in, Represents the feature value at the current moment; Update the feature vector using residuals and projections: Among them, residual This represents the portion of the current data that is not represented by the i-th feature vector; This represents the updated feature vector; Update the data vector to remove the influence of the current feature vector and prepare for the next feature vector estimation: After the iteration is complete, the updated feature vectors are arranged into a weight matrix column by column. W =[ w 1,..., w r This matrix is ​​an approximate estimate of the signal subspace.

8. The rapid harmonic detection system based on the improved PASTD-ESPRIT according to claim 7, characterized in that, The process of using rotational invariance techniques to solve for the frequencies of each harmonic is as follows: weight matrix W Perform block processing and remove W The last row yields the submatrix. W e ∈C (m 1)×r Remove W The first row yields the upper submatrix. W f ∈C (m 1)×r Constructing rotation operators S ∈C r×r It satisfies the rotation invariance relation: Solving for the rotation operator using the least squares method: Eigenvalue decomposition is performed on the rotation operator S to obtain the eigenvalues. z i , i =1,2,..., r Its mapping relationship with harmonic frequencies is as follows: The frequencies of each harmonic can be obtained by reverse calculation: In the formula, arg() represents the argument of a complex number, and its range is (- ). π , π The value is adjusted based on the harmonic frequency range of the power system.

9. The rapid harmonic detection system based on the improved PASTD-ESPRIT as described in claim 8, characterized in that, The process of jointly estimating the harmonic amplitude and phase using the least squares method is as follows: Obtaining harmonic frequencies f i Then, construct the observation matrix Φ∈C N×r : Organize the sampling sequence into an observation vector y =[ y (0), y (1),..., y ( N 1)] T Construct a least squares optimization problem: Where A=[ A 1 , A 2 ,..., A r ] T Given a complex amplitude vector, solve using canonical equations: For complex amplitude C i =A i Decompose the data to obtain the amplitude and phase: in, For the first r The initial phase of the subharmonic, Φ is the observation matrix, and Re() and Im() represent taking the real and imaginary parts, respectively.

10. The rapid harmonic detection system based on the improved PASTD-ESPRIT as described in claim 1, characterized in that, During the PASTD iteration process, the order adjustment is triggered based on the eigenvalue change pattern, and an online adaptive update mechanism for the subspace dimension is established: when the harmonic components change dynamically, the column pivot QR decomposition and diagonal second-order difference criterion are re-executed through the QR rank-revealing preprocessing module to obtain the updated effective order, and the updated order is transmitted to the PASTD-ESPRIT main detection module in real time to reinitialize the dimension of the weight matrix and the number of eigenvectors in the PASTD recursion process.

Images