JPWO2008136443A1 - Sinusoidal parameter estimation method - Google Patents

Abstract

An input signal is measured with respect to an unknown signal. Then, a technique for directly estimating the parameters (frequency, amplitude, phase) of the sine wave included in the input signal from the value is provided. Set up the governing equation of the analytic signal and introduce an arbitrary load function to replace it with the integral form. By solving it, the sine wave parameters are estimated. As a result, it has become possible to obtain a sine wave parameter more precisely than the conventional approximation.

Description

  The present invention relates to the fields of acoustic sensing and acoustic signal processing, and more particularly, to a method for estimating a sine wave parameter included in a signal waveform.

  Furthermore, in the present invention, the estimation of sinusoidal parameters is an extremely important technical matter in the field of sensing and signal processing of wave fields such as sound and light, and demodulation and decoding techniques in communication and data recording. It is also strongly related to the field.

  In the analysis of a signal waveform such as an audio signal, it is a basic and important matter to estimate the parameters of a sine wave included in the signal.

  Therefore, various techniques for estimating a sine wave parameter in a signal have been proposed. A widely used technique uses FFT (Fast Fourier Transform), and in particular, the method of estimating the sine wave parameter by quadratic interpolation of the obtained FFT spectrum is simple and highly accurate, so it is widely used. Has been.

Several prior arts that may be relevant to prior art literature sine wave parameter estimation are listed.

  Patent Document 1 below discloses a method for estimating the amplitude and frequency of a sine wave. The technique disclosed here is characterized in that a sine wave is estimated by interpolating between these signals based on a sampling signal.

  Patent Document 2 below discloses a spectrum analyzer that estimates an unknown signal. According to the spectrum analyzer disclosed here, the frequency and phase of the input signal can be estimated.

  Patent Document 3 below discloses a technique for extracting a sine wave of an input voice and synthesizing the voice using the sine wave.

  Japanese Patent Application Laid-Open No. 2004-228561 discloses a device that performs noise cancellation by estimating a sine wave during a noise period and interpolating an input signal when noise is mixed in an audio signal.

  Non-Patent Document 1 below discloses a study on a criterion for estimating a sine wave parameter by quadratic interpolation of FFT. In particular, research has been conducted on design parameters, window functions, window lengths, and the amount of zero padding when performing zero padding.

Japanese Patent Application Laid-Open No. 07-083964 JP 2000-258478 A JP 2004-109809 A JP 2006-135806 JP Abe, Smith, “Design criteria for sine wave parameter estimation based on FFT quadratic interpolation”, IEICE Technical Report, vol.104 No.304, pp.7-12, Nov. 2004.

  As described above, the conventional method of examining a sine wave parameter with respect to an unknown signal is an approximation using a technique such as interpolation, and it is considered that there is a certain limit in improving accuracy.

  Therefore, it is needless to say that a technique that can directly estimate the parameters (amplitude and phase) of the sine wave is desirable instead of approximation.

  Theoretically, from the Nyquist theorem, if the sampling period is sufficiently high, the original sinusoidal signal should be perfectly reproducible from the value of the sampled signal (except for the quantization error), such a technique Is desired.

  The present invention has been made in view of such a background, and is a technique for estimating a parameter of a sine wave included in an input signal directly from a value obtained by measuring the input signal, not an approximation such as interpolation. It aims to be realized.

  In order to solve the above-described problems, the present invention has a solution principle of directly estimating a sine wave parameter based on a finite interval load integral of a differential equation.

  In this patent, the parameters of the sine wave mean frequency, amplitude, and phase. In addition, in the case of a sine wave parameter, it may mean a DC bias in addition to the frequency, amplitude, and phase.

  In this patent, this DC bias is not directly handled as a parameter, but can be estimated as an amplitude in a special case where the frequency is DC.

  Specifically, the following means are adopted to solve the above problems.

  (1) In order to solve the above problems, the present invention includes a step of establishing a governing equation including a derivative of the signal in a method for estimating a parameter of a sine wave in a signal to be analyzed, and the derivative A sinusoidal parameter estimation method comprising: introducing an arbitrary load function into a governing equation and replacing it with an integral form of a finite interval; and solving an equation converted into the integral form.

  (2) In the method of estimating a parameter of a sine wave in a signal to be analyzed, the present invention constructs a governing equation including a derivative of the signal in an integral form of a finite section in which an arbitrary load function is introduced. A sinusoidal parameter estimation method comprising: a step; and a step of solving the governing equation in the integral form.

  (3) Further, in the sine wave parameter estimation method according to (1) or (2), the load function is a trigonometric function having a period of 1 / integer of the finite period. This is a sine wave parameter estimation method.

  By taking a fraction of an integer, as will be described later, it is possible to reduce the number of unknowns and counts without approximation and to leave the necessary number.

  (4) Further, in the sine wave parameter estimation method according to the above (1) or (2), the present invention introduces a load function having three different load frequencies in the step of solving the governing equation in the integral form. A sine wave parameter estimation method comprising: configuring simultaneous equations related thereto; and estimating a parameter including at least a frequency of a sine wave in the signal by solving the simultaneous equations group. is there.

  (5) Further, in the sine wave parameter estimation method according to the above (1) or (2), the present invention solves the integral-type governing equation by introducing a load function of two different load frequencies, A sine wave parameter estimation method comprising: configuring simultaneous equations related thereto; and estimating a parameter including at least a frequency of a sine wave in the signal by solving the simultaneous equations group. is there.

  (6) Further, in the sine wave parameter estimation method according to the above (1) or (2), the present invention solves the governing equation in the integral form by introducing a load function having three or more load frequencies, A sine wave parameter estimation method comprising: configuring simultaneous equations related thereto; and estimating a parameter including at least a frequency of a sine wave in the signal by solving the simultaneous equations group. is there.

  (7) Further, according to the present invention, in the sine wave parameter estimation method according to the above (1) or (2), the step of solving the governing equation in the integral form includes one load frequency and two or more different window functions. A set of simultaneous equations related to them, and estimating a parameter including at least the frequency of the sine wave in the signal by solving the set of simultaneous equations. This is a sine wave parameter estimation method.

  (8) In the sine wave parameter estimation method according to any one of (4) to (7), the step of estimating the parameter includes the estimated frequency, an integral boundary value term, and And a step of obtaining the amplitude and phase of the signal of the estimated frequency as a sine wave parameter. Here, the term “integral boundary value” refers to the difference between the signals in the integration interval and the values at both ends of the differentiation.

  (9) Further, in the sine wave parameter estimation method according to any one of (4) to (7), the step of estimating the parameter may include a phase and an amplitude with respect to the estimated frequency. A sine wave model including a step of constructing a sine wave model including parameters and fitting the sine wave model to an analysis signal to be analyzed to obtain an amplitude and a phase of the best fit sine wave model. This is a wave parameter estimation method.

  (10) Further, in the sine wave parameter estimation method according to any one of (1) to (9), the present invention provides a step or process of replacing the integral form with the integral boundary of the integral as the signal. It is a sine wave parameter estimation method characterized in that it is set in the middle of the sampling timing.

