CN106154035A - A kind of quickly harmonic wave and harmonic detection method - Google Patents

Abstract

Power network signal is quickly analyzed by the present invention first with adding five maximum side lobe attenuation window FFT interpolation algorithms, each harmonic wave of preliminary acquisition and m-Acetyl chlorophosphonazo frequency content；And then utilize the advantage that recurrence band filter single-point amount of calculation is little, the frequency range of measured harmonic wave and m-Acetyl chlorophosphonazo is analyzed, thus contained harmonic wave and the accurate amplitude-frequency characteristic of m-Acetyl chlorophosphonazo in trying to achieve this frequency range.So, the traversal value of recurrence band filter mid frequency is only limitted to this frequency band range, it is to avoid analysis to the frequency range without harmonic wave or m-Acetyl chlorophosphonazo.Greatly reduce unnecessary operand, ensure precision and real-time that algorithm realizes in DSP platform.

Description

A kind of quickly harmonic wave and harmonic detection method

Technical field

The present invention relates to the power network signal detection technique field containing harmonic wave Yu m-Acetyl chlorophosphonazo.

Background technology

Along with commercial scale constantly expands, containing increasing nonlinear-load in power system, just cause in electrical network String voltage is distorted, and has a strong impact on the quality of power supply, brings harm to commercial production.Sinusoidal voltage waveform is distorted and means Containing substantial amounts of harmonic wave (frequency is first-harmonic integral multiple) and m-Acetyl chlorophosphonazo (frequency is that first-harmonic is non-integral multiple).Harmonic wave and the inspection of m-Acetyl chlorophosphonazo Survey is harmonic wave control, the premise evaluating the quality of power supply.FFT Harmonic detection is that the application of current electric harmonic detection field is most A kind of method, can accurately detect out harmonic components when synchronized sampling.But FFT when non-synchronous sampling because there is spectrum leakage And fence effect, there is bigger error in the detection to harmonic wave and m-Acetyl chlorophosphonazo.Adding five maximum side lobe attenuation window FFT interpolation algorithms can Reduce the spectrum leakage and fence effect caused by non-synchronous sampling, under non-synchronous sampling, also can effectively detect each harmonic Point, but slightly worse to the accuracy of detection of m-Acetyl chlorophosphonazo.The recurrence band filter utilizing mid frequency variable carries out traversal filter to signal Ripple, can accurately detect harmonic wave and m-Acetyl chlorophosphonazo composition, and resolution is higher.But when the signal wider by its detection bandwidth, in Frequency of heart value is too much, and traversal calculates tediously long, and algorithm is complicated.Site Detection instrument uses DSP, complicated algorithm can make reality mostly Time property is deteriorated, and current the method is difficult in direct equipment the most at the scene apply.Therefore, for use DSP field apparatus, need into Harmonic wave and the m-Acetyl chlorophosphonazo detection algorithm that precision is taken into account is explored and studied to one step with real-time.

Summary of the invention

It is an object of the invention to study a kind of precision and the harmonic wave of real-time and m-Acetyl chlorophosphonazo detection algorithm taken into account, propose five Maximum side lobe attenuation window interpolation FFT combines with recurrence band filter algorithm, thus realize harmonic wave contained in electrical network and The quick detection of m-Acetyl chlorophosphonazo signal, this method is for realizing harmonic wave contained in electrical network and m-Acetyl chlorophosphonazo signal in DSP platform Accurately analyze and be significant.

For reaching above-mentioned purpose, the present invention by the following technical solutions:

A kind of quickly harmonic wave and harmonic detection method, comprise the following steps:

S101: sample to measured signal under conditions of meeting sampling thheorem, it is thus achieved that its discrete sample signals x (n), Wherein n is dis-crete sample values；

S102: described discrete sample signals is added five maximum side lobe attenuation windows and is analyzed, to the windowing signal obtained x_wN () utilizes fft analysis to obtain its discrete spectrum X (f), utilize interpolation algorithm to be modified described discrete spectrum X (f), Substantially radio-frequency component scope to harmonic wave and m-Acetyl chlorophosphonazo；

S103: set the mid frequency that performance metrics factor a of recurrence band filter, sampling period Ts and traversal calculate ω_nValue；According to harmonic wave to be obtained and m-Acetyl chlorophosphonazo Detection results, determine the bandwidth of recurrence band filter, and then setting is passed Return the value of performance metrics factor a of band filter；

