CN112964931A - Non-ideal multi-damping harmonic signal parameter measurement method based on double-channel undersampling - Google Patents

Non-ideal multi-damping harmonic signal parameter measurement method based on double-channel undersampling Download PDF

Info

Publication number
CN112964931A
CN112964931A CN202110103500.8A CN202110103500A CN112964931A CN 112964931 A CN112964931 A CN 112964931A CN 202110103500 A CN202110103500 A CN 202110103500A CN 112964931 A CN112964931 A CN 112964931A
Authority
CN
China
Prior art keywords
sampling
signal
frequency
channel
parameter
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110103500.8A
Other languages
Chinese (zh)
Other versions
CN112964931B (en
Inventor
黄国兴
倪安
卢为党
彭宏
张昱
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University of Technology ZJUT
Original Assignee
Zhejiang University of Technology ZJUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University of Technology ZJUT filed Critical Zhejiang University of Technology ZJUT
Priority to CN202110103500.8A priority Critical patent/CN112964931B/en
Publication of CN112964931A publication Critical patent/CN112964931A/en
Application granted granted Critical
Publication of CN112964931B publication Critical patent/CN112964931B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

A non-ideal multi-damping harmonic signal parameter measuring method based on double-channel undersampling comprises two sampling channels, in order to prevent frequency mixing and image frequency blurring, the two sampling channels are connected through a feedback channel, signal parameters can be recovered from non-ideal K harmonic MEDS signals through only 4K samples, and the method is also suitable for K unknown conditions. Compared with the method in the prior art, the method has higher reconstruction precision and robustness to noise.

Description

Non-ideal multi-damping harmonic signal parameter measurement method based on double-channel undersampling
Technical Field
The invention belongs to the technical field of signal processing, and particularly relates to a non-ideal multi-damping harmonic signal parameter measurement method based on dual-channel undersampling.
Background
Parameter estimation of multi-Damped harmonic signals (MEDS) is a problem that often occurs in the fields of magnetic resonance imaging, linear system recognition, speech analysis, and transient analysis. A related problem includes pole-zero modeling of empirically generated time series data, since the impulse response of a linear stable rational model is composed of exponentially damped sinusoids. A number of researchers have proposed methods for parameter estimation based on nyquist sampling theory. Fotinea et al propose a state space approach to spectral estimation that performs factor 2 decimation while utilizing all available data sets. Chen et al propose a subspace-based parameter estimator that is efficient and accurate in estimating the parameters of an exponentially decaying sinusoidal signal sum. However, according to the nyquist sampling theorem, the above system needs to satisfy the requirement that the sampling rate is twice the signal bandwidth. Therefore, for wideband signals, these systems face a common problem in that not only a high sampling rate but also a complicated process is required to complete parameter estimation.
The continuous-time med signal is generally composed of a linear combination of a plurality of damped complex exponential signals, but since the actual signal cannot be perfectly matched with such a signal, the present invention takes such a model matching error into account and reduces the influence of such an error through an optimization algorithm. Also, due to the periodicity of the trigonometric function, the under-nyquist sampling results in frequency aliasing and image frequency aliasing.
In recent years, parameter measurement methods based on sub-nyquist sampling have received increasing attention. However, when measuring parameters using the under-nyquist sampling method, a key problem to be solved is frequency aliasing. Scholars of r.venkataramani and a.bourdoux, et al, propose multi-rate asynchronous sub-nyquist sampling schemes to solve the frequency aliasing problem, which typically limit the number of frequency components, also depending on the number of channels. However, these methods often require random sampling, resulting in a complex hardware structure, which greatly reduces the practicality. Furthermore, the estimation accuracy depends on the spectral grid density, which is generally not high in view of computational complexity. In the p.pal and s.qin et al systems, a two-channel sampling scheme of a relatively prime undersampling rate is employed to estimate the line spectrum. Huang et al, however, show that frequency estimates sometimes cannot be uniquely determined. To solve this problem, s.huang et al propose a three-channel under-Nyquist scheme with pairwise co-prime undersampling rate to ensure that the frequency aliasing problem can be successfully solved. However, these methods require a large number of samples to determine the correct frequency, and these systems require that the number of frequencies be known and that they are not specifically designed for the MEDS signal. Therefore, the measurement problem of the MEDS signal parameters needs to be solved urgently.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a non-ideal multi-damping harmonic signal parameter measurement method based on double-channel undersampling, aiming at the problems of frequency aliasing and image frequency aliasing and parameter measurement caused by the undersampling of non-ideal MEDS signals, the method comprises two sampling channels which are respectively a main sampling channel and an auxiliary sampling channel, and a feedback channel is arranged between the two channels, and the feedback channel can avoid the occurrence of frequency aliasing and image frequency aliasing; then, an intelligent optimization method is also provided to solve the problems of signal nonideal and insufficient prior information.
The technical scheme adopted by the invention for solving the technical problems is as follows:
a non-ideal multi-damping harmonic signal parameter measurement method based on dual-channel undersampling comprises the following steps:
firstly, initializing to generate a non-ideal MEDS signal, wherein the signal has the following form:
Figure BDA0002916479670000031
wherein T is belonged to (0, T), ck(K-1, 2, … K) is the complex amplitude, sk＝rk+j2πfkIs the complex frequency, K is 1,2, … K, rkIs the damping factor, fkThe method comprises the following steps that pole frequency is adopted, K is the number of unknown frequency components, sigma (t) is a model matching error signal, and the model matching error signal sigma (t) is introduced into a model because a real signal and an MEDS signal cannot be completely matched;
step two, controlling a signal control switch so as to shunt the MEDS signals, firstly, dialing the control switch to a branch where a main sampling channel is located through the control signal, at the moment, sampling the obtained shunt MEDS signals through the main sampling channel, then dialing the control switch to a branch where an auxiliary sampling channel is located through the control signal, and at the moment, sampling the obtained shunt MEDS signals through the auxiliary sampling channel;
step three, the control signal dials the switch to the branch where the main sampling channel is located, the main sampling channel samples the obtained shunt MEDS signal, the sampling clock CLK1 uniformly samples the shunt signal, and the sampling rate is fsThe sample values are represented as follows:
Figure BDA0002916479670000032
where K is the unknown number of components, ckWhich is indicative of the complex amplitude of the signal,
Figure BDA0002916479670000033
sk＝rk+j2πfkrepresenting complex frequency, rkRepresenting the damping factor, fkRepresenting the pole frequency, m ∈ Z+，σ[n]Represents the sample value after σ (t) is sampled by the clock CLK 1;
and step four, establishing an optimization model to reduce the influence of model matching errors. Since the number of components K is unknown and the model matching error sigma (t) ≠ 0, the optimal component K and the optimal parameter are found by minimizing the energy of sigma (t)
Figure BDA0002916479670000034
Step five, sampling value x [ n ]]Input to the optimal estimation algorithm of the number of amplitude and frequency components, thereby obtaining the optimal K and
Figure BDA0002916479670000035
p (P ≦ K) represents a valid estimate of the threshold exceeded
Figure BDA0002916479670000041
The number of (2);
step six, estimating the parameter
Figure BDA0002916479670000042
Input to a feedback sample rate generator module to calculate the sample rate f 'of the secondary sample channel'sAnd f'sSatisfies the following formula:
Figure BDA0002916479670000043
wherein,
Figure BDA0002916479670000044
a,b∈{0,1,…,K-1}，m∈{-(2Q-2),-(2Q-3),…,2Q-2}，
Figure BDA0002916479670000045
fsis the sampling rate, 0 < fs＜fmax，fmaxFor signal frequency upper limit, when the sampling value is interfered by noise, in order to improve the robustness of the method, the sampling rate f 'generated by the sampling rate generator module is fed back'sThe following equation needs to be satisfied:
Figure BDA0002916479670000046
wherein rem { (. 1) } represents a remainder of division of (-) by 1, and ε represents a threshold determined by noise intensity;
and step seven, the control signal dials the switch to the auxiliary sampling channel, and the auxiliary sampling channel samples the obtained shunt signal. The signal is uniformly sampled by a sampling clock CLK2 at a sampling rate of f'sThe sample values are represented as follows:
Figure BDA0002916479670000047
wherein, c'kIn order to be a complex amplitude of the signal,
Figure BDA0002916479670000048
T′s＝1/f′sto adoptSample period, σ]n′]＝σ[n′T′s]For model matching the sampled values of the error signal, N 'is 0,1, …, N' -1, and then an optimal estimation algorithm of the number of amplitude and frequency components is run, the sampled values x 'N']By the algorithm, an estimated value is obtained
Figure BDA0002916479670000049
Step eight, estimating the parameter
Figure BDA00029164796700000410
And
Figure BDA00029164796700000411
input to a joint estimation parameter algorithm for complex frequency and complex amplitude
Figure BDA00029164796700000412
And
Figure BDA00029164796700000413
and (4) obtaining a joint estimation.
Further, the process of the step four is as follows:
step 4.1, the established optimization model is expressed as follows:
Figure BDA0002916479670000051
due to the fact that
Figure BDA0002916479670000052
(3) To convert to:
Figure BDA0002916479670000053
step 4.2, because the input non-ideal med signal x (t) is unknown, that is, the formula (4) cannot be solved, and the sampling value obtained in step (2) is used to replace x (t), then step (4) is converted into the following formula:
Figure BDA0002916479670000054
wherein,
Figure BDA0002916479670000055
number of frequency components K and parameter
Figure BDA0002916479670000056
Are unknown, then (5) is described by the following formula:
Figure BDA0002916479670000057
wherein,
Figure BDA0002916479670000058
(7) referred to as the optimization objective function.
Still further, the process of the fifth step is as follows:
step 5.1, initialization, first, obtaining the sampling value x [ n ] from the parameter measurement system](N-0, 1, …, N-1), setting the number of frequency components K' to 0, setting the optimum number of frequency components K to 0, and setting the estimation parameter PbestSetting the optimum value of equation (7) to fbest；
Step 5.2, updating the number K' +1 of the frequency components;
step 5.3, using the sampling value x [ n ]]And solving (2) by a spectrum estimation method to obtain the estimation parameters
Figure BDA0002916479670000059
As an initial bit for the next step;
step 5.4, the optimization problem of (6) is solved by the provided optimization algorithm, and P is setoptIs a global optimum position;
step 5.5, update the optimal solution
Figure BDA00029164796700000510
Step 5.6, the result is output, if
Figure BDA0002916479670000061
And (4) outputting the solution to be the optimal solution, otherwise, returning to the step 5.3 and repeating the steps.
The process of the step eight is as follows:
step 8.1, estimated normalized frequency due to the periodicity of the trigonometric function
Figure BDA0002916479670000062
And
Figure BDA0002916479670000063
by a parameter of 2 pi mkNamely:
Figure BDA0002916479670000064
further, the following formula is obtained:
Figure BDA0002916479670000065
wherein, the angle is the main value of the complex angle, and is more than or equal to 0 and less than 2 pi, mkE.g., Z, K e.g., {1,2, …, K }. Because of the fact that
Figure BDA0002916479670000066
fmaxFor the upper frequency limit of the MEDS signal, the following is obtained:
Figure BDA0002916479670000067
note the book
Figure BDA0002916479670000068
Then (13) is expressed as the following equation:
0≤mk≤(Q-1) (14)
step 8.2, according to (12), from
Figure BDA0002916479670000069
The argument of (c) is obtained from a frequency parameter solution set F and expressed in the form:
Figure BDA00029164796700000610
wherein f issIs the sampling rate of the main sampling channel, the angle is the amplitude main value of the complex number, the angle is more than or equal to 0 and less than 2 pi, mp∈Z，-(Q-1)≤nk≤Q-1，nkE.g. Z, so the true frequency
Figure BDA00029164796700000611
Meanwhile, the estimated value of the obtained damping factor is as follows:
Figure BDA00029164796700000612
step 8.3, in the subsampling channel, according to (12), from
Figure BDA00029164796700000613
The argument of (c) is obtained as a frequency parameter solution set F' and expressed in the form:
Figure BDA0002916479670000071
wherein, f'sIs the sampling rate of the secondary sampling channel, the angle is the amplitude primary value of a complex number, and is more than or equal to 0 and less than 2 pi, m'p′∈Z，-(Q′-1)≤n′k≤Q′-1，n′kE.g. Z, so the true frequency
Figure BDA0002916479670000072
Meanwhile, the estimated value of the obtained damping factor is as follows:
Figure BDA0002916479670000073
to this end, the frequency parameter is obtained as shown in the following formula:
Figure BDA0002916479670000074
and the damping factor parameter is as follows:
Figure BDA0002916479670000075
thereby obtaining a complex frequency skAs shown in the following formula:
Figure BDA0002916479670000076
step 8.4, once the frequency is found
Figure BDA0002916479670000077
We can calculate the parameters
Figure BDA0002916479670000078
And
Figure BDA0002916479670000079
first, by selecting
Figure BDA00029164796700000710
And
Figure BDA00029164796700000711
to determine parameters
Figure BDA00029164796700000712
Assume that the other elements are
Figure BDA00029164796700000713
And
Figure BDA00029164796700000714
then, the parameters
Figure BDA00029164796700000715
And
Figure BDA00029164796700000716
respectively from the set
Figure BDA00029164796700000717
And
Figure BDA00029164796700000718
is obtained, therefore, parameter ckEstimated as:
Figure BDA00029164796700000719
the invention has the following beneficial effects: obtaining the parameter complex frequency s of the MEDS signal to be measured by using the optimal estimation algorithm and the combined estimation parameter algorithm of the amplitude and frequency component numberkAnd complex amplitude ckA value of (d); the method has high reconstruction precision and robustness to noise.
Drawings
FIG. 1 is a system structure diagram of a non-ideal multi-damping harmonic signal parameter measurement method based on two-channel undersampling.
Fig. 2 is a flow chart of an optimal estimation algorithm for the number of amplitude and frequency components.
Fig. 3 is a result of parameter recovery when the frequency is not aliased.
Fig. 4 is a parameter recovery result at the time of frequency aliasing.
Figure 5 is a comparison of parameter recovery performance as signal-to-noise ratio increases.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
Referring to fig. 1 to 5, a method for measuring non-ideal multi-damping harmonic signal parameters based on dual-channel undersampling comprises the following steps:
firstly, initializing to generate a non-ideal MEDS signal, wherein the signal has the following form:
Figure BDA0002916479670000081
wherein T is belonged to (0, T), ck(K-1, 2, … K) is the complex amplitude, sk＝rk+j2πfkIs the complex frequency, K is 1,2, … K, rkIs the damping factor, fkIs the pole frequency, K is the unknown number of frequency components, and σ (t) is the model match error signal. Since the true signal and the med signal cannot be completely matched, we introduce a model matching error signal σ (t) in the model;
step two, controlling a signal control switch so as to shunt the MEDS signals, firstly, dialing the control switch to a branch where a main sampling channel is located through the control signal, at the moment, sampling the obtained shunt MEDS signals through the main sampling channel, then dialing the control switch to a branch where an auxiliary sampling channel is located through the control signal, and at the moment, sampling the obtained shunt MEDS signals through the auxiliary sampling channel;
step three, the control signal dials the switch to the branch where the main sampling channel is located, the main sampling channel samples the obtained shunt MEDS signal, the sample clock CLK1 uniformly samples the shunt signal, and the sampling rate is fsThe sample values are represented as follows:
Figure BDA0002916479670000091
where K is the unknown number of components, ckWhich is indicative of the complex amplitude of the signal,
Figure BDA0002916479670000092
sk＝rk+j2πfkrepresenting complex frequency, rkRepresenting the damping factor, fkRepresenting the pole frequency, m ∈ Z+，σ[n]Represents the sample value after σ (t) is sampled by the clock CLK 1;
step four, establishing an optimization model, reducing the influence of model matching errors, and finding the optimal component K and the optimal parameter by minimizing the energy of sigma (t) because the component number K is unknown and the model matching error sigma (t) is not equal to 0
Figure BDA0002916479670000093
The process is as follows:
step 4.1, the established optimization model is expressed as follows:
Figure BDA0002916479670000094
due to the fact that
Figure BDA0002916479670000095
(3) To convert to:
Figure BDA0002916479670000096
step 4.2, because the input non-ideal med signal x (t) is unknown, that is, the formula (4) cannot be solved, and the sampling value obtained in step (2) is used to replace x (t), then step (4) is converted into the following formula:
Figure BDA0002916479670000097
wherein,
Figure BDA0002916479670000098
number of frequency components K and parameter
Figure BDA0002916479670000099
Are unknown, then (5) is described by the following formula:
Figure BDA00029164796700000910
wherein,
Figure BDA00029164796700000911
(7) called the optimization objective function;
step five, sampling value x [ n ]]Input to the optimal estimation algorithm of the number of amplitude and frequency components, thereby obtaining the optimal K and
Figure BDA0002916479670000101
p (P ≦ K) represents a valid estimate of the threshold exceeded
Figure BDA0002916479670000102
The number of (2) is as follows:
step 5.1, initialization, first, the sampling values x [ n ] can be obtained from the parameter measurement system shown in FIG. 1](N-0, 1, …, N-1), setting the number of frequency components K' to 0, setting the optimum number of frequency components K to 0, and setting the estimation parameter PbestSetting the optimum value of equation (7) to fbest；
Step 5.2, updating the number K' +1 of the frequency components;
step 5.3, using the sampling value x [ n ]]And solving (2) by a spectrum estimation method to obtain the estimation parameters
Figure BDA0002916479670000103
As an initial bit for the next step;
step 5.4, the optimization problem of (6) is solved by the proposed optimization algorithm, and we assume PoptIs a global optimum position;
step 5.5, update the optimal solution
Figure BDA0002916479670000104
Step 5.6, the result is output, if
Figure BDA0002916479670000105
The output is the optimal solution, otherwise, the step 5.3 is returned to, and the steps are repeated;
step (ii) ofSixthly, estimating the parameters
Figure BDA0002916479670000106
Input to a feedback sample rate generator module that calculates the sample rate f 'of the secondary sampling channel'sAnd f'sSatisfies the following formula:
Figure BDA0002916479670000107
wherein,
Figure BDA0002916479670000108
a,b∈{0,1,…,K-1}，m∈{-(2Q-2),-(2Q-3),…,2Q-2}，
Figure BDA0002916479670000109
fsis the sampling rate, 0 < fs＜fmax，fmaxFor the upper limit of signal frequency, when the sampling value is interfered by noise, in order to improve the robustness of the method, the sampling rate f 'generated by the sampling rate generator module is fed back'sThe following equation needs to be satisfied:
Figure BDA00029164796700001010
wherein rem { (. 1) } represents a remainder of division of (-) by 1, and ε represents a threshold determined by noise intensity;
step seven, the switch is set to be in a sub-sampling channel by the control signal, the sub-sampling channel samples the obtained shunt signal, the sampling clock CLK2 samples the shunt signal uniformly, and the sampling rate is f'sThe sample values are represented as follows:
Figure BDA0002916479670000111
wherein, c'kIn order to be a complex amplitude of the signal,
Figure BDA0002916479670000112
T′s＝1/f′sis the sampling period, σ [ n']＝σ[n′T′s]For model matching the sampled values of the error signal, N 'is 0,1, …, N' -1, and then an optimal estimation algorithm of the number of amplitude and frequency components is run, the sampled values x 'N']Through the algorithm, an estimated value can be obtained
Figure BDA0002916479670000113
Step eight, estimating the parameter
Figure BDA0002916479670000114
And
Figure BDA0002916479670000115
input into the joint estimation parameter algorithm, the complex frequency and the complex amplitude can be used
Figure BDA0002916479670000116
And
Figure BDA0002916479670000117
and (4) obtaining a joint estimation. The algorithm process is as follows:
step 8.1, estimated normalized frequency due to the periodicity of the trigonometric function
Figure BDA0002916479670000118
And
Figure BDA0002916479670000119
by a parameter of 2 pi mkNamely:
Figure BDA00029164796700001110
further, the following formula is obtained:
Figure BDA00029164796700001111
wherein, the angle is the main value of the complex angle, and is more than or equal to 0 and less than 2 pi, mkE.g., Z, K e.g., {1,2, …, K }. Because of the fact that
Figure BDA00029164796700001112
fmaxFor the upper frequency limit of the MEDS signal, the following is obtained:
Figure BDA00029164796700001113
note the book
Figure BDA0002916479670000121
Then (13) is expressed as the following equation:
0≤mk≤(Q-1) (14)
step 8.2, according to (12), from
Figure BDA0002916479670000122
The argument of (c) is obtained from a frequency parameter solution set F and expressed in the form:
Figure BDA0002916479670000123
wherein f issIs the sampling rate of the main sampling channel, the angle is the amplitude main value of the complex number, the angle is more than or equal to 0 and less than 2 pi, mp∈Z，-(Q-1)≤nk≤Q-1，nkE.g. Z, so the true frequency
Figure BDA0002916479670000124
Meanwhile, the estimated value of the obtained damping factor is as follows:
Figure BDA0002916479670000125
step 8.3, in the subsampling channel, according to (12), from
Figure BDA0002916479670000126
The argument of (c) is obtained as a frequency parameter solution set F' and expressed in the form:
Figure BDA0002916479670000127
wherein, f'sIs the sampling rate of the secondary sampling channel, the angle is the amplitude primary value of a complex number, and is more than or equal to 0 and less than 2 pi, m'p′∈Z，-(Q′-1)≤n′k≤Q′-1，n′kE.g. Z, so the true frequency
Figure BDA0002916479670000128
Meanwhile, the estimated value of the obtained damping factor is as follows:
Figure BDA0002916479670000129
to this end, the frequency parameter is obtained as shown in the following formula:
Figure BDA00029164796700001210
and the damping factor parameter is as follows:
Figure BDA00029164796700001211
thereby obtaining a complex frequency skAs shown in the following formula:
Figure BDA0002916479670000131
step 8.4, once the frequency is found
Figure BDA0002916479670000132
Calculating parameters
Figure BDA0002916479670000133
And
Figure BDA0002916479670000134
first, by selecting
Figure BDA0002916479670000135
And
Figure BDA0002916479670000136
to determine parameters
Figure BDA0002916479670000137
Let other elements be
Figure BDA0002916479670000138
And
Figure BDA0002916479670000139
then, the parameters
Figure BDA00029164796700001310
And
Figure BDA00029164796700001311
respectively from the set
Figure BDA00029164796700001312
And
Figure BDA00029164796700001313
is obtained, therefore, parameter ckEstimated as:
Figure BDA00029164796700001314
example (c): in a first experiment, it was verified that the method of the present invention can reconstruct the original signal with a small number of sample values under a noise-free condition, and compare it with a clock-interleaved sampling system. In this experiment, the number of frequency components K is set to 5, and the frequency f is set tok＝[0.4,3,5.2,6.8,9.1]KHz。The method of the invention can recover the parameters of the original signal by using 4K to 20 sampling values, and the clock staggered sampling system can also accurately recover the parameters of the original signal, and the recovery result is shown in figure 3.
To further verify the validity of the method of the invention, we set the frequency fk＝[0.4,3,5.5,6.8,9.1]Due to the fact that
Figure BDA00029164796700001315
This will cause aliasing at 3KHz and 5.5KHz, and as can be seen from fig. 4, the feedback sampling system provided by the present invention can accurately recover the parameters of the original signal, whereas the clock interleaving sampling system cannot accurately recover the parameters of the original signal.
In a second experiment, we compared the method of the present invention with several other methods in the presence of noise, including in particular, a clocked interleaved sampling system, a time delayed sampling system, and a three channel co-prime sampling system. In the experiment process, we set the number of frequency components K to 5, and the frequency fkRandomly distributing the signals on (0,10) KHz, and then adding Gaussian noise to the MEDS signals to be measured, wherein the signal-to-noise ratio SNR of the noise is defined as follows:
Figure BDA00029164796700001316
wherein, PsignalAnd PnoiseRespectively representing the power of the signal and the noise, and the SNR ranges from [ -20,50 [ -20 [)]dB, the number of samples Num obtained by all systems is 300. The experiment is carried out for 2000 times totally, the experimental result is shown in fig. 5, and it can be seen from fig. 5 that the performance difference between the clock staggered sampling system and the time delay sampling system is not large, when the SNR is less than or equal to 25dB, the performance of the method of the invention is similar to that of the three-channel co-prime sampling system, but with the increase of the SNR, the method of the invention has higher accuracy than that of other systems. Experimental results show that the method has better robustness.
The embodiments described in this specification are merely illustrative of implementations of the inventive concepts, which are intended for purposes of illustration only. The scope of the present invention should not be construed as being limited to the particular forms set forth in the examples, but rather as being defined by the claims and the equivalents thereof which can occur to those skilled in the art upon consideration of the present inventive concept.

Claims (4)

1. A non-ideal multi-damping harmonic signal parameter measurement method based on double-channel undersampling is characterized by comprising the following steps:
firstly, initializing to generate a non-ideal MEDS signal, wherein the signal has the following form:
Figure FDA0002916479660000011
wherein T is belonged to (0, T), ck(K-1, 2, … K) is the complex amplitude, sk＝rk+j2πfkIs the complex frequency, K is 1,2, … K, rkIs the damping factor, fkThe method comprises the following steps that pole frequency is adopted, K is the number of unknown frequency components, sigma (t) is a model matching error signal, and the model matching error signal sigma (t) is introduced into a model because a real signal and an MEDS signal cannot be completely matched;
step two, controlling a signal control switch so as to shunt the MEDS signals, firstly, dialing the control switch to a branch where a main sampling channel is located through the control signal, at the moment, sampling the obtained shunt MEDS signals through the main sampling channel, then dialing the control switch to a branch where an auxiliary sampling channel is located through the control signal, and at the moment, sampling the obtained shunt MEDS signals through the auxiliary sampling channel;
step three, the control signal dials the switch to the branch where the main sampling channel is located, the main sampling channel samples the obtained shunt MEDS signal, the sampling clock CLK1 uniformly samples the shunt signal, and the sampling rate is fsThe sample values are represented as follows:
Figure FDA0002916479660000012
where K is the unknown number of components, ckWhich is indicative of the complex amplitude of the signal,
Figure FDA0002916479660000013
sk＝rk+j2πfkrepresenting complex frequency, rkRepresenting the damping factor, fkRepresenting the pole frequency, m ∈ Z+，σ[n]Represents the sample value after σ (t) is sampled by the clock CLK 1;
step four, establishing an optimization model, reducing the influence of model matching errors, and finding the optimal component K and the optimal parameter by minimizing the energy of sigma (t) because the component number K is unknown and the model matching error sigma (t) is not equal to 0
Figure FDA0002916479660000014
Step five, sampling value x [ n ]]Input to the optimal estimation algorithm of the number of amplitude and frequency components, thereby obtaining the optimal K and
Figure FDA0002916479660000015
p (P ≦ K) represents a valid estimate of the threshold exceeded
Figure FDA0002916479660000016
The number of (2);
step six, estimating the parameter
Figure FDA0002916479660000017
Input to a feedback sample rate generator module to calculate the sample rate f 'of the secondary sample channel'sAnd f'sSatisfies the following formula:
Figure FDA0002916479660000018
wherein,
Figure FDA0002916479660000019
a,b∈{0,1,…,K-1}，m∈{-(2Q-2),-(2Q-3),…,2Q-2}，
Figure FDA0002916479660000021
fsis the sampling rate, 0 < fs＜fmax，fmaxFor signal frequency upper limit, when the sampling value is interfered by noise, in order to improve the robustness of the method, the sampling rate f 'generated by the sampling rate generator module is fed back'sThe following equation needs to be satisfied:
Figure FDA0002916479660000022
wherein rem { (. 1) } represents a remainder of division of (-) by 1, and ε represents a threshold determined by noise intensity;
step seven, the switch is set to be in a sub-sampling channel by the control signal, the sub-sampling channel samples the obtained shunt signal, the sampling clock CLK2 samples the shunt signal uniformly, and the sampling rate is f'sThe sample values are represented as follows:
Figure FDA0002916479660000023
wherein, c'kIn order to be a complex amplitude of the signal,
Figure FDA0002916479660000024
T′s＝1/f′sis the sampling period, σ [ n']＝σ[n′T′s]For model matching the sampled values of the error signal, N 'is 0,1, …, N' -1, and then an optimal estimation algorithm of the number of amplitude and frequency components is run, the sampled values x 'N']By the algorithm, an estimated value is obtained
Figure FDA0002916479660000025
Step eight, estimating the parameter
Figure FDA0002916479660000026
And
Figure FDA0002916479660000027
input to a joint estimation parameter algorithm for complex frequency and complex amplitude
Figure FDA0002916479660000028
And
Figure FDA0002916479660000029
and (4) obtaining a joint estimation.
2. The method for measuring the parameters of the non-ideal multi-damping harmonic signal based on the two-channel undersampling as claimed in claim 1, wherein the process of the fourth step is as follows:
step 4.1, the established optimization model is expressed as follows:
Figure FDA00029164796600000210
due to the fact that
Figure FDA00029164796600000211
(3) To convert to:
Figure FDA00029164796600000212
step 4.2, because the input non-ideal med signal x (t) is unknown, that is, the formula (4) cannot be solved, and the sampling value obtained in step (2) is used to replace x (t), then step (4) is converted into the following formula:
Figure FDA00029164796600000213
wherein,
Figure FDA00029164796600000214
number of frequency components K and parameter
Figure FDA00029164796600000215
Are unknown, then (5) is described by the following formula:
Figure FDA00029164796600000216
wherein,
Figure FDA0002916479660000031
(7) referred to as the optimization objective function.
3. The method for measuring the parameters of the non-ideal multi-damping harmonic signal based on the two-channel undersampling as claimed in claim 1 or 2, characterized in that the process of the step five is as follows:
step 5.1, initialization, first, obtaining the sampling value x [ n ] from the parameter measurement system](N-0, 1, …, N-1), setting the number of frequency components K' to 0, setting the optimum number of frequency components K to 0, and setting the estimation parameter PbestSetting the optimum value of equation (7) to fbest；
Step 5.2, updating the number K' +1 of the frequency components;
step 5.3, using the sampling value x [ n ]]And solving (2) by a spectrum estimation method to obtain the estimation parameters
Figure FDA0002916479660000032
As an initial bit for the next step;
step 5.4, the optimization problem of (6) is solved by the provided optimization algorithm, and P is setoptIs a global optimum position;
step 5.5, update the optimal solution
Figure FDA0002916479660000033
Step 5.6, the result is output, if
Figure FDA0002916479660000034
And (4) outputting the solution to be the optimal solution, otherwise, returning to the step 5.3 and repeating the steps.
4. The method for measuring the parameters of the non-ideal multi-damping harmonic signal based on the two-channel undersampling as claimed in claim 1 or 2, characterized in that the procedure of the step eight is as follows:
step 8.1, estimated normalized frequency due to the periodicity of the trigonometric function
Figure FDA0002916479660000035
And
Figure FDA0002916479660000036
by a parameter of 2 pi mkNamely:
Figure FDA0002916479660000037
further, the following formula is obtained:
Figure FDA0002916479660000038
wherein, the angle is the main value of the complex angle, and is more than or equal to 0 and less than 2 pi, mkE.g., Z, K e {1,2, …, K }, since
Figure FDA0002916479660000039
fmaxFor the upper frequency limit of the MEDS signal, the following is obtained:
Figure FDA00029164796600000310
note the book
Figure FDA00029164796600000311
Then (13) is expressed as the following equation:
0≤mk≤(Q-1) (14)
step 8.2, according to (12), from
Figure FDA00029164796600000312
The argument of (c) is obtained from a frequency parameter solution set F and expressed in the form:
Figure FDA0002916479660000041
wherein f issIs the sampling rate of the main sampling channel, the angle is the amplitude main value of the complex number, the angle is more than or equal to 0 and less than 2 pi, mp∈Z，-(Q-1)≤nk≤Q-1，nkE.g. Z, so the true frequency
Figure FDA0002916479660000042
Meanwhile, the estimated value of the obtained damping factor is as follows:
Figure FDA0002916479660000043
step 8.3, in the subsampling channel, according to (12), from
Figure FDA0002916479660000044
The argument of (c) is obtained as a frequency parameter solution set F' and expressed in the form:
Figure FDA0002916479660000045
wherein f iss'is the sampling rate of the secondary sampling channel, and the angle is the amplitude primary value of a complex number, and is more than or equal to 0 and less than 2 pi, m'p′∈Z，-(Q′-1)≤n′k≤Q′-1，n′kE.g. Z, so the true frequency
Figure FDA0002916479660000046
Meanwhile, the estimated value of the obtained damping factor is as follows:
Figure FDA0002916479660000047
to this end, the frequency parameter is obtained as shown in the following formula:
Figure FDA0002916479660000048
and the damping factor parameter is as follows:
Figure FDA0002916479660000049
thereby obtaining a complex frequency skAs shown in the following formula:
Figure FDA00029164796600000410
step 8.4, once the frequency is found
Figure FDA00029164796600000411
Calculating parameters
Figure FDA00029164796600000412
And
Figure FDA00029164796600000413
first, by selecting
Figure FDA00029164796600000414
And
Figure FDA00029164796600000415
to determine parameters
Figure FDA00029164796600000416
Assume that the other elements are
Figure FDA00029164796600000417
And
Figure FDA00029164796600000418
then, the parameters
Figure FDA00029164796600000419
And
Figure FDA00029164796600000420
respectively from the set
Figure FDA00029164796600000421
And
Figure FDA00029164796600000422
is obtained, therefore, parameter ckEstimated as:
Figure FDA00029164796600000423
CN202110103500.8A 2021-01-26 2021-01-26 Non-ideal multi-damping harmonic signal parameter measurement method based on two-channel undersampling Active CN112964931B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110103500.8A CN112964931B (en) 2021-01-26 2021-01-26 Non-ideal multi-damping harmonic signal parameter measurement method based on two-channel undersampling

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110103500.8A CN112964931B (en) 2021-01-26 2021-01-26 Non-ideal multi-damping harmonic signal parameter measurement method based on two-channel undersampling

