JP2018097430A - Time series signal feature estimating apparatus and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a time series signal feature estimation technique capable of calculating a time series signal feature in a feature analysis of a time series signal based on a short-time Fourier transform in a form with less uncertainty.SOLUTION: A time series signal feature estimating apparatus comprises: a logarithmic amplitude time difference calculation part 110 for calculating a logarithmic amplitude time difference ΔlogA/Δtthat is a time difference in logarithmic amplitudes logAat discrete times tfor discrete angular frequencies ωfrom amplitudes Aat the discrete times tand the discrete angular frequencies ω; and a time series signal feature calculating part 120 for calculating, as a time series signal feature, a group delay characteristic GDat the discrete times tand the discrete angular frequencies ωusing a predetermined equation from the logarithmic amplitude time difference ΔlogA/Δt.SELECTED DRAWING: Figure 1

Description

  The present invention relates to a time series signal feature estimation technique, and more particularly to a technique for estimating a feature of a time series signal based on a short-time Fourier transform.

  One technique for analyzing time-series signals such as audio signals and acoustic signals is short-time Fourier transform (see Non-Patent Document 1). The time series signal to be analyzed is converted into a two-dimensional complex signal having indexes of time and angular frequency (time and frequency) as a result of short-time Fourier transform.

  Instead of using the real and imaginary parts of this complex signal as they are for analysis, convert them into amplitude components and phase components, and pay attention to the characteristics of both components or only one of them. Sometimes. In addition, the group delay characteristic which is a value obtained by reversing the positive / negative of the differential with respect to the angular frequency of the phase component, and the instantaneous angular frequency which is the differential value with respect to time of the phase component (the value obtained by dividing the instantaneous angular frequency by 2π Analysis may be performed based on the instantaneous frequency.

M. Parchami, WP Zhu, B. Champagne, E. Plourde, “Recent Developments in Speech Enhancement in the Short-Time Fourier Transform Domain”, IEEE Circuits and Systems Magazine, Vol.16, No.3, pp.45-77 , 2016.

  Of the two components, the amplitude component is often employed because of the ease of physical interpretation of the obtained features and the simplicity of numerical calculations.

  On the other hand, the phase component, group delay characteristics and instantaneous angular frequency (instantaneous frequency) obtained from angular frequency differentiation and time differentiation are expected to provide useful characteristics, but are calculated directly from complex signals. Since the phase component is in a state where the value is wrapped between −π and π, it becomes a two-dimensional sequence having discontinuities in both the time axis direction and the frequency axis direction. For this reason, in order to calculate the group delay characteristics and instantaneous angular frequency values, it is necessary to know the characteristics of the unwrapped phase components by eliminating the discontinuity of the phase components. It does not have a unique solution and becomes an uncertain problem.

  SUMMARY OF THE INVENTION An object of the present invention is to provide a time-series signal feature estimation technique capable of calculating a time-series signal feature with less uncertainty in a time-series signal feature analysis based on a short-time Fourier transform. To do.

One aspect of the present invention is that a discrete time t _n of a discrete short-time Fourier transform of a time-series signal analyzed using a Gauss window with a standard deviation σ (σ is a positive real number) as a window function, a discrete angular frequency ω _k ( n is an integer representing a discrete time index, k is an integer representing a discrete angular frequency index)

(Where f _s is the sampling frequency, N is the discrete Fourier transform score), the amplitude at the discrete time t _n and the discrete angular frequency ω _k is A _{n, k} , and the amplitude A _{n, k} (n ₁ , n ₂ is an integer satisfying n ₂ −n ₁ ≧ 1, n = n ₁ , n ₁ +1,..., n ₂ , k is an integer), and the logarithmic amplitude logA _{n at the} discrete time t _{n with} respect to the discrete angular frequency ω _{k ,} logarithmic magnitude time differential delta _t logA _n is the time difference _k, and the logarithmic amplitude time difference calculating section for calculating the _k / delta] t _n, the logarithmic amplitude time difference delta _t logA _n, from _k / delta] t _n, the following equation A time series signal feature calculation unit for calculating the group delay characteristic GD _{n, k} at the discrete time t _n and the discrete angular frequency ω _k as a time series signal feature using

(Where T is the effective width of the window function, T = N / f _s ).

One aspect of the present invention is that a discrete time t _n of a discrete short-time Fourier transform of a time-series signal analyzed using a Gauss window with a standard deviation σ (σ is a positive real number) as a window function, a discrete angular frequency ω _k ( n is an integer representing a discrete time index, k is an integer representing a discrete angular frequency index)