S104: according to scale factor a, sampling period Ts and the mid frequency ω that set_nValue determine at different centers The coefficient δ of recurrence band filter at frequency_iAnd λ_jValue and the amplitude-versus-frequency curve result of recurrence band filter preserving, Wherein i=1,2,3,4,5,6, j=1,2,3,4,5,6,7；

S105: apply recurrence band filter algorithm that sampled value x (n) of measured signal is counted at each mid frequency Calculating, sets the seventh day of lunar month the output valve of recurrence band filter as zero, utilization is in the first six sampled value of this sampled point and filtering The first seven output valve of device calculates its output result at the recurrence band filter of this sampled point, calculates the real part of output result With imaginary part, obtain in instantaneous amplitude testing result A (n) of this sampled point and frequency detecting result f (n)；

S106: judge whether recurrence band filter completes the calculating to all sampled points；If being not fully complete, then proceed to step Rapid S105 continues to calculate；If completing, carry out step S107；

S107: judge amplitude and rate-adaptive pacemaker result stable moment, takes amplitude and reaches the instantaneous of a period of time after stable state The meansigma methods of frequency and instantaneous amplitude output result is as the testing result at this mid frequency；

S108: ask for the correction coefficient under set performance metrics factor a, is carried out required amplitude detection result Preliminary corrections；Whether determination frequency testing result deviates the mid frequency of setting, if deviation, according to the frequency detected, calculates This frequency and the difference of its place mid-band frequency, obtain amplitude accurately according to wave filter amplitude-versus-frequency curve (accompanying drawing 2) Correction coefficient, corrects further to amplitude detection result, it is thus achieved that the final detection result of amplitude；

S109: judged whether that the traversal to whole mid frequencyes calculates；If being not fully complete, then go to step S103 and again set The value of the frequency that centers, continues to calculate.

The quick harmonic wave proposed according to the present invention and harmonic detection method, in described step S102, five maximum secondary lobes decline The discrete time-domain expression formula subtracting window is:

$w\left(n\right)=\underset{1=0}{\overset{4}{\Σ}}{(-1)}^1{b}_{1}cos\left(\frac{2\mathrm{\π}nl}{N}\right)$

Wherein l=0,1,2,3,4, N is maximum number of samples, b_lBeing the coefficient of five maximum side lobe attenuation windows, value is respectively For: b0=0.2734375, b1=0.4375, b2=0.21875, b3=0.0625, b4=0.0078125；

Described windowing signal x_w(n)=x (n) * w (n), the expression formula of described discrete spectrum X (f) is:

$X\left(f\right)=\frac{A_k}{2}\underset{i=1}{\overset{n}{\Σ}}\[W(f-f_k)e^{-j\left(\mathrm{\π}\right(f-f_k)-{\mathrm{\φ}}_k)}+W(f+f_k)e^{-j\left(\mathrm{\π}\right(f+f_k)+{\mathrm{\φ}}_k)}\]$

A in formula_kFor the amplitude of kth subharmonic, ψ_kPhase place for kth subharmonic；W (f) is five maximum side lobe attenuation window letters The frequency spectrum of number；f_kFor kth subfrequency by the value after frequency resolution normalization, its expression formula is f_k=kf₀×N/f_s, f₀For Fundamental frequency.

The quick harmonic wave proposed according to the present invention and harmonic detection method, wherein, utilize described discrete spectrum X (f) The step that interpolation algorithm is modified includes:

The spectral line of described discrete spectrum X (f) is spaced apart Δ f=fs/N, and wherein fs is sample frequency；If treating measured frequency f_m= k_m* Δ f, wherein k_mFor coefficient of frequency, it it is not the most integer；If f_mNeighbouring peak value spectral line is k, and corresponding amplitude is X (k), draws Enter frequency departure q=k_m-k-0.5, then have-0.5≤q≤0.5；

If the amplitude corresponding to spectral line of the peak value spectral line left and right sides is respectively X (k-1), X (k+1), if there being X (k-1)≤X (k+1), i.e.-1 spectral line of kth is time maximum spectral line near crest frequency, and making it is q with the ratio of peak value spectral line₁, then this Time:

$q_1=\frac{\left|X(k-1)\right|}{\left|X\left(k\right)\right|}=\frac{\left|W(q-1)\right|}{\left|W\left(q\right)\right|}$

Can obtain the expression formula of frequency departure q:

$q=\frac{5q_1-4}{1+q_1}$

If having X (k+1)≤X (k-1), i.e.+1 spectral line of kth is time maximum spectral line near crest frequency, make itself and peak value The ratio of spectral line is q₂, then have

$q_2=\frac{\left|X(k+1)\right|}{\left|X\left(k\right)\right|}=\frac{\left|W(q+1)\right|}{\left|W\left(q\right)\right|}$

Can obtain the expression formula of frequency departure q:

$q=-\frac{5q_2-4}{1+q_2}$

Then its frequency correction formula is: f=(k-q) * fs/N.

The quick harmonic wave proposed according to the present invention and harmonic detection method, wherein, described step S105 includes:

Carry out traversal filtering, at each mid frequency, apply recurrence band filter algorithm centrifugal pump x to measured signal N () calculates, if instantaneous output result is w after the filtered device of discrete spectrum X (f)_η(a,ω_n, n), its expression formula is as follows:

$w_{\mathrm{\η}}(a,{\mathrm{\ω}}_n,n)=T_s\underset{i=1}{\overset{6}{\Σ}}{\mathrm{\δ}}_{i}x(n-i)-\underset{j=1}{\overset{7}{\Σ}}{\mathrm{\λ}}_j{w}_{\mathrm{\η}}(a,{\mathrm{\ω}}_n,n-j)$

Wherein δ_iWith λ_jFor constant, by performance metrics factor a given and ω_nTry to achieve, particularly as follows:

${\mathrm{\δ}}_1=\[\frac{1}{12}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4-\frac{1}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{1}{90}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A/\sqrt{a};$

${\mathrm{\δ}}_2=\[\frac{3}{4}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4-\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{19}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^2/\sqrt{a};$

${\mathrm{\δ}}_3=\[-\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4-{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{151}{45}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^3/\sqrt{a};$

${\mathrm{\δ}}_4=\[-\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4+{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{151}{45}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^4/\sqrt{a};$

${\mathrm{\δ}}_5=\[\frac{3}{4}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4+\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{19}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^5/\sqrt{a};$

${\mathrm{\δ}}_6=\[\frac{1}{12}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4+\frac{1}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{1}{90}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^6/\sqrt{a};$

λ₁=-7A, λ₂=21A², λ₃=-35A³, λ₄=35A⁴, λ₅=-21A⁵, λ₆=7A⁶, λ₇=-A⁷,

Wherein

Can be obtained fom the above equation recurrence band filter certain sampled point output result only by this sampled point the first six sampling The first seven value of value and wave filter output determines；

Calculating recurrence band filter, in the output result of each sampled point, calculates its real part w_R(n) and imaginary part w_IN (), asks Going out in instantaneous amplitude testing result A (n) of each sampled point and frequency detecting result f (n), expression formula is as follows:

$A\left(n\right)=C_f{\left(w_R^2\right(n)+w_I^2(n\left)\right)}^{1/2}$

Compared with prior art, electrical network is believed by the present invention first with adding five maximum side lobe attenuation window FFT interpolation algorithms Number quickly analyze, preliminary obtain each harmonic wave and m-Acetyl chlorophosphonazo frequency content；And then utilize recurrence band filter single-point to calculate Measure little advantage, the frequency range of measured harmonic wave and m-Acetyl chlorophosphonazo be analyzed, thus in trying to achieve this frequency range contained harmonic wave and The accurate amplitude-frequency characteristic of harmonic wave.So, the traversal value of recurrence band filter mid frequency is only limitted to this frequency band range, it is to avoid Analysis to the frequency range without harmonic wave or m-Acetyl chlorophosphonazo.Greatly reduce unnecessary operand, ensure that algorithm is in DSP platform The precision of upper realization and real-time.

Accompanying drawing explanation

Fig. 1 is quick harmonic wave and the general flow chart of harmonic detection method of the present invention；

Fig. 2 present invention is recurrence band filter algorithm amplitude-versus-frequency curve when for a=0.1；

Fig. 3 is mesoscale factor pair recurrence band filter bandwidth contributions analysis of the present invention.

Detailed description of the invention

Below in conjunction with the accompanying drawing in the embodiment of the present invention, the technical scheme in the embodiment of the present invention is carried out clear, complete Describe, it is clear that described embodiment is only a part of embodiment of the present invention rather than whole embodiments wholely.Based on Embodiment in the present invention, it is every other that those of ordinary skill in the art are obtained under not paying creative work premise Embodiment, broadly falls into the scope of protection of the invention.

The quick harmonic wave of the present invention and harmonic detection method mainly comprise the steps that

S101: (sample frequency is more than 2 times of measured signal highest frequency) are to be measured under conditions of meeting sampling thheorem Signal is sampled, it is thus achieved that discrete sampling result x (n) of its signal；:

S102: the signal obtaining sampling adds five maximum side lobe attenuation window interpolation FFTs and processes: first to adopting Sample signal adds five maximum side lobe attenuation windows and is analyzed, and the signal obtained utilizes fft analysis obtain its discrete spectrum line analysis As a result, utilize interpolation algorithm to be modified the above results, obtain the substantially radio-frequency component scope of harmonic wave and m-Acetyl chlorophosphonazo；

S103: set the mid frequency that recurrence band filter performance metrics factor a, sampling period Ts and traversal calculate ω_nValue.According to harmonic wave to be obtained and m-Acetyl chlorophosphonazo Detection results, determine the bandwidth of recurrence band filter, and then setting is passed Return the bandwidth relationship of the value of scale factor a of band filter, scale factor and recurrence band filter as shown in Figure 3.Adopt The sample cycle is the inverse of 101 setting sample frequencys.Mid frequency is and signal adds five maximum side lobe attenuation window interpolation FFTs The rate-adaptive pacemaker result of gained after process；

S104: determine at different mid frequencyes according to the value of the scale factor, sampling period and the mid frequency that set The coefficient δ of recurrence band filter_iAnd λ_jValue and the amplitude-versus-frequency curve result of recurrence band filter preserving；

S105: carry out traversal filtering, applies recurrence band filter algorithm to adopt measured signal at each mid frequency Sample value x (n) calculates, and sets the seventh day of lunar month the output valve of recurrence band filter as zero, utilizes in the first six of this sampled point The first seven output valve of individual sampled value and wave filter calculates it and should calculate in the output result of the recurrence band filter of sampled point The real part of output result and imaginary part, obtain its amplitude at sampled point and the instantaneous output result of frequency；

S106: judge whether recurrence band filter completes the calculating to all sampled points；If being not fully complete, then turn S105 Continue to calculate；

S107: if completing, judges amplitude and rate-adaptive pacemaker result stable moment, takes amplitude when reaching after stable state one section Between instantaneous frequency and instantaneous amplitude output result meansigma methods as the testing result at this mid frequency；

S108: ask for the correction coefficient under set scale factor, at the beginning of carrying out amplitude detection result required above Step correction；Whether determination frequency testing result deviates the mid frequency of setting, if deviation, according to the frequency detected, calculating should Frequency and the difference of its place mid-band frequency, obtain amplitude school accurately according to wave filter amplitude-versus-frequency curve (such as Fig. 2) Positive coefficient, corrects further to amplitude detection result, it is thus achieved that the final detection result of amplitude；

S109: judged whether that the traversal to whole mid frequencyes calculates；If being not fully complete, then turn 103 and reset The value of frequency of heart, continues to calculate, if completing i.e. to obtain whole amplitudes to measured signal and frequency detecting result.

One of ordinary skill in the art will appreciate that: accompanying drawing is the schematic diagram of an embodiment, module in accompanying drawing or Flow process is not necessarily implemented necessary to the present invention.

Last it is noted that above example is only in order to illustrate technical scheme, it is not intended to limit；Although With reference to previous embodiment, the present invention is described in detail, it will be understood by those within the art that: it still may be used So that the technical scheme described in previous embodiment to be modified, or wherein portion of techniques feature is carried out equivalent；And These amendments or replacement, do not make the essence of appropriate technical solution depart from spirit and the model of embodiment of the present invention technical scheme Enclose.

Claims (4)

1. a quick harmonic wave and harmonic detection method, it is characterised in that comprise the following steps:

S101: sample to measured signal under conditions of meeting sampling thheorem, it is thus achieved that its discrete sample signals x (n), wherein N is dis-crete sample values；

S102: described discrete sample signals is added five maximum side lobe attenuation windows and is analyzed, to the windowing signal x obtained_w(n) Utilize fft analysis to obtain its discrete spectrum X (f), utilize interpolation algorithm to be modified described discrete spectrum X (f), obtain harmonic wave And the substantially radio-frequency component scope of m-Acetyl chlorophosphonazo；

S103: set the mid frequency ω that performance metrics factor a of recurrence band filter, sampling period Ts and traversal calculate_n's Value；According to harmonic wave to be obtained and m-Acetyl chlorophosphonazo Detection results, determine the bandwidth of recurrence band filter, and then it is logical to set recurrence band The value of performance metrics factor a of wave filter；

S104: according to scale factor a, sampling period Ts and the mid frequency ω that set_nValue determine at different mid frequencyes The coefficient δ of recurrence band filter_iAnd λ_jValue and the amplitude-versus-frequency curve result of recurrence band filter preserving, wherein i= 1,2,3,4,5,6, j=1,2,3,4,5,6,7；

S105: apply recurrence band filter algorithm that sampled value x (n) of measured signal is calculated at each mid frequency, if The seventh day of lunar month the output valve determining recurrence band filter is zero, utilizes before the first six sampled value and wave filter of this sampled point Seven output valves calculate its output result at the recurrence band filter of this sampling point, calculate real part and the imaginary part of output result, Obtain in instantaneous amplitude testing result A (n) of this sampled point and frequency detecting result f (n)；

S106: judge whether recurrence band filter completes the calculating to all sampled points；If being not fully complete, then proceed to step S105 continues to calculate；If completing, carry out step S107；

S107: judge amplitude and rate-adaptive pacemaker result stable moment, takes amplitude and reaches the instantaneous frequency of a period of time after stable state And the meansigma methods of instantaneous amplitude output result is as the testing result at this mid frequency；

S108: ask for the correction coefficient under set performance metrics factor a, is carried out tentatively required amplitude detection result Correction；Whether determination frequency testing result deviates the mid frequency of setting, if deviation, according to the frequency detected, calculates this frequency Rate and the difference of its place mid-band frequency, obtain amplitude rectification accurately according to wave filter amplitude-versus-frequency curve (accompanying drawing 2) Coefficient, corrects further to amplitude detection result, it is thus achieved that the final detection result of amplitude；

S109: judged whether that the traversal to whole mid frequencyes calculates；If being not fully complete, then go to step during S103 resets The value of frequency of heart, continues to calculate.

Quick harmonic wave the most according to claim 1 and harmonic detection method, it is characterised in that in described step S102 five The discrete time-domain expression formula of the maximum side lobe attenuation window of item is:

$w\left(n\right)=\underset{l=0}{\overset{4}{\Σ}}{(-1)}^l{b}_l\cos \left(\frac{2\mathrm{\π}nl}{N}\right)$

Wherein 1=0,1,2,3,4, N is maximum number of samples, b_lBeing the coefficient of five maximum side lobe attenuation windows, value is respectively as follows: B0=0.2734375, b1=0.4375, b2=0.21875, b3=0.0625, b4=0.0078125；

Described windowing signal x_w(n)=x (n) * w (n), the expression formula of described discrete spectrum X (f) is:

$X\left(f\right)=\frac{A_k}{2}\underset{i=1}{\overset{n}{\Σ}}\[W(f-f_k)e^{-j\left(\mathrm{\π}\right(f-f_k)-{\mathrm{\φ}}_k)}+W(f+f_k)e^{-j\left(\mathrm{\π}\right(f+f_k)+{\mathrm{\φ}}_k)}\]$

A in formula_kFor the amplitude of kth subharmonic, ψ_kPhase place for kth subharmonic；W (f) is five maximum side lobe attenuation window functions Frequency spectrum；f_kFor kth subfrequency by the value after frequency resolution normalization, its expression formula is f_k=kf₀×N/f_s, f₀For first-harmonic Frequency.

Quick harmonic wave the most according to claim 2 and harmonic detection method, it is characterised in that to described discrete spectrum X F step that () utilizes interpolation algorithm to be modified includes:

The spectral line of described discrete spectrum X (f) is spaced apart △ f=fs/N, and wherein fs is sample frequency；If treating measured frequency f_m=k_m*△ F, wherein k_mFor coefficient of frequency, it it is not the most integer；If f_mNeighbouring peak value spectral line is k, and corresponding amplitude is X (k), introduces frequency Rate deviation q=k_m-k-0.5, then have-0.5≤q≤0.5；

If the amplitude corresponding to spectral line of the peak value spectral line left and right sides is respectively X (k-1), X (k+1), if there being X (k-1)≤X (k+ 1), i.e.-1 spectral line of kth is time maximum spectral line near crest frequency, and making it is q with the ratio of peak value spectral line₁, the most now:

$q_1=\frac{\left|X(k-1)\right|}{\left|X\left(k\right)\right|}=\frac{\left|W(q-1)\right|}{\left|W\left(q\right)\right|}$

Can obtain the expression formula of frequency departure q:

$q=\frac{5q_1-4}{1+q_1}$

If having X (k+1)≤X (k-1), i.e.+1 spectral line of kth is time maximum spectral line near crest frequency, make itself and peak value spectral line Ratio be q₂, then have

$q_2=\frac{\left|X(k+1)\right|}{\left|X\left(k\right)\right|}=\frac{\left|W(q+1)\right|}{\left|W\left(q\right)\right|}$

Can obtain the expression formula of frequency departure q:

$q=-\frac{5q_2-4}{1+q_2}$

Then its frequency correction formula is: f=(k-q) * fs/N.

Quick harmonic wave the most according to claim 3 and harmonic detection method, it is characterised in that described step S105 bag Include:

Carry out traversal filtering, apply recurrence band filter algorithm that centrifugal pump x (n) of measured signal is entered at each mid frequency Row calculates, if instantaneous output result is w after the filtered device of discrete spectrum X (f)_η(a,ω_n, n), its expression formula is as follows:

$w_{\mathrm{\η}}(a,{\mathrm{\ω}}_n,n)=T_s\underset{i=1}{\overset{6}{\Σ}}{\mathrm{\δ}}_{i}x(n-i)-\underset{j=1}{\overset{7}{\Σ}}{\mathrm{\λ}}_j{w}_{\mathrm{\η}}(a,{\mathrm{\ω}}_n,n-j)$

Wherein δ_iWith λ_jFor constant, by performance metrics factor a given and ω_nTry to achieve, particularly as follows:

${\mathrm{\δ}}_1=\[\frac{1}{12}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4-\frac{1}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{1}{90}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A/\sqrt{a};$

${\mathrm{\δ}}_2=\[\frac{3}{4}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4-\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{19}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^2/\sqrt{a};$

${\mathrm{\δ}}_3=\[-\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4-{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{151}{45}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^3/\sqrt{a};$

${\mathrm{\δ}}_4=\[-\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4+{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{151}{45}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^4/\sqrt{a};$

${\mathrm{\δ}}_5=\[\frac{3}{4}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4+\frac{5}{6}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{19}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^5/\sqrt{a};$

${\mathrm{\δ}}_6=\[\frac{1}{12}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^4+\frac{1}{30}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^5+\frac{1}{90}{\left(\frac{{\mathrm{cT}}_s}{a}\right)}^6\]A^6/\sqrt{a};$

λ₁=-7A, λ₂=21A², λ₃=-35A³, λ₄=35A⁴, λ₅=-21A⁵, λ₆=7A⁶, λ₇=-A⁷,

Wherein $A=e^{-\frac{{\mathrm{cT}}_s}{a}+{\mathrm{j\ω}}_n{T}_s};$

Can be obtained fom the above equation recurrence band filter certain sampled point output result only by this sampled point the first six sampled value and The first seven value of wave filter output determines；

Calculating recurrence band filter, in the output result of each sampled point, calculates its real part w_R(n) and imaginary part w_IN (), obtains Instantaneous amplitude testing result A (n) of each sampled point and frequency detecting result f (n), expression formula is as follows:

$A\left(n\right)=C_f{\left(w_R^2\right(n)+w_I^2(n\left)\right)}^{1/2}$

$f\left(n\right)=\frac{1}{2\mathrm{\π}}\frac{{\tan }^{-1}\left(w_1\right(n)/w_R(n\left)\right)-{\tan }^{-1}\left(w_1\right(n-1)/w_R(n-1\left)\right)}{T_s}.$