Publications (2)

Publication Number Publication Date
CN112964931A true CN112964931A (en) 2021-06-15
CN112964931B CN112964931B (en) 2022-07-15

Family

ID=76272648

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110103500.8A Active CN112964931B (en) 2021-01-26 2021-01-26 Non-ideal multi-damping harmonic signal parameter measurement method based on two-channel undersampling

Country Status (1)

Country Link
CN (1) CN112964931B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115276799A (en) * 2022-07-27 2022-11-01 西安理工大学 Decision threshold self-adapting method for undersampling modulation and demodulation in optical imaging communication
CN115436702A (en) * 2022-09-02 2022-12-06 浙江工业大学 Non-ideal multi-damping harmonic signal multi-channel undersampling method

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101713795A (en) * 2009-09-09 2010-05-26 中国科学院国家授时中心 Method of digitalized measuring frequency in restriction of sampling rate
KR20150047109A (en) * 2013-10-23 2015-05-04 삼성전자주식회사 Magnetic resonance imaging apparatus and imaging method for magnetic resonance image thereof
EP3301461A1 (en) * 2016-09-28 2018-04-04 Siemens Aktiengesellschaft Method for detection of harmonics of a univariate signal
CN111224672A (en) * 2020-01-16 2020-06-02 哈尔滨工业大学 Multi-harmonic signal undersampling method based on multi-channel time delay

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101713795A (en) * 2009-09-09 2010-05-26 中国科学院国家授时中心 Method of digitalized measuring frequency in restriction of sampling rate
KR20150047109A (en) * 2013-10-23 2015-05-04 삼성전자주식회사 Magnetic resonance imaging apparatus and imaging method for magnetic resonance image thereof
EP3301461A1 (en) * 2016-09-28 2018-04-04 Siemens Aktiengesellschaft Method for detection of harmonics of a univariate signal
CN111224672A (en) * 2020-01-16 2020-06-02 哈尔滨工业大学 Multi-harmonic signal undersampling method based on multi-channel time delay

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
黄山等: "三通道欠采样频率估计", 《华中科技大学学报(自然科学版)》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115276799A (en) * 2022-07-27 2022-11-01 西安理工大学 Decision threshold self-adapting method for undersampling modulation and demodulation in optical imaging communication
CN115436702A (en) * 2022-09-02 2022-12-06 浙江工业大学 Non-ideal multi-damping harmonic signal multi-channel undersampling method