(Where f _s is the sampling frequency, N is the discrete Fourier transform score), the amplitude at the discrete time t _n and the discrete angular frequency ω _k is A _{n, k} , and the amplitude A _{n, k} (n is an integer) , k _1, k ₂ is an integer satisfying _{k 2 -k 1 ≧ 1, k} = k 1, k 1 + 1, ..., a k _2), the logarithmic amplitude logA _n in a discrete angular frequency omega _k for discrete time t _{n ,} logarithmic magnitude angular frequency difference is the angular frequency difference between the _{k Δ} ω logA _n, and the logarithmic amplitude angular frequency difference calculation unit for calculating the _k / [delta] [omega _k, the logarithmic amplitude angular frequency difference delta _omega logA _n, from _k / [delta] [omega _k A time series signal feature calculation unit for calculating the instantaneous frequency IF _{n, k} at the discrete time t _n and the discrete angular frequency ω _k as a time series signal feature using the following equation:

including.

One aspect of the present invention is that a discrete time t _n of a discrete short-time Fourier transform of a time-series signal analyzed using a Gauss window with a standard deviation σ (σ is a positive real number) as a window function, a discrete angular frequency ω _k ( n is an integer representing a discrete time index, k is an integer representing a discrete angular frequency index)

(Where f _s is the sampling frequency, N is the discrete Fourier transform score), the amplitude at the discrete time t _n and the discrete angular frequency ω _k is A _{n, k} , and the amplitude A _{n, k} (n ₁ , n ₂ is an integer satisfying n ₂ −n ₁ ≧ 1, n = n ₁ , n ₁ +1,..., n ₂ , k = 0, 1,. , N-2 group delay for calculating the group delay characteristic integral value Σ (-GD) _{n, k} , which is the integral value of the sign inversion value of the group delay characteristic from discrete angular frequency ω _k to ω _{k + 1} A characteristic integral value calculation unit, and k ₀ is an integer between ₀ and N−1, the initial phase value is φ _{n, k_0} , and the group delay characteristic integral value Σ (−GD) _{n, k} (k = k ₀ , k ₀ +1, ..., N-2), the phase φ _{n, k} (k = k ₀ +1, ..., N-1) in ascending order from k = k ₀ +1 Calculated as a sequence signal feature,

From the group delay characteristic integral value Σ (−GD) _{n, k} (k = 0, 1,..., K ₀ −1), using the following equation, the phase φ _{n, k in} descending order from k = k ₀ −1 a time series signal feature calculating unit for calculating (k = 0, 1,..., k ₀ -1) as a time series signal feature;

including.

One aspect of the present invention is that a discrete time t _n of a discrete short-time Fourier transform of a time-series signal analyzed using a Gauss window with a standard deviation σ (σ is a positive real number) as a window function, a discrete angular frequency ω _k ( n is an integer representing a discrete time index, k is an integer representing a discrete angular frequency index)

(Where f _s is the sampling frequency, N is the discrete Fourier transform score), the amplitude at the discrete time t _n and the discrete angular frequency ω _k is A _{n, k} , and the amplitude A _{n, k} (n ₀ is Integer, m is a positive integer, n = n ₀ , n ₀ +1, ..., n ₀ + m, k ₁ , k ₂ are integers satisfying k ₂ -k ₁ ≧ 1, k = k ₁ , k ₁ + 1, ..., k ₂ ) to n = n ₀ , n ₀ +1, ..., n ₀ + m-1 for each n, the integral value of 2π times the instantaneous frequency from discrete time t _n to t _{n + 1} Instantaneous frequency integrated value Σ (2πIF) _{n, k} , and an initial phase value φ _{n_0, k} , and the instantaneous frequency integrated value Σ (2πIF) _{n, k} (n = n ₀ , n ₀ +1, ..., n ₀ + m-1), using the following formula, the phase φ _{n, k} (n = n ₀ +1, ..., n ₀ + in ascending order from n = n ₀ +1) m) a time series signal feature calculation unit for calculating time series signal features;

including.

  According to the present invention, by calculating the phase component, group delay characteristic, and instantaneous frequency based on the amplitude component of the short-time Fourier transform of the time series signal, it is possible to calculate these values with less uncertainty. It becomes possible.

The figure which shows an example of a structure of the time series signal feature estimation apparatus. The figure which shows an example of operation | movement of the time series signal feature estimation apparatus. The figure which shows an example of a structure of the time series signal feature estimation apparatus. The figure which shows an example of operation | movement of the time series signal feature estimation apparatus. The figure which shows an example of a structure of the time series signal feature estimation apparatus 300. FIG. The figure which shows an example of operation | movement of the time series signal feature estimation apparatus 300. The figure which shows an example of a structure of the time series signal feature estimation apparatus 400. The figure which shows an example of operation | movement of the time series signal feature estimation apparatus 400. The figure which shows the estimated value of the true phase and group delay characteristic in the discrete angular frequency index k = 0.

  Hereinafter, embodiments of the present invention will be described in detail. In addition, the same number is attached | subjected to the structure part which has the same function, and duplication description is abbreviate | omitted.

<Notation>
_ (Underscore) represents a subscript. For example, ^xy_z represents that _yz is a superscript to x, and _{xy_z} represents that _yz is a subscript to x.

<Principle of time series signal feature estimation based on short-time Fourier transform>
First, an explanation will be given of the group delay characteristics, the instantaneous angular frequency, and the equation satisfied by the phase, which are the principle of time series signal feature estimation based on the short-time Fourier transform. By using these equations, the phase component of a two-dimensional complex signal obtained by short-time Fourier transform of a time-series signal can be estimated from the amplitude component of the complex signal instead of directly calculating from the complex signal value itself. Is possible. Further, the group delay characteristic and the instantaneous angular frequency can be estimated from the amplitude component.

[Equation satisfied by group delay characteristics, instantaneous angular frequency, and phase]
First, for a continuous time signal x (τ) that is a time series signal, a short-time Fourier transform X _{t, ω} at the time t and the angular frequency ω given by Equation (1) is considered.

However, w (t) is the window function that takes the maximum value at t = 0, and T is the effective width of the window function. J is an imaginary unit.

At this time, when the short-time Fourier transform X _{t, ω} is partially differentiated with _respect to time t, equation (2) is obtained.

Similarly, Equation (3) is obtained by partial differentiation of the short-time Fourier transform X _{t, ω} with the angular frequency ω.

  Here, a Gaussian window given by Equation (4) is selected as the window function.

However, the standard deviation σ of the Gaussian window is a real number (positive real number) greater than 0.

  At this time, if the window function w (τ-t) is partially differentiated with respect to time t, the following equation (5) is obtained.

  At this time, when Equation (4) is applied to Equation (3), Equation (5) is applied to Equation (2), and the common term is eliminated, Equation (6) is obtained.

If the amplitude of the short-time Fourier transform X _{t, ω} is At _{, ω} and the phase is φ _{t, ω} , Equation (7) is established. Further, when the logarithm of equation (7) is taken, equation (8) is obtained.

Here, when logX _{t, ω} is partially differentiated with respect to time t, equation (9) is obtained.

  Therefore, Expression (10) is established from Expression (8) and Expression (9).

Similarly, when logX _{t, ω} is partially differentiated with respect to the angular frequency ω, Expression (11) is obtained.

  Therefore, Expression (12) is established from Expression (8) and Expression (11).

  Substituting Equation (10) and Equation (12) into Equation (6) yields Equation (13).

  Here, since the real part and the imaginary part on the left side of Expression (13) are independently 0, Expression (14) and Expression (15) are obtained.

  The left side of Expression (14) is the definition of the group delay characteristic itself, and it can be seen from the right side of Expression (14) that the group delay characteristic can be calculated only from the logarithmic amplitude that is amplitude information without using the phase. Also, the left side of equation (15) is the definition of instantaneous angular frequency itself, and it can be seen from the right side of equation (15) that the instantaneous angular frequency can also be calculated from only the logarithmic amplitude without using the phase. Since the instantaneous frequency is obtained by dividing the instantaneous angular frequency by 2π, it can be seen that the instantaneous frequency can be calculated from only the logarithmic amplitude.

  Further, when equation (14) is integrated with respect to the angular frequency, equation (16) is obtained.

However, C ₁ is a constant of integration.

  Similarly, when equation (15) is integrated over time, equation (17) is obtained.

However, C ₂ is a constant of integration.

The integration constants C ₁ and C ₂ need to be determined separately.

It can be seen that the left side of the equations (16) and (17) is the phase φ _{t, ω} , and the phase φ _{t, ω} can also be calculated only from the logarithmic amplitude.

  From the above, by selecting a Gaussian window as the window function of the short-time Fourier transform, the group delay characteristic is expressed by equation (14), the instantaneous angular frequency (instantaneous frequency) is expressed by equation (15), and the phase is expressed by equation (16) or It can be seen that each can be calculated from the logarithmic amplitude using equation (17).

  The group delay characteristic and the instantaneous angular frequency (instantaneous frequency) are calculated by using the amplitude component instead of the wrapped phase component, thereby making it less susceptible to the uncertainty problem. In addition, the calculation of the phase itself can be obtained only in a wrapped state when calculating directly from a complex signal, but it is easy to understand physical characteristics by using Equation (16) or Equation (17). It can be calculated in an unwrapped state.

  The above-described equations satisfying the group delay characteristics, instantaneous angular frequency, and phase, which are the principle of feature estimation of the time series signal based on the short-time Fourier transform, describe the short-time Fourier transform for the continuous-time signal. However, when actually calculating with a computer, it is necessary to use the result of the short-time Fourier transform of the digitized digital signal.

  Therefore, in the following, an approximate expression for discretizing the time t and the angular frequency ω that take continuous values and estimating the group delay characteristic, instantaneous frequency, and phase will be described.

[Approximation formula for estimating group delay characteristics, instantaneous frequency, and phase]
The discrete Fourier transform is applied to the short-time Fourier transform X _{t, ω} using a Gaussian window whose standard deviation is σ (σ is a positive real number) of the continuous-time signal x (τ) that is a time series signal. Discrete short-time Fourier transform X _{n, k} obtained by application (n is an integer representing a discrete time index, k is an integer representing a discrete angular frequency index), that is, standard deviation σ (σ is a positive real number) as a window function Consider a discrete short-time Fourier transform X _{n, k} of a time-series signal analyzed using a Gaussian window.

As a result, the discrete time t _n obtained by discretizing the time t and the discrete angular frequency ω _{k obtained} by discretizing the angular frequency ω are expressed by Expression (18) and Expression (19), respectively.

Here, f _s is a sampling frequency, and N is a discrete Fourier transform (DFT) score.

  At this time, the effective width T of the window function is expressed by Expression (20).

With respect to the discretized digital signal x (τ _i ), the discrete short-time Fourier transform X _{n, k at} the discrete time t _n and the discrete angular frequency ω _k is expressed by Expression (21).

  Equation (21) corresponds to Equation (1).

(Approximate formula for estimating group delay characteristics)
By discretizing Equation (14), an approximate equation for estimating the group delay characteristic is obtained.

_Let A _{n, k be the} amplitude at discrete time t _n and discrete angular frequency ω _k . In addition, a logarithmic amplitude time difference that is a time difference of the logarithmic amplitude logA _{n, k at} the discrete time t _{n with} respect to the discrete angular frequency ω _k is denoted by Δ _t logA _{n, k} / Δt _n .

When the central difference approximation is used as the logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _{n at the} discrete time t _n for the discrete angular frequency ω _k , Equation (22) is obtained.

However, δ ₁ and δ ₂ are non-zero numbers introduced for numerical stabilization.

From Expression (14) and Expression (22), Expression (23), which is an approximate expression for estimating the group delay characteristic GD _{n, k} at the discrete time t _n and the discrete angular frequency ω _k , is obtained.

Further, from the relationship of Equation (24), Equation (25) is obtained instead of Equation (22) as the logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n using center difference approximation.

However, δ ₃ is a non-zero number introduced for numerical stabilization. Note that _{Δ t A n, k / Δt} n is an amplitude time difference is a time difference between the amplitudes A _{n, k} at discrete time t _n for discrete angular frequency omega _k.

From Expression (14) and Expression (25), Expression (26) which is an approximate expression for estimating the group delay characteristic GD _{n, k} at the discrete time t _n and the discrete angular frequency ω _k is obtained.

Here, the logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n is calculated using center difference approximation, but forward difference approximation or backward difference approximation may be used instead of the center difference approximation. For example, forward difference approximation and backward difference approximation corresponding to Equation (22) are as shown in Equation (22a) and Equation (22b), respectively.

  Similarly, forward difference approximation and backward difference approximation corresponding to Equation (25) are as shown in Equation (25a) and Equation (25b), respectively.

  Here, with respect to the standard deviation σ of the Gaussian window, assuming that α times σ (where α is an integer of 1 or more) corresponds to half the effective width of the window function, the standard deviation σ of the Gaussian window is Equation (27) is satisfied.

  For example, when α = 4, the standard deviation σ corresponding to α can be obtained by Expression (28).

  When the above assumption is introduced, there is an advantage that it is easy to intuitively interpret the validity of the selected σ.

(Approximate expression for estimating instantaneous frequency)
The approximate expression for estimating the instantaneous frequency is obtained by discretizing the expression (15) and dividing by 2π.

The logarithmic amplitude angular frequency difference that is the angular frequency difference of the logarithmic amplitude logA _{n, k at} the discrete angular frequency ω _{k with} respect to the discrete time t _n is denoted by Δ _ω logA _{n, k} / Δω _k .

Log magnitude angular frequency difference delta _omega logA _n in a discrete angular frequency omega _k for discrete time t _n, the use of central difference approximation as _k / _{Δω k,} obtain equation (29).

However, δ ₄ and δ ₅ are non-zero numbers introduced for numerical stabilization.

From Expression (15) and Expression (29), Expression (30), which is an approximate expression for estimating the instantaneous frequency IF _{n, k} at the discrete time t _n and the discrete angular frequency ω _k , is obtained.

Further, from the relationship of Equation (31), Equation (32) is obtained instead of Equation (29) as the logarithmic amplitude angular frequency difference Δ _ω logA _{n, k} / Δω _k using center difference approximation.

However, δ ₆ is a non-zero number introduced for numerical stabilization. Note that _{Δ ω A n, k / Δω} k is the amplitude angular frequency difference is the angular frequency difference between the amplitudes A _{n, k} at discrete angular frequency omega _k for discrete time t _n.

From Expression (15) and Expression (32), Expression (33) which is an approximate expression for estimating the instantaneous frequency IF _{n, k} at the discrete time t _n and the discrete angular frequency ω _k is obtained.

Here, the logarithmic amplitude angular frequency difference delta _omega logA _n using central difference approximation _has been calculated _k / [Delta] [omega _k, instead of the central difference approximation, it may be used forward difference approximation or backward difference approximation. For example, forward difference approximation and backward difference approximation corresponding to Equation (29) are as shown in Equation (29a) and Equation (29b), respectively.

  Similarly, forward difference approximation and backward difference approximation corresponding to Equation (32) are as shown in Equation (32a) and Equation (32b), respectively.

(Approximation formula for estimating phase)
In order to obtain an approximate expression for estimating the phase, in the following, integral calculation is performed by trapezoidal approximation for Equation (23), Equation (26), Equation (30), and Equation (33), and obtained in the form of a recursive equation. .

  Expressions (23) and (26) are integral approximated in the frequency direction. The equation for integral calculation by trapezoidal approximation in the frequency direction is as shown in Equation (34).

Equation (34) is a discrete time t _n, the phase phi _{n, k} at discrete angular frequency omega _k is discrete time t _n, the phase φ _{n, k-1} in discrete angular frequency omega _k-1, a discrete angular frequency omega It shows that the integral value of the sign inversion value (see equation (14)) of the group delay characteristic from _k−1 to ω _k is added. The second term on the right side of Equation (34) corresponds to the area of the trapezoid.

The phase φ _{n, k at the} discrete time t _n and the discrete angular frequency ω _k is expressed by the approximate expression of the expression (35) by applying the expression (23) to the expression (34).

However, δ ₇ and δ ₈ are non-zero numbers introduced for numerical stabilization.

Similarly, the phase φ _{n, k} at the discrete time t _n and the discrete angular frequency ω _k is expressed by the approximate expression of Expression (36) by applying Expression (26) to Expression (34).

  Further, Expression (37) is obtained from Expression (34).

Equation (37) is a discrete time t _n, the phase phi _{n, k} at discrete angular frequency omega _k is discrete time t _n, from the phase φ _{n, k + 1} at discrete angular frequency omega _{k + 1,} the discrete angular frequency omega It shows that the integral value of the sign inversion value (see equation (14)) of the group delay characteristic from _k to ω _{k + 1} is subtracted. The second term on the right side of Equation (37) corresponds to the area of the trapezoid. Using equation (37), it can be seen that the phase can be obtained in such a way that the calculation proceeds in the negative direction with respect to the frequency direction.

The phase φ _{n, k at the} discrete time t _n and the discrete angular frequency ω _k is an approximate expression for obtaining the formula (35) and the formula (36) from the formula (35) and the formula (36), respectively. 39) is obtained.

  Next, with respect to Equation (30) and Equation (33), integral approximation is performed in the time direction. The formula for integral calculation by trapezoidal approximation in the time direction is as shown in Equation (40).

Equation (40) is a discrete time t _n, the phase phi _{n, k} at discrete angular frequency omega _k is discrete time t _n-1, the phase φ _{n-1, k} at discrete angular frequency omega _k, discrete time t _n 2π times the instantaneous frequency of the toward t _{n -1} indicates that a plus integral value (equation (15) refer). The second term on the right side of Equation (40) corresponds to the area of the trapezoid.

The phase φ _{n, k at the} discrete time t _n and the discrete angular frequency ω _k is expressed by the approximate expression of the expression (41) by applying the expression (30) to the expression (40).

However, δ ₉ and δ ₁₀ are non-zero numbers introduced for numerical stabilization.

Similarly, the phase φ _{n, k} at the discrete time t _n and the discrete angular frequency ω _k is expressed by the approximate expression of the expression (42) by applying the expression (33) to the expression (40).

  Recursive approximations for estimating the phases of Equations (35) and (36) (Equations (38) and (39)) and Equations (41) and (42) are obtained by integral calculation using trapezoidal approximation. However, an integral calculation method other than the trapezoidal approximation may be used. For example, using numerical integration formulas of integrands that are discretized at regular intervals, commonly known as Newton-Cotes formulas, such as the Simpson formula, formulas (35) and (36) (formula (38) and A recursive approximate expression for estimating a phase such as Expression (39)), Expression (41), and Expression (42) may be obtained.

In general, the integral value of the sign inversion value of the group delay characteristic from the discrete angular frequency ω _k-1 to ω _k (hereinafter, the integral value of the sign inversion value of the group delay characteristic is simply referred to as the group delay characteristic integral value) is represented by Σ (− Assuming that (GD) _{n, k−1} , the equation (43) is obtained as an equation corresponding to the equation (34). Here, the group delay characteristic integral value Σ (−GD) _{n, k−1} is a value calculated using the group delay characteristic GD _{n, k−1} and the group delay characteristic GD _{n, k} .

If you look at equation (23) or equation (26) (or equation (22), equation (22a), equation (22b), equation (25), equation (25a), equation (25b)) As can be seen, the group delay characteristic GD _{n, k} can generally be calculated using at least two values of amplitude A _{n−1, k} , amplitude An _{n, k} , and amplitude A _{n + 1, k.} . Therefore, the group delay characteristic integral value Σ (−GD) _{n, k−1} is at _{least one} of the amplitudes A _{n−1, k−1} , amplitudes A _{n, k−1} , and amplitudes A _{n + 1, k−1.} It is also a value calculated by using at least two values among the two values, amplitude A _{n−1, k} , amplitude A _{n, k} , and amplitude A _{n + 1, k} .

Similarly, when the integral value (group delay characteristic integral value) of the sign inversion value of the group delay characteristic from the discrete angular frequency ω _k to ω _{k + 1} is Σ (-GD) _{n, k} , it corresponds to Equation (37) Equation (44) is obtained as Here, the group delay characteristic integral value Σ (−GD) _{n, k} is a value calculated using the group delay characteristic GD _{n, k} and the group delay characteristic GD _{n, k + 1} .

The group delay characteristic integral value Σ (-GD) _{n, k} is an amplitude An _{n-1, k} , an amplitude An _{n, k} , an amplitude An _{n + 1, k} , and an amplitude _An It is also a value calculated using at least two values of _{−1, k + 1} , amplitude A _{n, k + 1} , and amplitude A _{n + 1, k + 1} .

Also, the integral value 2π times the instantaneous frequency from discrete time t _n-1 to t _n (hereinafter, the integral value 2π times the instantaneous frequency is simply referred to as the instantaneous frequency integral value) is Σ (2πIF) _{n-1, k} Then, Expression (45) is obtained as an expression corresponding to Expression (40). Here, the instantaneous frequency integral value Σ (2πIF) _{n−1, k} is a value calculated using the instantaneous frequency IF _n− 1k and the instantaneous frequency IF _{n, k} .

If you look at equation (30) or equation (33) (or equation (29), equation (29a), equation (29b), equation (32), equation (32a), equation (32b)) As can be seen, the instantaneous frequency IF _{n, k} can generally be calculated using at least two values of the amplitude An _{n, k−1} , the amplitude An _{n, k} , and the amplitude An _{n, k + 1} . Therefore, the instantaneous frequency integrated value Σ (2πIF) _{n−1, k} is at least two of the amplitudes A _{n−1, k−1} , amplitudes A _{n−1, k} , and amplitudes A _{n−1, k + 1} . It is also a value calculated using at least two values among the value and the amplitude _{An, k-1} , amplitude _{An, k} , and amplitude _{An, k + 1} .

<First embodiment>
Hereinafter, the time-series signal feature estimation apparatus 100 will be described with reference to FIGS. As illustrated in FIG. 1, the time-series signal feature estimation apparatus 100 includes a logarithmic amplitude time difference calculation unit 110, a time-series signal feature calculation unit 120, and a recording unit 190. The recording unit 190 is a component that appropriately records information necessary for processing of the time-series signal feature estimation device 100. For example, the recording unit 190 records in advance the sampling frequency f _s , the discrete Fourier transform score N, and the Gaussian window standard deviation σ. Of course, the sampling frequency f _s , the discrete Fourier transform score N, and the standard deviation σ of the Gaussian window may be given as inputs to the time-series signal feature estimation apparatus 100 instead of being recorded in the recording unit 190 in advance.

The amplitudes A _{n, k} (n ₁ , n ₂ are integers satisfying n ₂ -n ₁ ≧ 1, n = n ₁ , n ₁ +1,..., N ₂ , k, which are inputs to the time-series signal feature estimation apparatus 100 Is an integer). For example, in order to reduce memory and transmission capacity _, only the amplitude _{An, k} may be recorded in the recording unit 190 from the result of performing a short-time Fourier transform on the time series signal in advance. Further, the amplitude _{An, k} may be transmitted from the outside before the processing by the time-series signal feature estimation apparatus 100 is started.

The time-series signal feature estimation apparatus 100 calculates a time series at discrete time t _n and discrete angular frequency ω _k from amplitudes A _{n, k} (n = n ₁ , n ₁ +1,..., N ₂ , k is an integer). The group delay characteristic GD _{n, k} which is a signal feature is estimated and output.

The operation of the time-series signal feature estimation apparatus 100 will be described with reference to FIG. The logarithmic amplitude time difference calculation unit 110 calculates the logarithmic amplitude at the discrete time t _{n with} respect to the discrete angular frequency ω _k from the amplitudes A _{n, k} (n = n ₁ , n ₁ +1,..., N ₂ , k is an integer). logA _n, which is a time difference of _k logarithmic amplitude time difference Δ _t logA _n, to calculate the _k / _{Δt n} (S110). Logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n is center difference approximation for n satisfying n ₁ <n <n ₂ and forward difference approximation or backward for n = n ₁ or n = n ₂ which is an end point What is necessary is just to calculate by a difference approximation. Of course, for n satisfying n ₁ <n <n ₂ , forward difference approximation or backward difference approximation may be used instead of the center difference approximation. Logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n , for n satisfying n ₁ <n <n ₂ , central difference approximation, n = n ₁ forward difference approximation, n = n ₂ backward difference approximation Can be calculated using, for example, Equation (22), Equation (22b), and Equation (22b).

The time series signal feature calculation unit 120 calculates the group delay characteristic GD _{n, k} at the discrete time t _n and the discrete angular frequency ω _k from the logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n calculated at S110. It is calculated as a feature (S120). Specifically, the standard deviation of the Gaussian window sigma, with an effective width T = N / f _s of the window function

And calculate. For example, using equation (23), the logarithmic amplitude time difference is calculated by the central difference approximation for n satisfying n ₁ <n <n ₂ , the forward difference approximation for n = n _{1 and the} backward difference approximation for n = n ₂ When calculating Δ _t logA _{n, k} / Δt _n

It becomes. In addition, using equation (26), logarithmic amplitude time difference is calculated by the central difference approximation for n satisfying n ₁ <n <n ₂ , the forward difference approximation for n = n _{1 and the} backward difference approximation for n = n ₂ When calculating Δ _t logA _{n, k} / Δt _n

It becomes.

  According to the invention of this embodiment, in the signal analysis based on the short-time Fourier transform, the amplitude information (logarithmic amplitude) that can be stably calculated instead of calculating the estimated value of the group delay characteristic using the phase that is the original definition. And amplitude) can be used to avoid the problem of phase unwrapping uncertainty.

  In addition, since the group delay characteristic can be estimated only from the amplitude information, the complete complex number can be obtained even in the situation where the phase information is missing and only the amplitude information is given out of the complex signal information obtained as a result of the short-time Fourier transform. It is possible to estimate a group delay characteristic equivalent to that when a signal is given.

<Second embodiment>
Hereinafter, the time-series signal feature estimation apparatus 200 will be described with reference to FIGS. As illustrated in FIG. 3, the time-series signal feature estimation apparatus 200 includes a logarithmic amplitude angular frequency difference calculation unit 210, a time-series signal feature calculation unit 220, and a recording unit 190. The recording unit 190 is a component that appropriately records information necessary for processing of the time-series signal feature estimation apparatus 200.

Amplitudes A _{n, k} (n is an integer, k ₁ and k ₂ are integers satisfying k ₂ −k ₁ ≧ 1, k = k ₁ , k ₁ +1,. k ₂ ) shall be obtained in advance.

The time-series signal feature estimation apparatus 200 calculates a time series at discrete time t _n and discrete angular frequency ω _k from amplitude A _{n, k} (n is an integer, k = k ₁ , k ₁ +1,..., K ₂ ). The instantaneous frequency IF _{n, k} that is a signal feature is estimated and output.

The operation of the time-series signal feature estimation apparatus 200 will be described with reference to FIG. The logarithmic amplitude angular frequency difference calculation unit 210 calculates the logarithm at the discrete angular frequency ω _k for the discrete time t _n from the amplitude A _{n, k} (n is an integer, k = k ₁ , k ₁ +1,..., K ₂ ). amplitude logA _n, logarithmic amplitude angular frequency difference is the angular frequency difference between the _{k Δ ω} logA _n, calculates the _k / _{Δω k} (S210). The logarithmic amplitude angular frequency difference Δ _ω logA _{n, k} / Δω _k is the central difference approximation for k satisfying k ₁ <k <k ₂ , or the forward difference approximation for k = k ₁ or k = k ₂ that is an endpoint. What is necessary is just to calculate by reverse difference approximation. Of course, for k satisfying k ₁ <k <k ₂ , forward difference approximation or backward difference approximation may be used instead of the center difference approximation. Logarithmic amplitude angular frequency difference Δ _ω logA _{n, k} / Δω _k , center difference approximation for k satisfying k ₁ <k <k ₂ , forward difference approximation for k = k ₁ , backward difference for k = k ₂ When calculating by approximation, for example, it can be calculated using Expression (29), Expression (29b), and Expression (29b).

The time-series signal feature calculation unit 220 calculates the instantaneous frequency IF _{n, k} at the discrete time t _n and the discrete angular frequency ω _k from the logarithmic amplitude angular frequency difference Δ _ω logA _{n, k} / Δω _k calculated at S210 as a time-series signal. It is calculated as a feature (S220). Specifically, using the standard deviation σ of the Gaussian window

And calculate. For example, using equation _{(30), k 1 <k} < central difference approximation for k satisfying k _2, k = k ₁ forward difference approximation for, for k = k ₂ are logarithmic amplitude angular frequency by the reverse difference approximation When calculating the difference Δ _ω logA _{n, k} / Δω _k

It becomes. Further, using equation _{(33), k 1 <k} < central difference approximation for k satisfying k _2, k = k ₁ forward difference approximation for, for k = k ₂ are logarithmic amplitude angular frequency by the reverse difference approximation When calculating the difference Δ _ω logA _{n, k} / Δω _k

It becomes.

  According to the invention of this embodiment, in the signal analysis by the short-time Fourier transform, instead of calculating the estimated value of the instantaneous frequency using the phase that is the original definition, the amplitude information (logarithmic amplitude and By calculating using (amplitude), it is possible to avoid the problem of uncertainty in the phase unwrapping process.

  In addition, since the instantaneous frequency can be estimated from only the amplitude information, the complete complex signal can be obtained even in the situation where the phase information is missing and only the amplitude information is given out of the complex signal information obtained as a result of the short-time Fourier transform. It is possible to estimate the instantaneous frequency equivalent to the case where is given.

<Third embodiment>
Hereinafter, the time-series signal feature estimation apparatus 300 will be described with reference to FIGS. As illustrated in FIG. 5, the time-series signal feature estimation apparatus 300 includes a group delay characteristic integral value calculation unit 310, a time-series signal feature calculation unit 320, and a recording unit 190. The recording unit 190 is a component that appropriately records information necessary for processing of the time-series signal feature estimation apparatus 300.

The amplitudes A _{n, k} (n ₁ , n ₂ are integers satisfying n ₂ -n ₁ ≧ 1, n = n ₁ , n ₁ +1,..., N ₂ , k, which are inputs to the time-series signal feature estimation apparatus 300 = 0, 1, ..., N-1) is assumed to have been obtained in advance.

The time-series signal feature estimation apparatus 300 has discrete times t _n and discrete _values from amplitudes A _{n, k} (n = n ₁ , n ₁ +1,..., N ₂ , k = 0, 1,..., N−1). A phase φ _{n, k} that is a time-series signal characteristic at the angular frequency ω _k is estimated and output.

The operation of the time-series signal feature estimation apparatus 300 will be described with reference to FIG. The group delay characteristic integral value calculation unit 310 calculates k = 0,1 from the amplitudes A _{n, k} (n = n ₁ , n ₁ +1,..., N ₂ , k = 0, 1,..., N−1). , ..., N-2 calculates the group delay characteristic integral value Σ (-GD) _{n, k} which is the integral value of the sign inversion value of the group delay characteristic from the discrete angular frequency ω _k to ω _{k + 1} for each k of N-2 (S310). For example, the group delay characteristic integral value calculation unit 310 uses the time-series signal feature estimation apparatus 100 to determine the amplitudes A _{n, k} (n = n ₁ , n ₁ +1,..., N ₂ , k = 0, 1, …, N-1), the group delay characteristic GD _{n, k} is calculated for each k of k = 0, 1,…, N-1, and then the group delay characteristic integrated value Σ (-GD) _{n, k} is It is sufficient to calculate.

As another example, the group delay characteristic integral value calculation unit 310 calculates n satisfying n ₁ <n <n ₂ using the following equation:

For n = n ₁ , calculate using the following formula:

For n = n _2, it may be calculated using the following equation.

As another example, the group delay characteristic integral value calculation unit 310 calculates n satisfying n ₁ <n <n ₂ using the following equation:

For n = n ₁ , calculate using the following formula:

For n = n _2, it may be calculated using the following equation.

The time-series signal feature calculation unit 320 sets the initial phase value to φ _n, k — ₀ (k ₀ is an integer between 0 and N−1) and the group delay characteristic integrated value Σ (−GD) _{n, k} calculated in S310. From (k = 0, 1, ..., N-2), the phase φ _{n, k} (k = 0, 1, ..., N-1, except k = k ₀ ) is calculated as a time-series signal feature. (S320). The integer k ₀ and the initial phase value φ _n, k — ₀ may be recorded in advance in the recording unit 190 or may be given as inputs to the time-series signal feature estimation apparatus 300. Specifically, the phase φ _{n, k} (k = 0, 1,..., N−1, except k = k ₀ ) is calculated by the following procedure.

First, k = k ₀ +1, and the initial phase value φ _{n, k_0} and the group delay characteristic integral value Σ (−GD) _{n, k} (k = k ₀ , k ₀ +1,..., N−2) The phase φ _{n, k — 0 + 1} , phase φ _{n, k — 0 + 2} ,..., Phase φ _{n, N−1} are calculated in ascending order from the following formula (formula (43)).

Next, k = k ₀ -1 and using the initial phase value φ _{n, k_0} and the group delay characteristic integral value Σ (-GD) _{n, k} (k = 0, 1, ..., k ₀ -1) From the following equation (Equation (44)), the phase φ _{n, k — 0-1} , phase φ _{n, k — 0-2} ,..., Phase φ _{n, 0} are calculated in descending order.

In order to calculate the phase by the time-series signal feature calculation unit 320, it is necessary to provide an initial value φ _{n, k — 0} of the phase. For each n, the phase φ _{n, 0} when k ₀ = 0 can be set as an initial value. When the signal to be analyzed is a real number, the value of the phase φ _{n, 0} is limited to either 0 or π. The discrete angular frequency index k = 0 estimated from the equation (23) or the equation (26) at the point where the value of the phase φ _{n, 0} changes from 0 to π or from π to 0 according to the change of the discrete time index n. A large spike characteristic change as shown in FIG. 9B is observed in the estimated value of the group delay characteristic in the surrounding discrete angular frequency index. By detecting such a change _, the value of the phase φ _{n, 0} can be switched between 0 and π and an appropriate initial value can be given.

Note that the starting point of the recursive formula used in the time-series signal feature calculation unit 320 is not limited to k ₀ = 0, so that an appropriate value can be obtained as an initial value for an arbitrary k ₀ . It may be the initial value. In this case, as described above, the estimated value of the phase can be obtained not only by calculating the recursive formula for the value of k in ascending order but also by calculating the recursive formula in descending order.

  According to the invention of this embodiment, the conventional method for obtaining the phase from the complex signal obtained as a result of the short-time Fourier transform can calculate only when the phase value is wrapped between -π and π. By using the amplitude information, the phase value in the unwrapped state can be estimated.

  In addition, since the phase can be estimated only from the amplitude information, even in a situation where the phase information is missing and only the amplitude information is given out of the complex signal information obtained as a result of the short-time Fourier transform, the complete complex signal is It is possible to estimate the phase equivalent to the given case.

  Furthermore, by restoring the estimated phase value under the condition where only the amplitude information is given, it is possible to reconstruct a complex signal in the absence of some information.

<Fourth embodiment>
Hereinafter, the time-series signal feature estimation apparatus 400 will be described with reference to FIGS. As illustrated in FIG. 7, the time-series signal feature estimation apparatus 400 includes an instantaneous frequency integral value calculation unit 410, a time-series signal feature calculation unit 420, and a recording unit 190. The recording unit 190 is a component that appropriately records information necessary for processing of the time-series signal feature estimation device 400.

Amplitudes A _{n, k} (n ₀ is an integer, m is a positive integer, n = n ₀ , n ₀ +1,..., N ₀ + m, k ₁ , k ₂ input to the time-series signal feature estimation apparatus 400 Is an integer satisfying k ₂ −k ₁ ≧ 1, and k = k ₁ , k ₁ +1,..., K ₂ ) is obtained in advance.

The time-series signal feature estimation apparatus 400 is discrete from the amplitudes A _{n, k} (n = n ₀ , n ₀ +1,..., N ₀ + m, k = k ₁ , k ₁ +1,..., K ₂ ). time t _n, to estimate the phase phi _{n, k} is a time-series signal, wherein at discrete angular frequency omega _k, and outputs.

The operation of the time-series signal feature estimation apparatus 400 will be described with reference to FIG. The instantaneous frequency integrated value calculation unit 410 calculates _{n from} the amplitudes A _{n, k} (n = n ₀ , n ₀ +1,..., N ₀ + m, k = k ₁ , k ₁ +1,..., K ₂ ). = n ₀ , n ₀ +1, ..., n ₀ + m-1 instantaneous frequency integral value Σ (2πIF), which is the integral value of 2π times the instantaneous frequency from discrete time t _n to t _{n + 1} for each n _{n, k} is calculated (S410). For example, the instantaneous frequency integrated value calculation unit 410 uses the time-series signal feature estimation apparatus 200 to determine the amplitudes A _{n, k} (n = n ₀ , n ₀ +1,..., N ₀ + m, k = k ₁ , k ₁ +1, ..., k ₂ ), the instantaneous frequency IF _{n, k} is calculated for each n of n = n ₀ , n ₀ +1, ..., n ₀ + m, and the frequency integral value Σ (2πIF ) Calculate _{n, k} .

As another example, the instantaneous frequency integral value calculation unit 410 calculates k satisfying k ₁ <k <k ₂ using the following equation:

k = k ₁ is calculated using the following formula:

k = k ₂ may be calculated using the following equation.

As another example, the instantaneous frequency integral value calculation unit 410 calculates k satisfying k ₁ <k <k ₂ using the following equation:

k = k ₁ is calculated using the following formula:

k = k ₂ may be calculated using the following equation.

The time-series signal feature calculation unit 420 includes the initial phase value φ n — ₀ , _k and the instantaneous frequency integration value Σ (2πIF) _{n, k} (n = n ₀ , n ₀ +1,..., N ₀ + m) calculated in S410. −1), phase φ _{n, k} (n = n ₀ +1,..., N ₀ + m) is calculated as a time-series signal feature (S420). The integer n ₀ and the initial phase value φ n — _{0 , k} may be recorded in advance in the recording unit 190 or may be given as inputs to the time-series signal feature estimation apparatus 400. Specifically, the phase φ _{n, k} (n = n ₀ +1,..., N ₀ + m) is calculated by the following procedure.

n = n ₀ +1, using initial phase value φ _{n_0, k} and instantaneous frequency integration value Σ (2πIF) _{n, k} (n = n ₀ , n ₀ +1, ..., n ₀ + m-1) From the following equation (equation (45)), the phase φ _{n — 0 + 1, k} , the phase φ _{n — 0 + 2, k} ,..., The phase φ _{n — 0 + m, k} are calculated in ascending order.

In order to calculate the phase by the time-series signal feature calculation unit 420, it is necessary to give the initial value φ _{n — 0, k} of the phase. For each k, the phase φ _{0, k} when n ₀ = 0 can be set as an initial value. Specifically, an appropriate initial value can be given by assuming that the signal is silent when the discrete time index n = 0.

  According to the invention of this embodiment, the conventional method for obtaining the phase from the complex signal obtained as a result of the short-time Fourier transform can calculate only when the phase value is wrapped between -π and π. By using the amplitude information, the phase value in the unwrapped state can be estimated.

  In addition, since the phase can be estimated only from the amplitude information, even in a situation where the phase information is missing and only the amplitude information is given out of the complex signal information obtained as a result of the short-time Fourier transform, the complete complex signal is It is possible to estimate the phase equivalent to the given case.

  Furthermore, by restoring the estimated phase value under the condition where only the amplitude information is given, it is possible to reconstruct a complex signal in the absence of some information.

<Modification>
The present invention is not limited to the above-described embodiment, and it goes without saying that modifications can be made as appropriate without departing from the spirit of the present invention. The various processes described in the above embodiment may be executed not only in time series according to the order of description, but also in parallel or individually as required by the processing capability of the apparatus that executes the processes or as necessary.

  For example, the difference between the estimated values of the phases of adjacent indexes is calculated by using the recursive equation in the frequency direction such as Equation (43) or Equation (44) in the third embodiment and Equation (45) in the fourth embodiment. ) Or an estimated value obtained by a recursive equation in the time direction, or a weighted average.

<Supplementary note>
The apparatus of the present invention includes, for example, a single hardware entity as an input unit to which a keyboard or the like can be connected, an output unit to which a liquid crystal display or the like can be connected, and a communication device (for example, a communication cable) capable of communicating outside the hardware entity. Can be connected to a communication unit, a CPU (Central Processing Unit, may include a cache memory or a register), a RAM or ROM that is a memory, an external storage device that is a hard disk, and an input unit, an output unit, or a communication unit thereof , A CPU, a RAM, a ROM, and a bus connected so that data can be exchanged between the external storage devices. If necessary, the hardware entity may be provided with a device (drive) that can read and write a recording medium such as a CD-ROM. A physical entity having such hardware resources includes a general-purpose computer.

  The external storage device of the hardware entity stores a program necessary for realizing the above functions and data necessary for processing the program (not limited to the external storage device, for example, reading a program) It may be stored in a ROM that is a dedicated storage device). Data obtained by the processing of these programs is appropriately stored in a RAM or an external storage device.

  In the hardware entity, each program stored in an external storage device (or ROM or the like) and data necessary for processing each program are read into a memory as necessary, and are interpreted and executed by a CPU as appropriate. . As a result, the CPU realizes a predetermined function (respective component requirements expressed as the above-described unit, unit, etc.).

  The present invention is not limited to the above-described embodiment, and can be appropriately changed without departing from the spirit of the present invention. In addition, the processing described in the above embodiment may be executed not only in time series according to the order of description but also in parallel or individually as required by the processing capability of the apparatus that executes the processing. .

  As described above, when the processing functions in the hardware entity (the apparatus of the present invention) described in the above embodiments are realized by a computer, the processing contents of the functions that the hardware entity should have are described by a program. Then, by executing this program on a computer, the processing functions in the hardware entity are realized on the computer.

  The program describing the processing contents can be recorded on a computer-readable recording medium. As the computer-readable recording medium, for example, any recording medium such as a magnetic recording device, an optical disk, a magneto-optical recording medium, and a semiconductor memory may be used. Specifically, for example, as a magnetic recording device, a hard disk device, a flexible disk, a magnetic tape or the like, and as an optical disk, a DVD (Digital Versatile Disc), a DVD-RAM (Random Access Memory), a CD-ROM (Compact Disc Read Only). Memory), CD-R (Recordable) / RW (ReWritable), etc., magneto-optical recording medium, MO (Magneto-Optical disc), etc., semiconductor memory, EEP-ROM (Electronically Erasable and Programmable-Read Only Memory), etc. Can be used.

  The program is distributed by selling, transferring, or lending a portable recording medium such as a DVD or CD-ROM in which the program is recorded. Furthermore, the program may be distributed by storing the program in a storage device of the server computer and transferring the program from the server computer to another computer via a network.

  A computer that executes such a program first stores, for example, a program recorded on a portable recording medium or a program transferred from a server computer in its own storage device. When executing the process, the computer reads a program stored in its own recording medium and executes a process according to the read program. As another execution form of the program, the computer may directly read the program from a portable recording medium and execute processing according to the program, and the program is transferred from the server computer to the computer. Each time, the processing according to the received program may be executed sequentially. Also, the program is not transferred from the server computer to the computer, and the above-described processing is executed by a so-called ASP (Application Service Provider) type service that realizes the processing function only by the execution instruction and result acquisition. It is good. Note that the program in this embodiment includes information that is used for processing by an electronic computer and that conforms to the program (data that is not a direct command to the computer but has a property that defines the processing of the computer).

  In this embodiment, a hardware entity is configured by executing a predetermined program on a computer. However, at least a part of these processing contents may be realized by hardware.

Claims (8)

Discrete time t _n of discrete short-time Fourier transform of time series signal analyzed using Gaussian window with standard deviation σ (σ is a positive real number) as window function, discrete angular frequency ω _k (n is discrete time index) Integer, k is an integer representing the discrete angular frequency index)

(Where f _s is the sampling frequency and N is the discrete Fourier transform score),
_Let A _{n, k be the} amplitude at discrete time t _n and discrete angular frequency ω _k ,
From the amplitude A _{n, k} (n ₁ , n ₂ are integers satisfying n ₂ −n ₁ ≧ 1, n = n ₁ , n ₁ +1,..., N ₂ , k is an integer), the discrete angular frequency ω _k A logarithmic amplitude time difference calculating unit for calculating a logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n that is a time difference of the logarithmic amplitude logA _{n, k} at discrete time t _n ,
The logarithmic amplitude time difference delta _t logA _n, from _k / Delta] t _n, the time-series signal for calculating discrete time t _n using the following formula, the group delay characteristic GD _n in a discrete angular frequency omega _k, a _k as time series signals, wherein With the feature calculator

(Where T is the effective width of the window function, T = N / f _s )
A time-series signal feature estimation device. The time-series signal feature estimation device according to claim 1,
The logarithmic amplitude time difference calculation unit calculates the logarithmic amplitude time difference Δ _t logA _{n, k} / Δt _n for n satisfying n ₁ <n <n ₂ by center difference approximation,
Let δ ₁ and δ ₂ be non-zero numbers,
The time-series signal feature calculation unit calculates the group delay characteristic GD _{n, k} for n satisfying n ₁ <n <n _{2 according} to the following equation:

A time-series signal feature estimation apparatus. Discrete time t _n of discrete short-time Fourier transform of time series signal analyzed using Gaussian window with standard deviation σ (σ is a positive real number) as window function, discrete angular frequency ω _k (n is discrete time index) Integer, k is an integer representing the discrete angular frequency index)

(Where f _s is the sampling frequency and N is the discrete Fourier transform score),
_Let A _{n, k be the} amplitude at discrete time t _n and discrete angular frequency ω _k ,
The amplitude A _{n, k} (n is an integer, k _1, k ₂ is an integer satisfying _{k 2 -k 1 ≧ 1, k} = k 1, k 1 + 1, ..., k 2) from the discrete time t _n A logarithmic amplitude angular frequency difference calculating unit for calculating a logarithmic amplitude angular frequency difference Δ _ω logA _{n, k} / Δω _k that is an angular frequency difference of the logarithmic amplitude logA _{n, k at} the discrete angular frequency ω _k ;
The logarithmic amplitude angular frequency difference delta _omega logA _n, from _k / [Delta] [omega _k, the time-series signal for calculating discrete time t _n using the following equation, the instantaneous frequency IF _n in a discrete angular frequency omega _k, a _k as time series signals, wherein With the feature calculator

A time-series signal feature estimation device. The time-series signal feature estimation device according to claim 3,
The logarithmic amplitude angular frequency difference calculation unit calculates the logarithmic amplitude angular frequency difference Δ _ω logA _{n, k} / Δω _k for k satisfying k ₁ <k <k ₂ by center difference approximation,
Let δ ₄ and δ ₅ be non-zero numbers,
The time-series signal feature calculator calculates the instantaneous frequency IF _{n, k} for k satisfying k ₁ <k <k _{2 according} to the following equation:

A time-series signal feature estimation apparatus. Discrete time t _n of discrete short-time Fourier transform of time series signal analyzed using Gaussian window with standard deviation σ (σ is a positive real number) as window function, discrete angular frequency ω _k (n is discrete time index) Integer, k is an integer representing the discrete angular frequency index)

(Where f _s is the sampling frequency and N is the discrete Fourier transform score),
_Let A _{n, k be the} amplitude at discrete time t _n and discrete angular frequency ω _k ,
The amplitudes A _{n, k} (n ₁ and n ₂ are integers satisfying n ₂ −n ₁ ≧ 1, n = n ₁ , n ₁ +1,..., N ₂ , k = 0, 1,..., N−1 ) To k = 0,1,..., N-2, the group delay characteristic integral value Σ (− that is the integral value of the sign inversion value of the group delay characteristic from the discrete angular frequency ω _k to ω _{k + 1} for each k. GD) a group delay characteristic integral value calculation unit for calculating _{n, k} ;
k ₀ is an integer between 0 and N-1, and the initial phase value is φ _{n, k_0} ,
From the group delay characteristic integral value Σ (−GD) _{n, k} (k = k ₀ , k ₀ +1,..., N−2), the phase φ is ascending from k = k ₀ +1 using the following equation: _{n, k} (k = k ₀ +1, ..., N-1) are calculated as time series signal features,

From the group delay characteristic integral value Σ (−GD) _{n, k} (k = 0, 1,..., K ₀ −1), using the following equation, the phase φ _{n, k in} descending order from k = k ₀ −1 a time series signal feature calculating unit for calculating (k = 0, 1,..., k ₀ -1) as a time series signal feature;

A time-series signal feature estimation device. The time-series signal feature estimation apparatus according to claim 5,
The group delay characteristic integral value calculation unit uses the time-series signal feature estimation device according to claim 1 to calculate the amplitudes A _{n, k} (n = n ₁ , n ₁ +1,..., N ₂ , k = 0, 1,..., N-1), the group delay characteristic GD _{n, k} is calculated for each k of k = 0, 1,..., N−1 _, and from the group delay characteristic GD _{n, k} , k = The group delay characteristic integral value Σ (−GD) _{n, k} is calculated for each k of 0, 1,..., N−2. Discrete time t _n of discrete short-time Fourier transform of time series signal analyzed using Gaussian window with standard deviation σ (σ is a positive real number) as window function, discrete angular frequency ω _k (n is discrete time index) Integer, k is an integer representing the discrete angular frequency index)

(Where f _s is the sampling frequency and N is the discrete Fourier transform score),
_Let A _{n, k be the} amplitude at discrete time t _n and discrete angular frequency ω _k ,
The amplitude A _{n, k} (n ₀ is an integer, m is a positive integer, n = n ₀ , n ₀ +1,..., N ₀ + m, k ₁ and k ₂ satisfy k ₂ −k ₁ ≧ 1. Integer, k = k ₁ , k ₁ +1, ..., k ₂ ) and n = n ₀ , n ₀ +1, ..., n ₀ + m-1 for each n of discrete times t _n to t _{n + 1} An instantaneous frequency calculation unit for calculating an instantaneous frequency integral value Σ (2πIF) _{n, k} that is an integral value of 2π times the instantaneous frequency
The initial value of the phase is φ _{n_0, k}
From the instantaneous frequency integrated value Σ (2πIF) _{n, k} (n = n ₀ , n ₀ +1,..., N ₀ + m−1), using the following formula, phase is increased from n = n ₀ +1 in ascending order a time-series signal feature calculator that calculates φ _{n, k} (n = n ₀ +1,..., n ₀ + m) as a time-series signal feature;

A time-series signal feature estimation device.   A program for causing a computer to function as the time-series signal feature estimation apparatus according to any one of claims 1 to 7.

Images