Means of Reference Invention Hereinafter, means of the reference invention related to the present invention will be described.

  The sine wave parameter estimation method described above is preferably implemented using the following hardware.

  (11) This reference invention is an apparatus for supporting the load integration of the sine wave parameter estimation method according to any one of (1) to (10), wherein N multiplications for multiplying the analysis signal by the load are performed. And N delay means for delaying each output signal of the multiplier means by N periods, and for each of the N delay means, an output after the delay. N addition means for obtaining a difference between the signal and the signal before the delay, and accumulating the difference in an accumulated value; and N accumulation means provided for each of the addition means and accumulating the output signal of the addition means. A sine wave parameter estimation method support apparatus. Here, N is a positive integer.

  (12) Further, means for constructing the governing equation of the signal in an integral form of a finite interval into which an arbitrary load function is introduced in order to implement the sine wave parameter estimation method described above, and means for solving the equation It is preferable to use a sine wave parameter estimation support device characterized by including:

  As described above, according to the present invention, the parameters of the sine wave included in the analysis signal can be directly estimated.

It is a flowchart of the flow of a process of the iterative method 2 in a present Example. It is a graph showing the estimation result of the simulated frequency in a present Example. It is a graph of the integration function (sine integration) of the sinc function in this example, and the sinc integration function (ω = 2π / 4) of the complex sine wave load. It is a graph of load distribution and error distribution when the integration range in the present embodiment is shifted to the middle position of the sampling points with respect to the sampling mesh. It is a block diagram of the structure which performs the movement finite Fourier transform in a present Example. In a present Example, it is a graph which shows the example of the window function of another solution of an equation, and a load function.

Explanation of symbols

DESCRIPTION OF SYMBOLS 10 Multiplication means 20 Delay part 22 Unit delay means 30 Accumulation part 32 Addition means

Estimation of frequency etc. of the first analysis signal
a) A signal to be analyzed by the governing equation is called an analysis signal. First, the governing equation will be described.

の一般解は、

となる。ここで、

であるから、その結果、下記式（２）

が得られる。すなわち、上で述べた式（１）（数１）が支配方程式である。複素正弦波の周波数は、この式中のａの虚部に相当し、その振幅変化率はａの実部に表現される。この支配方程式が、請求の範囲の「微分を含む支配方程式」の好適な一例に相当する。The following equation (1) for the analysis signal

The general solution of is

It becomes. here,

As a result, the following formula (2)

Is obtained. That is, the equations (1) and (Equation 1) described above are governing equations. The frequency of the complex sine wave corresponds to the imaginary part of a in this equation, and the amplitude change rate is expressed in the real part of a. This governing equation corresponds to a preferred example of the “dominating equation including differentiation” in the claims.

を導入すると、下記の同値関係式（３）によって、積分形式に置き換えることができる。 b) It is assumed that the load integration observation equation observation period is [−T / 2, T / 2] (that is, the total time is T), and the signal frequency is constant during the period. In other words, it is assumed that the differential equation described in a) above is uniformly established. This condition is an arbitrary load function

Can be replaced with the integral form by the following equivalence relation (3).

この支配方程式が、請求の範囲の「積分形式に変換した支配方程式」の好適な一例に相当する。この式（３）において、ｗ（ｔ）は、具体的には、区間［−Ｔ／２，Ｔ／２］で完備な関数系

を荷重関数として用いればよい。

This governing equation corresponds to a preferable example of “the governing equation converted into the integral form” in the claims. In this equation (3), w (t) is specifically a complete function system in the interval [−T / 2, T / 2].

May be used as a load function.

となる。ここで、

と選ぶと、これらは区間［−Ｔ／２，Ｔ／２］で完備な直交系をなし、

のように、荷重の微分による荷重積分が原関数の荷重積分に比例する。さらに、

のように、全ての

に関する積分境界値項が符号を除いて同一となる。ここで、積分境界値項とは積分区間の信号やその微分の両端の値の差を言う。以下、同様である。さてこの結果、

及び、荷重積分を、下記式（４）に示すように

と置くことによって、題意の荷重積分式は

を未知変数とする線形方程式（下記式（５））

に帰着させることができる。ａを任意の複素数としても未知数は４個であり、方程式は上記の実部と虚部の２個である。したがって、上記の線形方程式は、任意の２個の周波数で連立することによって解くことが可能になる。In order to eliminate the differentiation under this equation (3), if the subsequent integration of equation (3) is executed,

It becomes. here,

These form a complete orthogonal system in the interval [−T / 2, T / 2],

As shown, the load integral based on the load differential is proportional to the load integral of the original function. further,

Like all

The integral boundary value terms for are the same except for the sign. Here, the term “integral boundary value” refers to the difference between the signals in the integration interval and the values at both ends of the differentiation. The same applies hereinafter. Now, as a result,

And the load integral is as shown in the following formula (4).

By putting

Is an unknown variable (equation (5) below)

Can be reduced to Even if a is an arbitrary complex number, there are four unknowns, and the equations are two of the real part and the imaginary part. Therefore, the above linear equation can be solved by simultaneous operation at any two frequencies.

より、下記の式（６）

と解かれる。 c) Simultaneous method with DC component Now, when the equation of ω = 0 and the above equation are made simultaneous,

From the following formula (6)

It is solved.

にはこのｆ（ｔ）のヒルベルト変換のフーリエ変換が加わり、

も虚部を有する。解析信号化が上記の荷重積分で同時になされると考えることも可能ではあるが、それは何を解析信号であると認識するかの問題でもある。It should be noted here that f (t) is an analysis signal to be analyzed. Therefore, f (t) has an imaginary part,

Is added with the Fourier transform of the Hilbert transform of f (t),

Also has an imaginary part. Although it is possible to think that analytic signal conversion is performed simultaneously by the above-mentioned load integration, it is also a problem of what is recognized as an analytic signal.

は２倍されるだけであろうかという疑問が生じるかもしれないが、正しいヒルベルト変換ではないので２倍にはならない点にも注意すべきである。Even if the Fourier transform of the Hilbert transform is added

It may be wondered if is only doubled, but it should not be doubled because it is not a correct Hilbert transform.

となる。したがって、

となり、よって、下記（７）式、（８）式

を得る。 d) Simultaneous method using two adjacent frequency components If the two frequencies to be used are ω ₁ and ω _{2 and} s _t at these frequencies are s ₁ and s ₂ , the simultaneous equations are

It becomes. Therefore,

Therefore, the following equations (7) and (8)

Get.

とおくと（ここでＡは複素振幅を表す）、

より下記（９）式

が求められる。ここで、ａＴ＝２ｎπのときはｓｉｎ（ａＴ／２）＝０となって適用不可能となるが、この場合は、単一周波数スペクトルとなって、振幅も位相もこのスペクトル値から決定することができる。 e) Estimation of amplitude and phase Now, as the estimated frequency a

(Where A represents the complex amplitude)

From the following formula (9)

Is required. Here, when aT = 2nπ, sin (aT / 2) = 0, which is not applicable. In this case, a single frequency spectrum is obtained, and the amplitude and phase are determined from this spectrum value. Can do.

の一般解は、周知のように、下記式（１１）

となる。ここで、ＡやＢは符号も含めて任意である。すなわち、上で述べた式（１０）（数２８）が支配方程式であり、上記式（１１）の一般解の形の意味でのａの推定が周波数を推定することに相当する。この支配方程式が、請求の範囲の「微分を含む支配方程式」の好適な一例に相当する。ここでａの符号は決定しない、又は意味を有しないと言い換えてもよい。 Estimation of frequency of second real signal
a) Regarding the real signal to be analyzed by the governing equation , first, the governing equation will be explained in the same manner as described above. The following formula (10) for the actual signal

As is well known, the general solution of is the following equation (11)

It becomes. Here, A and B are arbitrary including symbols. That is, Equation (10) (Equation 2 8 ) described above is the governing equation, and estimation of a in the sense of the general solution form of Equation (11) corresponds to estimating the frequency. This governing equation corresponds to a preferred example of the “dominating equation including differentiation” in the claims. Here, the symbol a may not be determined or may be rephrased as meaningless.

を導入すると、下記の同値関係式（１２）によって、積分形式に置き換えることができる。 b) The load integration observation equation observation period is set to [−T / 2, T / 2] (that is, the total time is T) as described above, and the signal frequency is constant during the period, that is, the above a Assume that the differential equation (domination equation) described in (1) holds uniformly. This condition is the same as described above for any load function.

Can be replaced with an integral form by the following equivalence relation (12).

この支配方程式が、請求の範囲の「積分形式に変換した支配方程式」の好適な一例に相当する。この式（１２）（数３１）において、上述した荷重関数は、具体的には、区間［−Ｔ／２，Ｔ／２］で完備な関数系

を用いればよい。一般に、

であるので、前節で述べたように、

とおき、

とおくと、前節で述べたのと同様の論法によって、同値条件を構成する周波数ごとの方程式として、下記式（１３）

が得られる。

This governing equation corresponds to a preferable example of “the governing equation converted into the integral form” in the claims. In the equation (12) (Equation 31), the above-described load function is specifically a complete function system in the section [−T / 2, T / 2].

May be used. In general,

So, as mentioned in the previous section,

Toki,

By the same reasoning as described in the previous section, the following equation (13) is used as an equation for each frequency that constitutes the equivalence condition:

Is obtained.

  It should be noted that making the governing equation into an integral form using a load function is a mechanical work as described above. In particular, the fact that the period of the load function is set to an integral fraction of the integration interval also contributes to the realization of simple work. Therefore, it is preferable to implement on a computer or the like. If actually measured data is introduced into the formula constructed as described above, an equation including a sine wave parameter as a variable is completed.

  In this process, it is preferable to first determine the governing equation, and then perform a process of introducing and deforming the load function on the computer, and preparing an integral governing equation in advance, It is also preferable to configure the equation by setting the coefficient according to the load function.

  The operation that “establishes” an equation in the scope of the claims means the processing for constructing the equation as described above, and is preferably executed by a computer for quick processing.

  Further, the period of the load function may be particularly called a load period in order to distinguish it from the period of the signal (sine wave) to be processed. Similarly, the frequency of the load function may be particularly called a load frequency in order to distinguish it from the frequency of the signal of interest.

積分境界値が１個消去され、関係式は

となる。この方程式は複素変数の方程式であり、

の実部に注意し、上記方程式の実部のみを取り出すと、

が得られる。よって、以下の式（１４）のように、

と求められる。また、積分境界値は、

のように書くことができる（式（１５）及び式（１６））。 c) Simultaneous method including DC component Now, using the following relationship for ω = 0,

One integral boundary value is deleted, and the relational expression is

It becomes. This equation is a complex variable equation,

If you take out only the real part of the above equation,

Is obtained. Therefore, like the following formula | equation (14),

Is required. The integral boundary value is

(Equation (15) and Equation (16)).

であり、よってこの前段の第１式にｓ_２ω_２を乗じ、後段の第２式にｓ_１ω_１を乗じて差し引くことによって、

のように、残りの積分境界値も消去され、よって、下記式（１７）

が得られる。対称性が良く、従来の周波数成分の補間のほうほうと対応が取りやすいという利点を有する。 d) Another simultaneous method for solving the simultaneous load integral observation equation using three frequencies including a DC component is a method in which _two frequencies ω ₁ and ω ₂ are further provided in addition to the DC. By setting s _t at each frequency as s ₁ and s ₂ ,

Therefore, by multiplying the first expression of the preceding stage by s ₂ ω ₂ and multiplying the second expression of the latter stage by s ₁ ω ₁ and subtracting,

The remaining integral boundary value is also deleted as shown in FIG.

Is obtained. It has the advantage that the symmetry is good and it is easier to deal with the conventional interpolation of frequency components.

を得る。この方程式の解は下記式（１８）、式（１９）、式（２０）のようになる。 e) Simultaneous method using three adjacent frequency components Three frequencies are set as ω ₁ , ω ₂ , and ω ₃ , and s _t at that time is set as s ₁ , s ₂ , and s ₃ , respectively. By

Get. The solution of this equation is represented by the following formula (18), formula (19), and formula (20).

ただし、

である（式（２１））。なお、ここで、２×２と３×３の行列の逆行列の公式を用いている（下記数４９参照）。

However,

(Formula (21)). Here, the formula of the inverse matrix of 2 × 2 and 3 × 3 matrices is used (see Equation 49 below).

また、念のため、ω_３＝０、ｓ_３＝１の場合について求めてみると

となり、前節の結果と一致することが理解されよう。

Also, just in case, if you ask for the case of ω ₃ = 0 and s ₃ = 1

It will be understood that this is consistent with the results in the previous section.

から、振幅・位相を求めることが考えられる。すなわち、ａ＞０を推定された周波数として、

とおく。すると、

であるので、

すなわち、下記式（２２）

を求めることができ、この結果から、振幅と位相が求められる。なお、数５３の導出には、下記三角関数の和と積の公式が用いられている。 f-1) Estimation of amplitude and phase The estimated frequency and the integral boundary value term estimated at the same time

From this, it is conceivable to obtain the amplitude and phase. That is, let a> 0 be the estimated frequency,

far. Then

So

That is, the following formula (22)

From this result, the amplitude and phase can be obtained. Note that the formula of the sum and product of the following trigonometric functions is used for deriving Equation 53.

さて、ａＴ＝２ｎπのときには、ｓｉｎ（ａＴ／２）＝０となってこの方法を適用することが困難であるが、この場合は、ωＴ＝２ｎπの中の単一周波数のみにスペクトルが立つ（表れる）ので、周波数も位相もこのスペクトルから決定可能である。

Now, when aT = 2nπ, sin (aT / 2) = 0 and it is difficult to apply this method, but in this case, a spectrum stands only at a single frequency in ωT = 2nπ ( Frequency and phase can be determined from this spectrum.

  Therefore, it is preferable to use a determination step for switching or an estimation algorithm integrated until switching as the algorithm. Since it is relatively easy to detect the case where only such a spectrum stands by examining the spectrum, a person skilled in the art can easily switch the algorithm based on the determination.

へ正弦波モデル

をフィッティングする方法もある。この方法は、上述した単一の周波数のみにスペクトルが立つ場合でも、条件判断をしてアルゴリズムを切り換える必要がないので便利である。また、全ての周波数成分の影響を受ける積分境界値を用いないため、複数の周波数が混合する信号を取り扱う場合にも高い精度を維持できると考えられる。 f-2a) Other estimation methods of amplitude and phase As another method , frequency components used for estimation

Hesine wave model

There is also a method of fitting. This method is convenient because it is not necessary to judge the condition and switch the algorithm even when the spectrum stands only at the single frequency described above. Further, since the integration boundary value affected by all frequency components is not used, it is considered that high accuracy can be maintained even when a signal in which a plurality of frequencies are mixed is handled.

を展開すると、

であり、Ａで微分して零とおくと、

よって、

となる。この式では、Ａは常に実数であるが、正の値だけでなく負にもなりえる点には注意すべきである。次に、φで微分して零とおくと、

すなわち、

すなわち、

を得る。この式では、位相が［−π／２，π／２］の範囲でしか求まらないが、上記Ａの符号を正に固定し、符号を位相φに含ませると、全範囲で位相を推定することが可能である。In other words, the least squares evaluation function of the title model

Expand

And when it is differentiated by A and set to zero,

Therefore,

It becomes. Note that in this equation, A is always a real number, but it can be negative as well as positive. Next, if you differentiate it with φ and set it to zero,

That is,

That is,

Get. In this equation, the phase can be obtained only in the range of [−π / 2, π / 2]. However, if the sign of A is fixed positive and the sign is included in the phase φ, the phase is changed over the entire range. It is possible to estimate.

をＡ、Ａ^＊で微分して零とおくと、

よって、

となる。 f-2b) The other solution complex amplitude of f-2 is represented by A, and the least squares evaluation function

Is differentiated by A and A ^* and set to zero,

Therefore,

It becomes.

g) Complex sine wave simultaneous estimation method (multiple estimation method)
A case where a plurality of sine waves are included in the observation signal (also referred to as “analysis signal” or “real signal”, which means a signal to be analyzed) will be described. It is assumed that the number N of sine waves included in this observation signal is known. In many cases, the N pieces can be known in advance. For example, N can be obtained in advance from the number of peaks in the discrete spectrum distribution.

  The calculation of N is also preferably performed on a computer. Since calculating peaks from the spectrum distribution and counting the number is known as general signal processing by a computer, such a technique may be used.

とおく。この記法は、本節ｇ）においてのみ採用する記法であり、他の節では他の記法を採用する場合がある。また、正弦波を

とおくと、

と記述できる。これらによって表される各信号が満たすべき方程式は、

である。よって、これに対応する荷重積分式は、

である。この方程式は複素方程式であり、１個のωに対して２個の方程式を与える。一方、観測式としては全ての正弦波の混合信号に対する荷重積分

しか得られない。この式も複素方程式であり、同様に２個以上のωの組に対して２個の方程式を与える。これに対して、未知数は、１個の正弦波に対して、周波数

の１個と、実数である２種類の積分境界値である。したがって、この境界積分値の１＋１＝２個を足し、計３個であり、これに周波数Ｎ個を乗じて、合計３Ｎ個の未知数がある。ここで、未知数を周波数と振幅と位相で考えた場合も同様の結果になる。Now, these N frequencies are

far. This notation is used only in this section g), and other notations may be used in other sections. Also, sine wave

After all,

Can be described. The equation that each signal represented by these must satisfy is

It is. Therefore, the corresponding load integral formula is

It is. This equation is a complex equation and gives two equations for one ω. On the other hand, as an observation formula, weight integration for all sinusoidal mixed signals

Can only be obtained. This equation is also a complex equation, and similarly two equations are given for two or more sets of ω. On the other hand, the unknown is a frequency for one sine wave.

And two types of integral boundary values that are real numbers. Therefore, 1 + 1 = 2 of the boundary integral values are added to obtain a total of three, and this is multiplied by N frequencies to obtain a total of 3N unknowns. Here, the same result is obtained when the unknown is considered in terms of frequency, amplitude, and phase.

  Now, since these unknowns and load integrals are uniquely connected, two equations are established here.

となる。Ｎ＝１の場合は、この式（２３）から、Ｍ＞３／２と求まる。すなわち、２個の周波数が必要で、前節の２周波数連立法の結果と一致することが理解されよう。Therefore, if the number of ω used for estimation is M, the necessary condition for obtaining these as unknowns is the following equation (23):

It becomes. In the case of N = 1, M> 3/2 is obtained from this equation (23). That is, it will be understood that two frequencies are required and are consistent with the results of the two-frequency simultaneous method in the previous section.

  It should be noted that these equations are preferably constructed on a computer. Since obtaining an equation from a predetermined ω is a mechanical operation as described above, it can be easily executed on a general-purpose computer.

  In this patent, there is a part that says "Establish" the equation, but the work that says "Establish" in this patent is a work that "prepares" the equation in a fixed form mechanically as described above. Therefore, it can be sufficiently executed by a general-purpose computer.

  By ignoring the slight error added by the calculation of the Fourier coefficient from the sample value, the finite Fourier transform can be executed at high speed by the discrete Fourier transform. At this time, for example, if it is assumed that 256-point FFT is used, if M (number of ω used for estimation) = 127 (here, for the sake of simplicity, the actual DC component and the maximum frequency are used). N <(2M + 2) /3=257/3=85.33 is obtained. Accordingly, up to 85 sine wave parameters can be estimated in the calculation. Similarly, it is possible to estimate 42 sine wave parameters for 128-point FFT, 21 for 64 points, 10 for 32 points, 5 for 16 points, 2 for 8 points, and 1 for 4 points. it can.

  Hereinafter, specific examples will be described.

  In the implementation of the present invention, it is considered difficult to perform a direct solution, so it is considered practical to use an iterative method for improving the estimation accuracy as much as necessary.

“Solution 1” (Iteration Method 1)
Step 1
First, assuming a single frequency, the frequency is estimated from all adjacent three frequency load integral values.

Step 2
Next, the frequency and phase of the main frequency component are estimated.

Step 3
The distribution of the load integral value by the above main frequency components (including frequency perturbation) is calculated.

Step 4
The estimated value is updated according to the error from the observed distribution of the load integral value.

Step 5
The above steps 3 and 4 are repeated as many times as necessary. If the amount to be updated is less than a certain value, it is determined that a predetermined accuracy has been obtained, and the process ends.

“Solution 2” (Iteration 2)
A flowchart of the process flow of this method is shown in FIG.

Step 1
First, assuming a single frequency, the frequency, amplitude, and phase are estimated from all adjacent three frequency load integrated values. This process corresponds to S1-1 in FIG.

Step 2
Next, main frequency components are extracted by detecting peaks (corresponding to step S1-2 in FIG. 1), and the frequency, amplitude, and phase of these frequency components are estimated (S1-3 in FIG. 1). Equivalent to

Step 3
The influence of other major frequency components is subtracted from the load integral value used for the estimation of the major frequency components. This process corresponds to S1-4 in FIG.

Step 4
A new frequency, amplitude, and phase are estimated from the corrected load integral value. This process corresponds to S1-5 in FIG.

Step 5
The latter half of step 2 (S1-3), step 3 (S1-4), and step 4 (S1-5) are repeated as many times as necessary. If the amount to be subtracted is less than a certain value, it is determined that a predetermined accuracy has been obtained, and the process ends. This iterative determination process corresponds to S1-6 in FIG.

  The result of actually simulating this algorithm on a computer is shown in FIG. FIG. 2 is a graph showing frequency estimation results for single and multiple sinusoidal signals with N = 32. In this graph, the horizontal axis is a given frequency (which is a 45-degree frequency component rising to the right), and the vertical axis is an estimated frequency.

  2A shows a case where one frequency is estimated, FIG. 2B shows a case where two frequencies are estimated, and FIG. 2 (4) shows a case where four frequencies are estimated. It is.

  All adjacent three load integrals are estimated independently assuming a single frequency, and those whose detA value is equal to or greater than a predetermined first threshold value are represented by large black dots, and second values smaller than the first threshold value. Things above the threshold are indicated by small black dots. Swell-like errors occur near the intersection of multiple frequency components.

は、ｓｉｎｃ関数補間された連続関数

と等価である。よって区間［−Ｔ／２，Ｔ／２］の荷重積分は、

となる。ただし、ωＴ＝２ｎπ（ｎ：整数）のように記述される。上式（数８２）中の積分のｔ＞０をグラフで表したものが、図３に示されている。h) Improving the accuracy of the load integral using the basic value series In general, the time series of the sample values sampled at the interval Δt satisfying the Nyquist condition

Is a sinc function interpolated continuous function

Is equivalent to Therefore, the load integral in the interval [−T / 2, T / 2] is

It becomes. However, it is described as ωT = 2nπ (n: integer). A graphical representation of the integration t> 0 in the above equation (Equation 82) is shown in FIG.

FIG. 3 shows a graph of an integration function (sinusoidal integration) of a sinc function and a sinc integration function (ω = 2π / 4) of a complex sine wave load. The solid line is Si (t), the broken line is the real part of Si _u (t), and the alternate long and short dash line is the imaginary part of Si _u (t).

  That is, the value is almost zero when t <0, half when t = 0, and almost one when t> 0. However, there is undulation near the origin, and the intersection with 1 at t> 0 and the intersection with 0 at t <0 are located at approximately the middle point of the integer point of t. For this reason, the accuracy can be remarkably improved by taking the integration range in the middle of the sampling points.

  In this patent, it is assumed that the signal to be analyzed is sampled at a predetermined sampling period.

The values on the integer point of the sine integral are respectively from the origin (t = 0),
0.000000 (t = 0)
0.589490 (t = 1)
0.451412 (t = 2)
0.533094 (t = 4)
0.474970 (t = 5)
0.520108 (t = 6)
0.483206 (t = 7)
It becomes. The value at the midpoint is from t = 0.5,
0.436328 (t = 0.5)
0.511961 (t = 1.5)
0.495237 (t = 2.5)
0.502519 (t = 3.5)
0.498452 (t = 4.5)
0.501047 (t = 5.5)
It is.

  Thus, FIG. 4 shows a graph of the load distribution and the error distribution when the integration range is shifted to the middle position of the sampling points with respect to the sampling mesh.

  FIG. 4A is a graph showing an error distribution when the integration boundary is set in the middle of the sampling mesh. The horizontal axis represents sampling points, and the vertical axis represents errors. FIG. 4B is a graph showing the load distribution when the integration boundary is in the middle of the sampling mesh. The horizontal axis represents the sampling point, and the vertical axis represents the load.

  As shown in FIG. 4, the variation in the load distribution is greatly reduced, and only a small error remains at the sample point adjacent to the boundary. This load is 0.5−0.436328 = 0.063672 at the outer adjacent point of the boundary and 0.5 + 0.436328 = 0.936328 at the inner adjacent point. The above compensation corresponds to sample values of −1, 0 and N−1, N, and application requires expansion of a total of two pixels (samples) before and after. The compensation will be described next.

のように位相シフト項

が乗じられる。よって、標本化間隔Δｔを用いて、補償量は、具体的には、

となることが理解されよう。 Compensation Processing Now, in order to compensate for the influence of such a shift of the integration range (shifting the integration range to the middle position of the sample point), the following processing is performed. First, if the shift amount for the sampling mesh in the integration region is ε, the integration boundary value term is

Phase shift term as

Is multiplied. Therefore, using the sampling interval Δt, the compensation amount is specifically:

It will be understood that

h) Examination of properties of estimation method As described above, various estimation methods have been examined. It is assumed that the estimation method has the following properties.

  1. A noiseless, single-frequency sine wave with an integral Fourier coefficient gives an exact solution.

  2. Even if the discrete Fourier transform using the sample value is used, it is possible to approach the exact solution as much as possible by sufficiently increasing the sampling density.

  3. When there is no noise, the number of included frequencies is known, and the Fourier coefficient by integration is used for the observation site, an exact solution can be obtained.

Application of Third Sine Wave Frequency Estimation As described above, according to the present embodiment, the parameters of the sine wave included in the signal can be estimated (frequency, amplitude, phase). Such estimation of a sine wave should be called basic processing of audio signal processing, and its application range is extremely wide. The fields, uses, and purposes that can be applied are as follows.

  ・ Measure exact frequency estimation independent of analysis interval.

・ Measure frequency of multiple sine waves that are not in harmony ・ General harmonic analysis ・ Time / frequency analysis of voice, instrument sound, abnormal sound ・ Rise and fall analysis of sound ・ Reverberation analysis by time response ・ Frequency accuracy and frequency resolution Clear conceptual separation-Demodulation of frequency modulation, differential demodulation of phase modulation-Doppler detection, vibration detection by reflected waves, vibration detection of sound source The present invention can be applied not only to analysis of audio signals but also to various other signals Needless to say. In this patent, the parameters of the sine wave mean frequency, amplitude, and phase.

  In addition, in the case of a sine wave parameter, it may mean a DC bias in addition to the frequency, amplitude, and phase.

  In this patent, this DC bias is not directly handled as a parameter, but can be estimated as an amplitude in a special case where the frequency is DC.

The fourth moving Fourier transform time-integration type sound source identity has a remarkable feature that the window function is not required in the load integral, that is, the finite Fourier transform. We will consider using this for numerical calculation and hardware signal processing.

を乗じ、これらを区間幅Ｎで移動平均することと備えられる。荷重が「１」で一定であるため、この移動平均は、各標本点で、区間の前端と後端の差を累積メモリに加算するだけで実行可能である。これを、ここでは、移動（有限）フーリエ変換と称することにする。The above-described features are particularly prominent in the case of calculation of dense short-time Fourier transform with overlapping. That is, this calculation is a complex exponential function at each sample point.

And moving average of these with a section width N. Since the load is constant at “1”, this moving average can be executed by simply adding the difference between the front end and the rear end of the section to the accumulation memory at each sample point. This is referred to herein as a moving (finite) Fourier transform.

If the initial phase shift corresponding to the analysis interval generated at this time is compensated later, there is no problem.

であり、完全に密な移動演算とすると、

である。これに対して上記の移動フーリエ変換法の乗算回数は、

であるので、ＦＦＴをそのまま用いるよりも少ない値となる。もし、窓関数が必要であれば、窓付きでＦＦＴを完全に密に適用するには、乗算回数

が必要である。これに対して上記の移動フーリエ変換法の場合は、

であり、やはり少なくなる。 4.1 Evaluation of Computation Amount The number of multiplications and the number of additions (complex numbers) of the moving Fourier transform are N + N = 2N at each sample point. When the total data length is M, the number of multiplications is FFT with no overlap.

And if it is a completely dense movement calculation,

It is. On the other hand, the number of multiplications of the above moving Fourier transform method is

Therefore, the value is smaller than when FFT is used as it is. If a window function is required, the number of multiplications can be applied to apply FFT with a window densely.

is required. On the other hand, in the case of the above moving Fourier transform method,

And it will be less.

倍の開きがある。Ｎ＝２５６の場合は、上記値は、２５６×８＝２０４８となり、２０４８倍の開きとなる。また、窓関数なしのＦＦＴと移動フーリエ変換法とでは

倍の開きがある。もし、Ｎ＝２５６の場合は、上記値は８となり、８倍の開きとなる。Next, in the comparison of the FFT with a dense window function and the moving Fourier transform method without the window function,

There is a double opening. In the case of N = 256, the above value is 256 × 8 = 2048, which is 2048 times larger. Also, with FFT without window function and moving Fourier transform method,

There is a double opening. If N = 256, the above value is 8, which is 8 times larger.

  That is, it is expected that the speed can be increased by three orders of magnitude compared with the case where the delay estimation is performed densely by the conventional method (for example, the Fourier phase method). In addition, a speed increase of at least one digit is expected even when the sound source identity is densely performed using FFT.

とすると、区間［−Ｔ／２、Ｔ／２］の移動平均は、

であり、

である。ところが、一方、

である。ここで、

は区間内及びその両端で連続かつ微分可能である。したがって、

では、

となる。すなわち、

のような漸化式による計算が可能である。この漸化式の右辺の中括弧は単に積分の境界における信号値の差であり、区間内の値を用いずに容易に計算できる。 4.2 Calculation method of moving average recurrence formula

Then, the moving average of the section [−T / 2, T / 2] is

And

It is. However, on the other hand,

It is. here,

Is continuous and differentiable within and across the interval. Therefore,

Then

It becomes. That is,

It is possible to calculate using a recurrence formula such as The curly braces on the right side of this recurrence formula are simply the difference in signal values at the boundary of integration, and can be easily calculated without using values in the interval.

の移動平均

は

であり、

とすると、

であり、

ならば

となるからである。なお、数値計算誤差の累積を防ぐ手段を講じることが好ましい。具体的には、漸化式の直流応答を有限なものとする等の手段を高じることが好ましい。このような手段は従来から知られている工夫であるから、当業者であれば容易に実行することができる。On the other hand, when discrete values are used, the above recurrence formula is a strict formula, and an exact solution is obtained. The reason is that the interval [0, N-1]

Moving average

Is

And

Then,

And

If

Because it becomes. It is preferable to take measures to prevent accumulation of numerical calculation errors. Specifically, it is preferable to increase the means such as making the recurrence DC response finite. Such means is a conventionally known device and can be easily executed by those skilled in the art.

4.3 Example of Implementation of Moving Finite Fourier Transform Now, a block diagram of a configuration for executing the moving finite Fourier transform based on the items described so far is shown in FIG.

に

を乗じてＮ系統、Ｎ段の遅延部に入力する。As shown in FIG. 5, input data (observation signal)

In

Is input to N systems and N stages of delay units.

  The N multiplication means 10 in FIG. 5 multiply each input data by the load mentioned above. The multiplication result is supplied to the delay unit 20. The delay unit 20 is configured by connecting unit delay means 22 in N stages.

  The difference between the head and the end of the delay unit 20 becomes the input of the accumulation unit 30, and the accumulated value at each time gives the STFT output at each time. Here, the addition means 32 adds the difference between the beginning and end of the delay unit 20 and the accumulated value so far. The leading data of the delay unit 20, that is, the data to be input is supplied to the adding unit 32 without changing the sign, and the terminal data of the delay unit 20 is supplied to the adding unit 32 after the sign is inverted. Further, the output signal of the accumulating unit 30 is further supplied to the adding means 32. With such a configuration, the adding means 32 subtracts the end data from the top data to obtain a “difference”, and “accumulates” this difference and the current accumulated value of the accumulating unit 30, and obtains the result. Store again in the accumulator 30. The accumulating unit 30 is constituted by a so-called accumulator. The output signal of this accumulating unit 30 gives the STFT output.

  With the configuration shown in FIG. 5, it is possible to perform load compensation near the front end and the rear end. It should be noted that the number of stages of multiplication and addition from the delay unit is increased by the number of derivatives of the non-zero window function in a window function other than a rectangle. It should be noted that the calculation amount is reduced because the window function is unnecessary (uniformly −1 except for the boundary).

である。これに対して、移動フーリエ変換で求められるのは、

であり、ｔ＋τを新たにｔと置き換えることによって、

のように互いに変換可能であることが理解されよう。この際の補償位相は分析の離散周波数ごとに異なることに注意するべきである。積分範囲における１／２標本点シフト（上述したようなサンプリング周期の１／２の期間のシフト）のために位相補償の乗算は元々必要である。したがって、これに含ませることによって、この演算は演算量の増大には結びつかない。 4.4 Compensation of the initial phase of the analysis interval

It is. On the other hand, what is required by the moving Fourier transform is

And by replacing t + τ with t anew,

It will be understood that they can be converted to each other. It should be noted that the compensation phase at this time differs for each discrete frequency of analysis. Phase compensation multiplication is inherently necessary because of the 1/2 sample point shift in the integration range (shift of 1/2 period of the sampling period as described above). Therefore, by including it, this calculation does not lead to an increase in the amount of calculation.

5th. Reference invention So far, a sine wave parameter estimation method has been described in which parameters can be determined strictly instead of approximation, but this method is preferably implemented using the following hardware.

5.1 In order to implement the sine wave parameter estimation method described on the apparatus , means for constructing a governing equation of the signal in an integral form of a finite interval with an arbitrary load function introduced, means for solving the equation, It is preferable to use a sine wave parameter estimation support device characterized by including

  Such means is preferably realized by, for example, a digital signal processor. For example, since a semiconductor chip that performs FFT, a cryptographic chip that performs encryption processing, and the like are known, it is possible to create a semiconductor chip that performs predetermined processing at high speed in the same manner as such a semiconductor chip. preferable.

  In order to support various applications, in order to support various sampling periods and various numbers of sound sources, it is also preferable to implement it with software and a general-purpose processor that executes the software. May be slow.

  It is also preferable to configure only the number of sound sources, data input / output parts, and sampling number setting parts by software, and only the arithmetic part for solving equations using a semiconductor chip. Thus, if a part is configured by hardware and the other part is configured by software, it is possible to cope with various applications and realize high-speed data processing (signal processing).