Also Published As

Publication number Publication date
CN112964931B (en) 2022-07-15

Similar Documents

Publication Publication Date Title
CN112964931B (en) Non-ideal multi-damping harmonic signal parameter measurement method based on two-channel undersampling
CN107102255B (en) Single ADC acquisition channel dynamic characteristic test method
CN111224672B (en) Multi-channel delay-based multi-harmonic signal undersampling method
CN109407501B (en) Time interval measuring method based on relevant signal processing
CN103941088A (en) Method for quickly measuring frequency of electric power system based on three-phase signals
CN106546817B (en) A kind of Frequency Estimation and energy state postulate with feedback function
Zou et al. Blind timing skew estimation using source spectrum sparsity in time-interleaved ADCs
CN112162152B (en) Sine wave coherent pulse train signal frequency estimation method based on phase straight line fitting
CN112859019A (en) Intra-pulse modulation type parameter extraction system and using method
Fu et al. Parameter Measurement of $ M $-Ary PSK Signals With Finite Rate of Innovation
Chen et al. Robust precise time difference estimation based on digital zero-crossing detection algorithm
CN116522269B (en) Fault diagnosis method based on Lp norm non-stationary signal sparse reconstruction
CN112595889B (en) under-Nyquist sampling and parameter measuring method for non-ideal multi-exponential decay sinusoidal signal
Huang et al. Sub-Nyquist sampling of multiple exponentially damped sinusoids with feedback structure
CN110808929A (en) Real-complex conversion type signal-to-noise ratio estimation algorithm of subtraction strategy
CN107870338B (en) A kind of satellite navigation carrier wave tracing method of low update frequency
Yun et al. Sub-Nyquist sampling and measurement of MPSK signal based on parameter matching
CN112883787B (en) Short sample low-frequency sinusoidal signal parameter estimation method based on spectrum matching
Kusuma et al. On the accuracy and resolution of powersum-based sampling methods
CN112129983B (en) Waveform recovery data processing method based on equivalent sampling at equal time intervals
CN103746699A (en) Signal reconstruction method based on rotation matrix error estimation for alternative sampling system
Tamim et al. Hilbert transform of FFT pruned cross correlation function for optimization in time delay estimation
CN107576842B (en) Broadband synchronous sampling method
CN112953468A (en) Multi-exponential decay sinusoidal signal feedback type under-sampling hardware implementation method
Wu et al. A faster method for accurate spectral testing without requiring coherent sampling

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant