JP2016090670A - Frequency analysis device and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a frequency analysis device and program that allow time-series digital data such as an acoustic signal and the like to be efficiently analyzed in short time.SOLUTION: A frequency analysis device performing a frequency analysis of a frame of a length having a prescribed frame length N for each prescribed hop length H with respect to input time-series digital data is configured to: add a prescribed phase difference to a result of the frequency analysis of an analyzed m-1-th frame to thereby acquire a first term of an analysis result of an m-th frame; apply a prescribed discrete Fourier transform processing to a difference between sample data on the m-1-th frame and new sample data on the m-th frame to thereby acquire a second term of an analysis result of the m-th frame; and add the second term to the first term to thereby acquire a result of the frequency analysis of the m-th frame, which is an analysis object this time.SELECTED DRAWING: Figure 2

Description

  The present invention relates to a frequency analysis apparatus and program for analyzing frequency components of input time-series digital data.

  Conventionally, Discrete Fourier Transform (DFT) is known as a technique for performing frequency analysis of time-series digital data such as acoustic signals. Normally, it is impossible to perform integration over an infinite interval, so the data of the predetermined frame length is extracted from the digital signal to be analyzed for each predetermined hop length and is subjected to a discrete Fourier transform by applying a window function. A so-called Short-Time Fourier Transform (STFT) technique is used in which processing is repeated and analyzed.

FIG. 7 shows a conceptual diagram of the STFT. Reference numeral 701 is a graph showing a time change of the digital signal x (t) which is time series data. Reference numeral 702 denotes a frame (data string) of a section having a frame length of N that is cut out first. A DFT operation is performed on this frame 702 to obtain the first analysis result X ₀ (f). The subscript of X ₀ (f) indicates a frame number that identifies the analyzed frame (starting with frame 0), and the parameter f indicates a number (or frequency) that identifies the frequency bin. The analysis result X ₀ (f) is a data string representing complex components of the signal to be analyzed in each frequency bin. Reference numeral 703 denotes the next frame 1. Frame 1 is obtained by cutting out a section of frame length N from a position advanced from frame 0 by hop length H. A DFT operation is performed on frame 1 to obtain an analysis result X ₁ (f). Similarly, for each hop length H, the section of the frame length N is sequentially cut out as a frame and the DFT operation is executed. In the DFT calculation, a Fast Fourier Transform (FFT) technique that can perform an operation at high speed is often used.

Furthermore, a technique called sliding discrete Fourier transform (Sliding-DFT: SDFT) using the circularity of DFT has also been proposed (see Non-Patent Documents 1 and 2). FIG. 8 shows a conceptual diagram of SDFT. Reference numeral 801 denotes a graph of the signal x (t) similar to 701. Reference numeral 802 denotes a frame of a section whose frame length is N that is cut out first. A DFT operation is performed on this frame 802 to obtain the first analysis result X ₀ (f). Next, the data to be processed in the frame 802 is slid by one sample. That is, the first sample 803 is deleted and a sample 805 next to the last sample of the frame 802 is added, and a frequency analysis result of a series of data indicated by the shaded areas 804 and 805 is obtained. This calculation can be obtained by correcting the analysis result X ₀ (f) of the frame 802. This correction calculation is performed for each frequency bin, but it is simpler and requires less calculation than the DFT calculation. Similarly, the analysis target data is slid as shown in 806 and 807, and the analysis result of each frequency bin is obtained by correction calculation. By sliding to the hop length H, the analysis result of the next frame 808 is obtained.

IEEE SIGNAL PROCESSING MAGAZINE March 2003, pages 74-80, dsp tips & tricks "The Sliding DFT", Eric Jacobsen and Richard Lyons IEEE SIGNAL PROCESSING MAGAZINE January 2004, pp. 110-111, dsp tips & tricks "An Update to the Sliding DFT", Eric Jacobsen and Richard Lyons

  The calculation of the short-time Fourier transform applied to the frequency analysis of the acoustic signal is required to be performed as efficiently and in a short time as possible.

In general, the calculation amount of FFT is O [N log ₂ (N)]. Based on this, the amount of calculation in each frame in the above-described STFT and SDFT will be considered. First, in the STFT of FIG. 7, since the FFT is performed once, the calculation amount is O [N log ₂ (N)]. In the SDFT of FIG. 8, since the slide process is performed only for the hop length H, the calculation amount is O (NH).

  An object of the present invention is to provide a frequency analyzer and a program that can efficiently perform frequency analysis of time-series digital data such as an acoustic signal in a short time.

  In order to achieve the above object, an invention according to claim 1 is a frequency analyzer for performing frequency analysis of a frame having a predetermined frame length N for each predetermined hop length H with respect to input time-series digital data. Means for obtaining the first term of the analysis result of the mth frame that is the current analysis object by adding a predetermined phase difference to the frequency analysis result of the m−1th frame that has already been analyzed; Means for taking the difference between the sample data of the first frame and the new sample data of the m-th frame that is the object of analysis this time, and applying the predetermined discrete Fourier transform processing to the difference, the current analysis Means for acquiring the second term of the analysis result of the m-th frame that is the object, and frequency analysis of the m-th frame that is the object of analysis this time by adding the second term to the first term Characterized in that it comprises a means for obtaining results.

を施すことにより今回の分析対象であるｍ番目のフレームの分析結果の第１項を取得する手段と、ｍ−１番目のフレームのサンプルデータx(n+mN)と今回の分析対象であるｍ番目のフレームの新たなサンプルデータx(n+mN-N)との差分を

にて求める手段と、前記差分に対して離散フーリエ変換の処理を施すことにより

でｍ番目のフレームの分析結果の第２項を取得する手段と、前記第１項に第２項を加えることにより今回の分析対象であるｍ番目のフレームの周波数分析結果X_m(k)を取得する手段とを備えることを特徴とする。 The invention according to claim 2 is a frequency analysis device that performs frequency analysis of a frame having a predetermined frame length N for each predetermined hop length H with respect to input time-series digital data, and has already been analyzed. -Calculation that adds a predetermined phase difference to the frequency analysis result X _m-1 (k) of the -1st frame

To obtain the first term of the analysis result of the m-th frame that is the analysis target of this time, the sample data x (n + mN) of the (m−1) -th frame, and m of the analysis target of the current analysis The difference from the new sample data x (n + mN-N) in the second frame

And by applying discrete Fourier transform to the difference

Means for obtaining the second term of the analysis result of the m-th frame and adding the second term to the first term to obtain the frequency analysis result X _m (k) of the m-th frame that is the current analysis target. Means for obtaining.

  According to a third aspect of the present invention, in the frequency analyzer according to the first or second aspect, the hop length H is a power of 2, and the discrete Fourier transform process for the difference is a fast Fourier transform process. Features.

  The invention according to claim 4 is the frequency analyzer according to any one of claims 1 to 3, wherein the hop length H and the frame length N satisfy HQ = N (H is a divisor of N). And the discrete Fourier transform processing for the difference is executed by Q times of discrete Fourier transform loop processing.

  The invention according to claim 5 is the frequency analyzer according to any one of claims 1 to 3, wherein the hop length H and the frame length N satisfy HQ = N (H is a divisor of N). In the discrete Fourier transform process for the difference, the discrete Fourier transform process for the H point difference is executed, and an N point discrete Fourier transform is performed by interpolation using the obtained discrete Fourier transform result of the H point. The conversion result is obtained.

  According to a sixth aspect of the present invention, in the frequency analyzer according to any one of the first to fifth aspects, an operation for convolving a frequency response of a predetermined window function with the obtained frequency analysis result of the mth frame. It is characterized by obtaining an optimized frequency analysis result by performing.

  The invention according to claim 7 is a frequency analysis program for realizing a frequency analysis device that performs frequency analysis of a frame having a predetermined frame length N for each predetermined hop length H with respect to input time-series digital data. The first term of the analysis result of the mth frame, which is the current analysis target, is added to the frequency analysis result of the m−1th frame that has already been analyzed by being executed by the processing device. Means for obtaining, means for obtaining a difference between the sample data of the (m-1) th frame and new sample data of the mth frame to be analyzed this time, and processing of a predetermined discrete Fourier transform for the difference Means for obtaining the second term of the analysis result of the mth frame that is the subject of analysis this time, and the present analysis by adding the second term to the first term A frequency analysis program to function as a frequency analyzer and means for obtaining frequency analysis results of the m-th frame is elephant.

  According to the present invention, frequency analysis for time-series digital data can be efficiently performed at a higher speed than before.

The block diagram which shows the hardware constitutions of the frequency analyzer which is embodiment of this invention Conceptual diagram of frequency analysis in this embodiment The figure which shows the numerical formula (the 1) for demonstrating the principle of the calculation of the frequency analysis in this embodiment The figure which shows the numerical formula (the 2) for demonstrating the principle of the calculation of the frequency analysis in this embodiment The figure which shows the numerical formula (the 3) for demonstrating the principle of the calculation of the frequency analysis in this embodiment The flowchart which shows the flow of the frequency analysis process in the frequency analyzer of this embodiment (the 1) The flowchart which shows the flow of the frequency analysis process in the frequency analyzer of this embodiment (the 2) Flowchart showing a modification of frequency analysis processing Conceptual diagram of short-time Fourier transform (STFT) Conceptual diagram of sliding discrete Fourier transform (SDFT)

  Hereinafter, embodiments of the present invention will be described with reference to the drawings.

  FIG. 1 is a block diagram showing a hardware configuration of a frequency analyzer according to an embodiment of the present invention. The CPU 101 is a processing device that controls the operation of the entire apparatus. The memory 102 is a storage device that stores various programs executed by the CPU 101, various data, and the like, and can be configured by appropriately combining RAM, ROM, flash memory, hard disk, and the like. The display 103 is a display for displaying various information provided on the operation panel of the apparatus. The operation element 104 is a variety of operation elements that are provided on the operation panel of the apparatus and are operated by the user. The signal processing unit 105 is a DSP, for example, and executes various signal processing programs based on instructions from the CPU 101. In particular, the signal processing unit 105 performs a frequency analysis process on the acoustic signal input via the waveform I / O 106 by executing a frequency analysis process program of FIGS. 5A and 5B described later, and outputs the analysis result to the CPU 101. . The bus 110 is a bus line that connects these parts, and is a generic term for a control bus, a data bus, and an address bus.

  In the apparatus of FIG. 1, the frequency analysis process is executed by the signal processing apparatus 105, but it may be configured to be executed by the CPU 101. In addition, the frequency analyzer can be realized by causing a general-purpose personal computer (PC) to execute the program of the apparatus of the present embodiment described below. Further, the apparatus of FIG. 1 is configured as a single frequency analysis apparatus. For example, the apparatus of FIG. A frequency analysis function may be realized.

FIG. 2 shows a conceptual diagram of frequency analysis in the apparatus of FIG. 201 is a graph showing temporal changes in the acoustic signal x (t) to be analyzed. Reference numeral 202 denotes a frame 0 that is first cut out as a processing target. For the first frame 0, DFT calculation is performed in the same manner as in the conventional technique 702 in FIG. 7 and 802 in FIG. 8 to obtain the analysis result X ₀ (k).

  Reference numeral 206 denotes the next frame 1 to be processed. In order to obtain the DFT calculation result for the frame 1, the present embodiment focuses on the difference between the previous frame 0 and the current frame 1. In the frames 0 and 1, the waveform samples in the shaded area 204 overlap. Further, if the waveform sample of the shaded 203 is deleted from the frame 0 and the waveform sample of the shaded 205 is added, the frame 1 is obtained. Therefore, in this embodiment, from the term obtained by adding a predetermined phase difference to the analysis result of the previous frame, and the correction term obtained based on the DFT calculation result for the difference between the previous frame and the current analysis target frame, This time, obtain the analysis result of the analysis target frame.

  Hereinafter, the principle of such an analysis method according to the present invention will be described. 3, 4A, and 4B are mathematical expressions for explaining the principle of frequency analysis calculation in the present embodiment.

Expression (1) represents an arithmetic expression of DFT applied to frame 0. In Expression (1), x (n) is an input signal and indicates a waveform sample value (represented by a complex number) at time n. e is the base of the natural logarithm, j is the imaginary unit, and π is the pi. N is a frame length, and k is a number (referred to as a bin number) that identifies a frequency bin. X ₀ (k) of the calculation result is a complex number indicating a component for each frequency bin, and from this, the amplitude and phase of the input signal in each frequency bin can be calculated.

  The following formula (2) shows a DFT calculation formula applied to the frame 1. Although it is an arithmetic expression similar to Expression (1), it is an operation for frame 1 and is an operation that takes a sum (Σ) from N = 0 to N−1. Therefore, the x (n) portion of Expression (1) is Replaced with x (n + H). Equation (3) divides the calculation of the sum of Equation (2) into a first term corresponding to the calculation for the shaded area 204 in FIG. 2 and a second term corresponding to the calculation for the shaded area 205 in FIG. Is. The first term of equation (4) is rewritten from n = 0 to NH-1 to n = H to N-1 in the first term of equation (3), and + H of x (n + H) The part of was adjusted. In the second term of equation (4), the constant of the exponent part of the second term of equation (3) is given before Σ. Formula (5) puts the constant of the exponent part of the 1st term of formula (4) before Σ. Equation (6) takes the form of the sum in the range of n = 0 to N−1 for the first term of Equation (5), and the excess of the sum is subtracted by the second term. In Expression (7), exp (−j2π / N) kN = exp (−j2πk) = 1 is added to the second term of Expression (6). Equation (7) can be transformed into (8) and (9) to finally obtain Equation (10).

The first term of the equation (10) is obtained by adding a phase difference to the DFT result X ₀ (k) of the frame 0 immediately before the frame 1 to be obtained. Since x (n + N) −x (n) included in the second term of Expression (10) is a difference obtained by subtracting the waveform sample value of Frame 0 from the waveform sample value of Frame 1, Expression (10) The second term of () can be said to be a correction term for obtaining the analysis result of frame 1 based on the analysis result of frame 0. This correction term can be said to correspond to deleting the shaded area 203 and adding the shaded area 205 in FIG.

  Equation (11) in FIG. 4A is obtained by extracting the Σ portion of the correction term of the second term of Equation (10). Using Y (k) in equation (11), equation (10) can be transformed into equation (12). Although Y (k) in equation (11) is a portion with a large calculation amount from the form of equation (10), equation (11) is a DFT operation on the difference and can be calculated together, so that the conventional method It can calculate faster.

  Expression (13) is a modification of Expression (11) on the assumption that HQ = N (H is a divisor of N). Since H is a divisor of N, k in equation (11) is replaced with ^ kQ + i (^ k is k with ^ (circumflex or hat symbol) at the top in Fig. 4A. Shall be shown). This equation (13) can be calculated by executing DFT Q times. Therefore, the amount of calculation when the calculation is performed using Expression (13) is O [Q * H * log (H)] = O [N * log (H)].

FIG. 4B shows representations of the (m-1) th frame and the mth frame. Expressions (15) and (16) are DFT arithmetic expressions for the (m-1) th and mth frames. Expression (17) is a modification of Expression (16) and corresponds to Expression (10) in FIG. Expression (18) is a modification of Expression (17) using Y _m (k) defined as in Expression (19). Expression (20) is an expression corresponding to Expression (13) when HQ = N (H is a divisor of N). When calculating the m-th frame based on the analysis result X _m-1 (k) of the _m-1 frame, Equation (18) may be applied. In addition, when HQ = N (H is a divisor of N), it can be calculated at high speed by applying equation (20).

  FIG. 5A and FIG. 5B are flowcharts showing the flow of frequency analysis processing in the frequency analyzer of the present embodiment. This processing program is executed by the signal processing unit 105 in FIG. Note that the processing in this flowchart is processing assuming that HQ = N (H is a divisor of N).

Steps 501 to 504 are parts for performing the calculation of the first frame 0. In step 501, N waveform samples (N is the frame length) of the first frame 0 are taken into the buffer. m is a variable for storing the frame number, and is set to an initial value m = 0 here. In step 502, DFT is performed on the N waveform samples in the buffer to obtain X ₀ (k). In step 503, the calculation result X ₀ (k) is output. The calculation result is output to the CPU 101 and used as a frequency analysis result. In step 504, the waveform sample x (n) of frame 0 and the DFT result X ₀ (k) are temporarily stored for processing of the next frame.

  Steps 505 to 517 are loop processing for obtaining the analysis result of the frame m based on the analysis result of the frame m−1 while stepping m from 1 in order. First, in step 506, a phase difference is added to the value of each frequency bin of the DFT result of the previous frame m−1 that has been temporarily stored. This calculation is the calculation of the first term of the equation (18) in FIG. 4B. In step 507, new H waveform samples of the frame m to be processed this time are acquired. In FIG. 2, this corresponds to the waveform sample in the shaded area 205. Next, in step 508, the difference between the old H waveform samples of the previous frame m-1 (corresponding to 203 in FIG. 2) and the new H waveform samples of frame m (corresponding to 205 in FIG. 2). Take. This corresponds to a portion that takes a difference in Σ in equation (17) of FIG. 4B.

The next steps 509 to 514 are loop processing in which the correction term of the second term of the equation (18) is obtained and sequentially added to the first term while i is stepped in the range of 0 to Q-1. First, in step 510, pre-processing for applying phase adjustment to the difference obtained in step 508 is performed. This is the calculation in brackets of Σ in equation (20) with one value of i. At step 511, DFT is applied to the pre-processed H samples that are the result of step 510. This is the calculation of equation (20). In step 512, the result of step 511 is adjusted. This is the calculation of the second term of Equation (18). In step 513, the result of step 512 is added to the result of adding the phase difference to the DFT result of the previous frame m-1 calculated (ie, the result of step 506). The above is processed while i is incremented, and when it is processed Q times, the loop is exited. The result of step 513 at this time is the analysis result X _m (k) of frame m.

In step 515, the analysis result X _m (k) is output. In step 516, the waveform sample and analysis result of the current frame are temporarily stored for processing of the next frame. In step 517, when the process reaches the last frame, the process exits from the loop and ends.

  In the above embodiment, it has been described that the DFT is calculated in Step 511 and Expression (20). However, if the hop length H is a power of 2, a fast Fourier transform (FFT) method can be applied. In this case, calculation can be performed with Q FFTs, and generally FFT can be performed faster than DFT. Therefore, in the first modification of the above embodiment, the hop length H is a power of 2, and the calculation in step 511 is FFT. Thereby, even faster analysis can be realized.

  A second modification of the above embodiment will be described. Usually, in the calculation of DFT or FFT, the analysis result can be optimized by applying a window function. However, the normal window function is to optimize the data sequence of the frame subjected to the DFT or FFT operation by applying the window function. However, in the method of the above embodiment, the analysis result of the current frame is the previous one. Since it is obtained using the analysis result of the frame, if the window function is multiplied by the previous frame, the value of the previous frame multiplied by the window function is restored when obtaining the difference from the current frame. The calculation becomes complicated because it must be taken into consideration. Therefore, in the present embodiment, the conventional window function method cannot be applied as it is.

In view of the fact that the multiplication of the window function in the time domain corresponds to the convolution operation in the frequency domain, the frequency response of the window function to be applied is obtained in advance, and the obtained result X _m (k) immediately before step 515. Is calculated by convolving the frequency response of the window function. In step 515, the result of the convolution operation is output as an analysis result, and in step 516, the calculation result before the convolution operation is temporarily stored for the next frame. By doing so, it is possible to output a frequency analysis result optimized by applying the window function.

  By selecting an appropriate window function, the amplitude characteristic of the frequency response of the window function can be concentrated on the main lobe and the side lobe can be ignored. Accordingly, if there is a range in which the frequency bins of each frequency bin are desired to be prominent, this can be realized by selecting and applying a window function having a frequency response that emphasizes the range.

  Moreover, although the said convolution may be calculated with respect to all the frequency bins, you may make it obtain | require the said convolution only about the frequency bin of the frequency band which the user is paying attention. For example, the convolution calculation can be omitted for a frequency band that is too low or too high that is unnecessary for the user.

  Further, the frequency response of the window function having different characteristics may be convoluted for each frequency bin.

A third modification of the above embodiment will be described. In the third modification, the above-described steps 509 to 514 are replaced with the processing of steps 601 to 604 in FIG. Other processes are the same as those in FIGS. 5A and 5B. In Steps 509 to 514 of the above embodiment, the difference between the H points is multiplied by W _H ^{(i / Q) n} every time, and the DFT calculation is performed Q times in total. On the other hand, in the present modification, one DFT calculation is performed on the difference between the H points to obtain the DFT result at the H point on the frequency axis, and interpolation is performed using the DFT result at the H point. Find the DFT result of a point.

Specifically, after obtaining the H point difference in step 508, the process proceeds to step 601, and DFT calculation is performed on the obtained H point difference. Next, in step 602, the DFT result is set to Y _m (^ kQ), and Y _m (^ kQ + i) consisting of N points is obtained by interpolating the result. In step 603, the obtained Y _m (^ kQ + i) = Y _m (k) is adjusted. In step 604, the adjusted Y _m (k), which is the result of step 603, is added to the previously calculated DFT result of the previous frame m−1 plus the phase difference. Thus, the analysis result X _m (k) of the frame m can be acquired, and is output in step 515.

  According to the third modified example, a rough DFT result can be obtained by a simple interpolation calculation instead of a strict calculation. For interpolation, a conventionally known interpolation method such as linear interpolation or four-point interpolation may be used.

  DESCRIPTION OF SYMBOLS 101 ... Central processing unit (CPU), 102 ... Memory, 103 ... Display, 104 ... Operator, 105 ... Signal processing part (DSP), 106 ... Waveform I / O.

Claims (7)

A frequency analysis device that performs frequency analysis of a frame having a predetermined frame length N for each predetermined hop length H with respect to input time-series digital data,
Means for obtaining the first term of the analysis result of the mth frame that is the current analysis object by adding a predetermined phase difference to the frequency analysis result of the m−1th frame that has already been analyzed;
means for taking a difference between the sample data of the (m-1) th frame and the new sample data of the mth frame to be analyzed this time;
Means for obtaining a second term of the analysis result of the m-th frame that is the current analysis target by performing a predetermined discrete Fourier transform process on the difference;
Means for obtaining a frequency analysis result of the m-th frame that is the current analysis object by adding the second term to the first term. A frequency analysis device that performs frequency analysis of a frame having a predetermined frame length N for each predetermined hop length H with respect to input time-series digital data,
An operation for adding a predetermined phase difference to the frequency analysis result X _m-1 (k) of the already analyzed m−1 th frame.

Means for obtaining the first term of the analysis result of the m-th frame that is the object of analysis this time,
The difference between the sample data x (n + mN) of the (m-1) th frame and the new sample data x (n + mN-N) of the mth frame to be analyzed this time

Means to find in
By applying discrete Fourier transform to the difference

And means for obtaining the second term of the analysis result of the mth frame,
Means for acquiring a frequency analysis result X _m (k) of the m-th frame that is the current analysis target by adding the second term to the first term. The frequency analyzer according to claim 1 or 2,
The hop length H is a power of 2,
The frequency analysis apparatus, wherein the discrete Fourier transform process for the difference is a fast Fourier transform process. In the frequency analyzer as described in any one of Claim 1 to 3,
The hop length H and the frame length N are set so as to satisfy HQ = N (H is a divisor of N), and the discrete Fourier transform process for the difference is executed by Q discrete loop transform processes. The frequency analyzer characterized by this. In the frequency analyzer as described in any one of Claim 1 to 3,
The hop length H and the frame length N are set so as to satisfy HQ = N (H is a divisor of N), and the process of the discrete Fourier transform for the difference executes the process of the discrete Fourier transform for the difference of the H points. A frequency analyzer that obtains the result of N-point discrete Fourier transform by interpolation using the obtained result of H-point discrete Fourier transform. In the frequency analyzer as described in any one of Claim 1 to 5,
An optimized frequency analysis result is obtained by performing an operation of convolving a frequency response of a predetermined window function on the obtained frequency analysis result of the m-th frame. A frequency analysis program for realizing a frequency analysis device for performing frequency analysis of a frame having a predetermined frame length N for each predetermined hop length H with respect to input time-series digital data,
By running it on the processing device,
Means for obtaining the first term of the analysis result of the mth frame that is the current analysis object by adding a predetermined phase difference to the frequency analysis result of the m−1th frame that has already been analyzed;
means for taking a difference between the sample data of the (m-1) th frame and the new sample data of the mth frame to be analyzed this time;
Means for obtaining a second term of the analysis result of the m-th frame that is the current analysis target by performing a predetermined discrete Fourier transform process on the difference;
A frequency analysis program that functions as a frequency analysis device comprising: adding a second term to the first term to obtain a frequency analysis result of an m-th frame that is a current analysis target.

Images