5.2 Specific Configuration In addition, as a reference invention, in an apparatus supporting the sine wave parameter estimation method described so far, N multiplication means for multiplying the analysis signal by the load, and the N multiplications N delay means provided for each of the means for delaying each output signal of the multiplication means by N periods, and for each of the N delay means, an output signal after the delay, a signal before the delay, And N number of addition means for accumulating the difference in the accumulated value, and N number of accumulation means each provided for each of the addition means for accumulating the output signals of the addition means. Is a sine wave parameter estimation method support device. Here, N is a positive integer.

  Such a configuration is as shown in FIG.

Other Solutions of Sixth Equation Various solutions of the equations have been described in the above-mentioned Chapter 2 (c) Simultaneous method including DC component, (d) Simultaneous method using three frequencies including DC component, and the like. However, there are other methods for solving simultaneous equations, and other methods will be described below.

において、

かつ、この窓関数ｐ（ｔ）を原点対称に選ぶと、

である。ゆえに、荷重積分観測方程式は

となる。ただし、

である。特に

となるように、窓関数を選ぶと、荷重積分観測方程式は、

となり、

のように求められる。ｐ（ｔ）の具体例としてＨａｎｎの窓関数を利用することが考えられる。Ｈａｎｎの窓関数を利用し、

とすると（ここで、

である）、原点対象で

となる。また、

である。これらのグラフが図６に示されている。図６のグラフは、横軸が時間を表し、縦軸が荷重を表す。図６（ａ）は、窓関数を示し、図６（ｂ）は荷重関数（Ｔ＝８）を表す。図６（ａ）では、実線が

を表し、破線が

を表し、一点鎖線が

をそれぞれ表す。図６（ｂ）は、これらに、ｎ＝４のｃｏｓ正弦波荷重関数を乗じたものである。a) Method using single load frequency and window function In this section as well, we will introduce a method to remove the integral boundary value from unknowns, with the knowledge that the type of load integral observation value increases. Load integral observation equation as a starting point

In

And if this window function p (t) is selected to be symmetrical with the origin,

It is. Therefore, the load integral observation equation is

It becomes. However,

It is. In particular

If you choose a window function so that

And

It is required as follows. As a specific example of p (t), it is conceivable to use Hann's window function. Using Hann's window function,

(Where

At the origin

It becomes. Also,

It is. These graphs are shown in FIG. In the graph of FIG. 6, the horizontal axis represents time, and the vertical axis represents load. FIG. 6A shows a window function, and FIG. 6B shows a load function (T = 8). In FIG. 6A, the solid line is

Represents a broken line

Represents an alternate long and short dash line

Respectively. FIG. 6B is obtained by multiplying these by a cos sine wave load function of n = 4.

Claims (10)

In a method for estimating a parameter of a sine wave in a signal to be analyzed,
Establishing a governing equation including a derivative of the signal;
Introducing an arbitrary load function into the governing equation including the derivative and replacing it with an integral form of a finite interval;
Solving a governing equation converted to the integral form;
A sine wave parameter estimation method comprising: In a method for estimating a parameter of a sine wave in a signal to be analyzed,
Constructing a governing equation including a derivative of the signal in an integral form of a finite interval with an arbitrary weight function introduced;
Solving the governing equation in integral form;
A sine wave parameter estimation method comprising: In the sine wave parameter estimation method according to claim 1 or 2,
The sine wave parameter estimation method, wherein the load function is a trigonometric function having a period of 1 / integer of the finite period. In the sine wave parameter estimation method according to claim 1 or 2,
Solving the governing equation in integral form includes
Introducing load functions of three different load frequencies and constructing simultaneous equations related to them,
Estimating a parameter including at least a frequency of a sine wave in the signal by solving the group of simultaneous equations;
A sine wave parameter estimation method comprising: In the sine wave parameter estimation method according to claim 1 or 2,
Solving the governing equation in integral form includes
Introducing load functions of two different load frequencies and constructing simultaneous equations related to them,
Estimating a parameter including at least a frequency of a sine wave in the signal by solving the group of simultaneous equations;
A sine wave parameter estimation method comprising: In the sine wave parameter estimation method according to claim 1 or 2,
Solving the governing equation in integral form includes
Introducing load functions of three or more load frequencies and constructing simultaneous equations related to them,
Estimating a parameter including at least a frequency of a sine wave in the signal by solving the group of simultaneous equations;
A sine wave parameter estimation method comprising: In the sine wave parameter estimation method according to claim 1 or 2,
Solving the governing equation in integral form includes
Introducing a load function of one load frequency and two or more different window functions and constructing simultaneous equations related to them,
Estimating a parameter including at least a frequency of a sine wave in the signal by solving the group of simultaneous equations;
A sine wave parameter estimation method comprising: In the sine wave parameter estimation method according to any one of claims 4 to 7,
Estimating the parameter comprises:
Obtaining the estimated frequency, the integral boundary value term, and the amplitude and phase of the signal of the estimated frequency as parameters of a sine wave;
A sine wave parameter estimation method comprising: Here, the term “integral boundary value” refers to the difference between the signals in the integration interval and the values at both ends of the differentiation. In the sine wave parameter estimation method according to any one of claims 4 to 7,
Estimating the parameter comprises:
By constructing a sine wave model including the phase and amplitude as parameters with respect to the estimated frequency and fitting the sine wave model to an analysis signal to be analyzed, the amplitude and phase of the best-fitting sine wave model can be obtained. The process of seeking
A sine wave parameter estimation method comprising: In the sine wave parameter estimation method according to any one of claims 1 to 9,
The process or process of replacing the integral form is as follows:
A method for estimating a sine wave parameter, wherein an integration boundary of the integration is set in the middle of a timing for sampling the signal.

Images