Time domain filter device, measured seismic intensity approximation device, measured seismic intensity approximation system, and time domain filter method

Abstract

To provide a time domain filter device, and time domain filter method that enable performing filter processing as it is without interpolation even in a low sampling frequency, and to provide an instrumental seismic intensity estimation device and instrumental seismic intensity estimation method that have a small estimation error while keeping immediacy, using the time domain filter device.SOLUTION: A time domain filter device, which performs filter processing of a time series of ground motion acceleration, comprises: a first filter that performs filter processing in an expression (11); a second filter that performs filter processing in an expression (12); a third filter that performs filter processing in an expression (13); a fourth filter that performs filter processing in an expression (14); an output series acquisition unit that acts these filter on an input series y(k) to acquire an output series y(k) in an expression (15); a gain adjustment unit that adjusts gains in an expression (16); and a computation unit that implements a computation 11 to a computation 17.SELECTED DRAWING: Figure 9

Description

[Technical field]

[0001] The present invention relates to a time domain filtering device and a time domain filtering method for filtering a time series of earthquake acceleration, and a seismic intensity estimating device and a seismic intensity estimating system for estimating seismic intensity from a time series of earthquake acceleration filtered using a time domain filtering device. [Background technology]

[0002] The current situation is that the measured seismic intensity indicated by the Japan Meteorological Agency in the Official Gazette No. 1831, Japan Meteorological Agency Notification No. 4, Notification Concerning the Enforcement Regulations of the Meteorological Services Act, has been used as an index for issuing warnings and controlling railways, etc., when an earthquake occurs. In order to calculate the measured seismic intensity, it is necessary to accumulate earthquake records for a certain period of time and then filter them using FFT, so it is not possible to calculate values ​​for earthquake records that are measured every moment, and it is not suitable for applications that require immediacy.

[0003] The measured seismic intensity is currently the most widely recognized strong motion index, but it is not suitable for applications that require rapid reporting because it requires frequency domain filtering for calculation. To solve this problem, the present inventor proposed a real-time calculation method in Patent Document 1 that obtains an approximation of the seismic intensity in real time by substituting a time domain approximation filter for the filter calculation in the frequency domain. In addition, the present applicant proposed a real-time calculation method in Patent Document 2 that maintains immediacy while further reducing conversion errors by improving the approximation filter in Patent Document 1.

[0004] The filter described in Patent Document 2 is configured so that the difference from the filter based on the Japan Meteorological Agency Notification is within 0.974 to 1.029 times in the main band of 0.1 Hz to 10 Hz. This can be read from the solid line C, which shows the relative amplitude-frequency characteristics. In this way, the filter described in Patent Document 2 was able to estimate the measured seismic intensity quickly and with a small conversion error. [Prior art documents] [Patent documents]

[0005] [Patent Document 1] Patent No. 4229337 [Patent Document 2] Patent No. 5946067 Summary of the Invention [Problem to be solved by the invention]

[0006] However, conventional estimation of seismic intensity was not able to construct a stable filter at low sampling frequencies (below about 77 Hz), so it was necessary to interpolate to a sampling frequency of 100 Hz and then use an approximation filter.

[0007] The present invention aims to provide a time domain filter device and a time domain filtering method that can process signals directly without interpolation even at low sampling frequencies, and to provide a seismic intensity measurement device and a seismic intensity measurement system that use a time domain filter device to maintain real-timeness and reduce conversion errors. [Means for solving the problem]

[0008] In order to solve the above problems, the time domain filter device of the present invention comprises: 1. A time domain filter apparatus for filtering a time series of ground acceleration, comprising: In the following formula (11), a first filter for filtering is In the following formula (12), a second filter for filtering is In the following formula (13), a third filter for filtering is provided, In the following formula (14), a fourth filter for filtering is In the following formula (15), an output time series acquisition unit applies these filters to the input time series x(k) to acquire an output time series y(k); A gain adjustment unit that adjusts the gain in the following equation (16); A calculation unit that executes calculations 11 to 17; Equipped with. · 1st filter a 0 =8 / (ΔT) 2 +(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 (11) a 1 =2h a1 oh a2 -16 / (ΔT) 2 a 2 =8 / (ΔT) 2 -(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 b 0 =4 / (ΔT) 2 +2h a2 / (ΔT) b 1 =-8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -2h a2 / (ΔT) oh a1 =2πf a1 , oh a2 =2πf a2 · 2nd filter a 0 =16 / (ΔT) 2 +17h a3 / (ΔT)+ ω a3 2 (12) a 1 =2h a3 2 -32 / (ΔT) 2 a 2 =16 / (ΔT) 2 -17h a3 / (ΔT)+ ω a3 2 b 0 =4 / (ΔT) 2 +8.5h a3 / (ΔT) + ω a3 2 b 1 =2h a3 2 -8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -8.5h a3 / (ΔT) + ω a3 2 oh a3 =2πf a3 · 3rd filter a 0 =1+ h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 (13) a 1 =-2+ω b 2 (ΔT) 2 (1-2c) b2 ) a 2 =1- h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 b 0 =1+ h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 b 1 =-2+ω b 2 (ΔT) 2 (1-2c) b1 ) b 2 =1- h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 ω b =2πf b Fourth filter α 0 =1+ h c ω c ΔT +γ c2 ω c 2 (ΔT) 2 (14) α 1 =-2+ω c 2 (ΔT) 2 (1-2γ c2 ) α 2 =1- h c ω c ΔT +γ c2 ω c 2 (ΔT) 2 β 0 =γ c1 ω c 2 (ΔT) 2 β 1 =ω c 2 (ΔT) 2 (1-2γ c1 ) β 2 =γ c1 ω c 2 (ΔT) 2 ω c =2πf c Output time series acquisition section y(k)=[-α 1 y(k-1) -α 2 y(k-2)+β 0 x(k) +β 1 x(k-1) +β 2 x(k-2)] / α 0 (15) however, x(k) is the input time series of the filter. y(k) is the output time series of the filter, k is the number of time steps. ΔT is the sampling interval. π is the ratio of the circumference of a circle to its circumference. Let us assume that. Gain adjustment section y(k)=g d x(k) (16) However, g d is the gain. Operation 11:f a1 = f 0 , f a2 = f 1 Then, the calculation of equations (11) and (15) Operation 12:f a3 = f 1 Then, the calculation of equations (12) and (15) Operation 13:h b1 = h 2a , h b2 = h 2b , f b = f 2 , γ b1 = γ 2a , γ b2 = γ 2b Then, the calculation of equations (13) and (15) Operation 14:h c = h 3 , f c = f 3 , γ c1 = γ 3a , γ c2 = γ 3b Then, the calculation of equations (14) and (15) Operation 15:h c = h 4 , f c = f 4 , γ c1 = γ 4a , γ c2 = γ 4b Then, the calculation of equations (14) and (15) Operation 16:h c = h 5 , f c = f 5 , γ c1 = γ 5a , γ c2 = γ 5bThen, the calculation of equations (14) and (15) Operation 17:g d = g, and calculate equation (16). however, f 0 is the parameter that defines the cutoff frequency of the filter, f 1 is a parameter that specifies the position of the frequency domain where the signal is attenuated at 0.5 order. f 2 , h 2a , h 2b , are parameters that define the characteristics of the correction filter, f 3 , h 3 are the parameters that define the cutoff frequency and damping of the filter, f 4 , h 4 are the parameters that define the cutoff frequency and damping of the filter, f 5 , h 5 are the parameters that define the cutoff frequency and damping of the filter, g is the gain adjustment parameter, Gamma 2a、 Gamma 2b、 Gamma 3a、 Gamma 3b, Gamma 4a、 Gamma 4b、 Gamma 5a、 Gamma 5b is the approximate method for the quadratic integral of the filter. The parameters to be specified, and The following equations (21) to (24) are satisfied. (1 / 4-1 / (2πf 2 ΔT) 2 )≦γ 2b (twenty one) (1 / 4-1 / (2πf 3 ΔT) 2 )≦γ 3b (twenty two) (1 / 4-1 / (2πf 4 ΔT) 2 )≦γ 4b (twenty three) (1 / 4-1 / (2πf 5ΔT) 2 )≦γ 5b (twenty four) Effect of the Invention

[0009] With such a time domain filter device and time domain filter method, even low sampling frequencies can be processed directly without interpolation. Also, with the seismic intensity estimation device and seismic intensity estimation system, the time domain filter device can be used to reduce conversion errors while maintaining immediacy. [Brief description of the drawings]

[0010] [Figure 1] FIG. 1 is a block diagram of a seismic intensity approximation device according to an embodiment of the present invention. [Diagram 2] FIG. 2 is a block diagram of a measured seismic intensity approximation unit according to the present embodiment. [Diagram 3] FIG. 2 is a flowchart for calculating an approximate measured seismic intensity value according to the present embodiment. [Figure 4] FIG. 2 is a diagram illustrating discretization according to the present embodiment. [Diagram 5] FIG. 13 is a diagram illustrating the concept of counting a duration according to the present embodiment. [Figure 6] FIG. 4 is a flowchart of a duration count according to the present embodiment. [Figure 7] FIG. 11 is a diagram showing a second concept of duration counting. [Figure 8] FIG. 11 is a flow chart of a second duration count. [Figure 9] FIG. 1 is a block diagram of a time domain filter device according to an embodiment of the present invention. [Figure 10] FIG. 1 is a diagram showing an earthquake waveform, a waveform obtained by a conventional interpolation method, and a waveform obtained by using the seismic intensity approximation device of this embodiment. [Figure 11] FIG. 1 is a diagram showing a seismic intensity estimation system according to an embodiment of the present invention. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0011] An embodiment of the present invention will be described with reference to the drawings. Fig. 1 is a diagram showing the main configuration of a measured seismic intensity approximation device according to an embodiment of the present invention. In the figure, 1 is a measured seismic intensity approximation device, 2 is a processing unit, 3 is a control / calculation unit, 4 is a measured seismic intensity calculation unit, 5 is an SI value calculation unit, 6 is a response value calculation unit, 7 is a measured seismic intensity approximation unit, 8 is a power supply unit, and 9 is a measurement device.

[0012] The measured seismic intensity approximation device 1 estimates seismic intensity by inputting acceleration and the like measured by a measuring device 9 to a control and calculation unit 3. The processing unit 2 has a control and calculation unit 3 driven by a power supply unit 8. In order to achieve higher accuracy, the control and calculation unit 3 of this embodiment is configured by adding a measured seismic intensity approximation unit 7 to a current measuring seismic intensity meter, but in cases where greater importance is placed on immediacy, the measured seismic intensity calculation unit 4, SI value calculation unit 5, and response value calculation unit 6 may not be included. Furthermore, the measuring device 9 may be provided within the processing unit 2.

[0013] 2 is a diagram showing the main components of the measured seismic intensity approximation unit 7 according to the embodiment of the present invention. In the figure, 71 is a ground acceleration time series acquisition unit, 72 is a time domain filter device, 73 is a calculation unit, 74 is an amplitude time series conversion unit, 75 is a seismic intensity converted amplitude time series conversion unit, 76 is a discretization unit, and 77 is a duration counting unit.

[0014] The ground acceleration time series acquisition unit 71 acquires the time series of ground acceleration of three orthogonal components and sends it to the time domain filter device 72. The time domain filter device 72 filters the time series of ground acceleration acquired by the ground acceleration time series acquisition unit 71. The calculation unit 73 calculates an approximate value of the measured seismic intensity from the time series of ground acceleration filtered by the time domain filter device 72.

[0015] The calculation unit 73 has an amplitude time series conversion unit 74, a seismic intensity converted amplitude time series conversion unit 75, a discretization unit 76, and a duration counting unit 77. The amplitude time series conversion unit 74 converts the time series of ground acceleration filtered by the time domain filter device 72 into an amplitude time series, and the seismic intensity converted amplitude time series conversion unit 75 converts the amplitude time series into a seismic intensity converted amplitude time series. The discretization unit 76 discretizes the seismic intensity converted amplitude time series obtained by the seismic intensity converted amplitude time series conversion unit, and the duration counting unit 77 executes a duration count for each discretized seismic intensity converted amplitude value obtained by the discretization unit.

[0016] Fig. 3 is a flow chart for determining an approximate measured seismic intensity value according to this embodiment. Fig. 4 is a diagram showing discretization according to this embodiment.

[0017] Next, a method for estimating measured seismic intensity by the seismic intensity estimating device 1 will be described with reference to the flowchart shown in Fig. 3. First, in step 1, earthquake motion is observed (ST1).

[0018] Next, in step 2, the three orthogonal components of the ground acceleration time series xa[k], ya[k], za[k] are obtained, sampled at time intervals dt (for example, 0.01 second intervals). Here, k is the number of time steps (ST2). The input required for this method is the time series of ground acceleration, but it does not necessarily have to be obtained by an accelerometer, as long as it can be converted to ground acceleration by performing characteristic correction of the seismometer and differential calculation of the velocity. In addition, the observation of earthquake motion is performed for two horizontal components in the east-west direction and the north-south direction, and a vertical component, respectively, to obtain a time series of three components of ground acceleration. However, this method can also be applied to the time series of any two components or one component of ground acceleration. If the observation of earthquake motion is not performed for two horizontal components and a vertical component, the necessary conversion is performed to obtain a time series of three components of ground acceleration. This conversion may be performed after the filtering process in step 3 described below.

[0019] Next, in step 3, the three orthogonal components of the ground acceleration time series (xa[k], ya[k], za[k]) are subjected to a filter process defined in operation A described later to obtain the three-component filtered ground acceleration time series (xf[k], yf[k], zf[k]) (ST3). The characteristics of the filter will be described later. If the number of components in the earthquake observation is less than three, the non-existent (xf[k], yf[k], zf[k]) can be set to 0 and the processes from step 3 onwards can be carried out. For example, since the order of x, y, and z is irrelevant, if the number of components in the earthquake observation is two, one of xf[k], yf[k], and zf[k] is set to 0, and if the number of components in the earthquake observation is one, any two of xf[k], yf[k], and zf[k] is set to 0.

[0020] Then, in step 4, we use the filtered three-component ground acceleration time series as af[k]=(xf[k]×xf[k]+yf[k]×yf[k]+zf[k]×zf[k]) 1 / 2 (A) The amplitude time series af[k] is then converted into the amplitude time series af[k] (ST4).

[0021] According to the original definition in the Japan Meteorological Agency Notification No. 4, the maximum amplitude level a whose integrated duration is 0.3 seconds or more within a certain time interval is 0 Calculate the conversion formula I = 2 × log(a 0 ) + 0.94 to convert to the equivalent measured seismic intensity, but for the sake of efficiency, the accumulated duration is processed after converting to the equivalent measured seismic intensity.

[0022] First, in step 5, the amplitude time series obtained in step 4 is If[k]=2×log(af[k])+0.94 ···(B) The seismic intensity converted amplitude time series If[k] is converted into seismic intensity using the following conversion formula (ST5).

[0023] Next, in step 6, If[k] is discretized with a discretization width dI of the required accuracy (ST6), as shown in Fig. 4. In this embodiment, since the measured seismic intensity is calculated in increments of 0.1, dI is discretized with 0.01, which is one-tenth of 0.1, in order to maintain sufficient accuracy in the final result.

[0024] Specifically, the value obtained by subtracting the lower limit value Ii (-8 in this embodiment) from If[k] is divided by the discretization width dI to convert it to an integer. If the integer converted value is negative, it is replaced with 0. The time series of discretized seismic intensity converted amplitude values ​​obtained in this way is called Id[k]. Next, the value of Id[k] is limited to be less than the upper limit value N. In this embodiment, N=1800, and if the value of Id[k] is N or more, it is replaced with N-1. As a result, the value of Id[k] ranges from 0 to N-1. In this embodiment, the seismic intensity converted amplitude value after discretization ranges from 0 to 1799, which means that the value before discretization was limited to values ​​equivalent to -8 to 9.99.

[0025] The limit range of values ​​must include the range of the measured seismic intensity to be obtained. The upper limit of the range must be sufficiently larger than the measured seismic intensity of 6.5, which is the lower limit of seismic intensity 7. In the embodiment, it is set to 9.99. This corresponds to the case where the measured acceleration exceeds 10,000 gal, and is a value that exceeds the measurement upper limit of the seismic wave sensor itself, so it is sufficiently large. The lower limit must be set to a value smaller than the measured seismic intensity of 0.4, which is the upper limit of seismic intensity 0. In the embodiment, it is set to -8. This value is set as a value smaller than the measured seismic intensity calculated by regarding the weak shaking of the ground when there is no earthquake as seismic motion, and is a sufficiently small value.

[0026] Next, in step 7, the duration of each discrete seismic intensity converted amplitude value is counted (ST7).

[0027] Here, the basic idea of ​​the first duration counting method will be explained using Fig. 5. Fig. 5 shows the relationship between the number of time steps k and the discrete seismic intensity converted amplitude value, where Fig. 5(a) shows the case where the count range of the integrated duration is set to the number of time steps 2 to 4, and Fig. 5(b) shows the case where the count range of the integrated duration is set to the number of time steps 3 to 5. At this time, the discrete seismic intensity converted amplitude value corresponding to the number of time steps 2 falls outside the count range, and the discrete seismic intensity converted amplitude value corresponding to the number of time steps 5 is added to the count range. Therefore, the discrete seismic intensity converted amplitude value for the number of time steps 3 to 4 does not change even if the number of time steps is increased by 1, so it is left as it is.

[0028] Using this concept, the duration count is executed as a subroutine. Note that the duration mentioned here is clear from the original definition in the Japan Meteorological Agency Notification No. 4, and does not mean the time during which the time series of the discrete seismic intensity converted amplitude value continues to be equal to a certain amplitude, but the time during which the amplitude remains equal to or exceeds a certain amplitude.

[0029] Figure 6 is a flowchart of a subroutine showing the first duration counting method. A memory area for counting duration is secured in advance for each of N discrete seismic intensity converted amplitude values. This array is initialized to 0, and this memory area is called Ic[i]. The subscript i of the array ranges from 0 to N-1.

[0030] First, in step 101, the duration is counted using the above Id[k] and Ic[i] by incrementing all values ​​of the count memory area Ic[i] having an index i less than or equal to the value of Id[k] by +1 (ST101).

[0031] Next, in step 102, it is necessary to reset the counted values ​​for those outside the time interval range for counting the integrated duration. This can be efficiently performed by decrementing all values ​​of Ic[i] with subscript i less than or equal to the value of Id[km] by 1 (ST102). Here, m is the number of samples corresponding to the time interval range for counting the integrated duration. In this embodiment, the time interval range for counting the integrated duration is 10 seconds with 100 Hz sampling. In this case, m is 1000. To perform this, it is necessary to secure m memory areas and constantly store m values ​​from Id[k] to Id[km]. This counting method makes it possible to avoid processing m durations at each time step.

[0032] Next, in step 103, if the maximum subscript of Ic[i] that is equal to or greater than the total count number M (30 in the case of 100 Hz sampling as in the embodiment) corresponding to an integration duration of 0.3 seconds is taken as Ic[p], p is found (ST103) and the subroutine is terminated. In this way, p is found by executing the subroutine, and the process returns to the flowchart of FIG.

[0033] The basic concept of the second duration counting method will be explained using Figure 7. Figure 7(a) shows the seismic intensity converted amplitude values ​​in the time interval for counting the integrated duration rearranged in descending order of value, with a sampling frequency of 10 Hz and a time interval for counting the integrated duration of 1 second. Note that the numbers on the top row in the figure indicate the value of the seismic intensity converted amplitude value, and the numbers on the bottom row indicate the number of time steps at which that value was obtained. In this case, the seismic intensity converted amplitude value for which the integrated duration is 0.3 x 10 Hz = 3 is the third largest value, 6.

[0034] Next, as shown in Fig. 7(b), when the time step advances by one, the discretized seismic intensity converted amplitude value with time step number 0 in Fig. 7(a) moves out of the time interval for counting the integrated duration. Then, as shown in Fig. 7(c), the discretized seismic intensity converted amplitude value with time step number 10 newly enters the time interval for counting the integrated duration. In this case, the seismic intensity converted amplitude value with an integrated duration of 0.3 × 10 Hz = 3 is the third largest value from the top, 5.

[0035] FIG. 8 is a flowchart of a subroutine showing the second duration counting method.

[0036] First, in step 201, it is determined whether the array Is[j] has been initialized (ST201).

[0037] If it is determined in step 201 that initialization has not been performed, initialization is performed in step 202 (ST202). The initialization will be described. First, an array Is[j] having m memory areas is prepared. The subscript j of the array ranges from 1 to m. Here, m is the number of samples corresponding to the time interval range for counting the integrated duration, as in the first method used in step 102 of FIG. 6. In the case of the second method, the sampling is 100 Hz, and the time interval range is 10 seconds, so it is 1000. To initialize the array Is[j], go back m steps from the time step one step before the latest time step (the number of time steps is k), and copy the values ​​Id[km] to Id[k-1] of the discretized seismic intensity converted amplitude value to the array Is[j]. Next, Is[j] is sorted in descending order of value. At this time, prepare m arrays It[n] (the array subscript n is from 1 to m) and record which element of the sorted Is[j] is the value of which time step. This completes the initialization.

[0038] Next, for time step number k, operations are performed only on values ​​Id[km] that fall outside the time interval range for counting the integrated duration and on newly added values ​​Id[k], and Is[j] is kept in an array sorted in descending order of value.

[0039] Specifically, if it is determined in step 201 that initialization has been performed, or if initialization has been performed in step 202, then in step 203, q for which It[q]=km is found, and Is[q] corresponding to Id[km], a value outside the time interval range, is deleted from array Is[j]. The deleted portion is filled up to make array Is[j] an array of m-1 elements. Similarly, It[q] is deleted from array It[n], and the deleted portion is filled up (ST203).

[0040] Next, in step 204, the value of Id[k] is compared with the value of array Is[j], a location is found where Is[j] is sorted in descending order, and an element is inserted there, making Is[j] an array of m elements again. An element is also inserted into the same location of array It[n], with its value set to k, and the number of time steps corresponding to the added value is stored (ST204).

[0041] When the above operations are completed, Is[j] will be an array sorted in ascending order, and It[n] will also remember which element of the updated Is[j] corresponds to which time step.

[0042] Then, in step 205, the discretized seismic intensity converted amplitude value Is[M] corresponding to the integrated duration of 0.3 seconds, which is the Mth largest value, is obtained (ST205), and the subroutine ends. Is[M] corresponds to p in the first method used in step 103 of FIG. 6. Also, M is 30 in the second method, as in the first method. In this way, p is obtained by executing the subroutine, and the process returns to the flowchart of FIG. 3.

[0043] In the case of the second duration counting method, the counting process is based on a comparison of the magnitude of the discrete seismic intensity converted amplitude value. Therefore, the process can be performed using the seismic intensity converted amplitude value itself, rather than the integer-valued discrete seismic intensity converted amplitude value.

[0044] Furthermore, the method of counting the duration may be other than the first and second methods as long as the same effect can be obtained. Furthermore, the data does not have to be stored in an array.

[0045] After the subroutine is completed, the process returns to the flowchart in Fig. 3, and finally, in step 8, an approximate value of the measured seismic intensity is calculated (ST8). The value of p multiplied by the discretization width dI and added to the lower limit value Ii becomes the approximate value Ia[k] of the measured seismic intensity at that time.

[0046] The above process is performed for each time step. Therefore, each time the three components of ground acceleration (xa[k], ya[k], za[k]) are obtained, the approximate seismic intensity Ia[k] at that time can be obtained. As shown in Figure 9, the approximate seismic intensity Ia[k] obtained by this method is obtained as a time series that changes over time while the earthquake motion continues. The maximum value Ix of this time-varying value is defined as the approximate measured seismic intensity in the original definition, where one value is determined for the section where the earthquake motion continues. If it is difficult to define the section where the earthquake motion continues, it can be defined as the maximum value for one minute after the earthquake detection, for example. Since Ia[k] changes over time in this way, it is also possible to issue an immediate warning if the alert level is exceeded before reaching the maximum value Ix.

[0047] The order of the steps described above is not limited to this, and the order of the steps may be changed as long as the estimated measured seismic intensity value does not change.

[0048] Next, we describe a method for designing a time domain filter that approximates the filter defined in the original seismic intensity definition announced by the Japan Meteorological Agency.

[0049] The filter proposed in Patent Document 2 was a series combination of eight filters expressed by the Laplace transform below and a gain adjustment that multiplies the final output by g.

[0050] 1 / (ω 0 s -1 +1) (1') However, ω 0 =2πf 0 (ω 1 s -1 +1) / (ω 1 s -1 +2) (2') (4h 1 s -1 +1) / (ω 1 s -1 +8) (3') (0.25h 1 s -1 +1) / (ω 1 s -1 +0.5) (4') However, oh 1 =2πf 1 (1+2h 2a oh 2 s -1 +oh 2 2 s -2 ) / (1+2h 2b oh 2 s -1 +oh 2 2 s -2 ) (5') However, oh 2 =2πf 2 oh 3 2 s -2 / (1+2h 3 oh 3 s -1 +oh 3 2 s -2 ) (6') However, oh 3 =2πf 3 oh 4 2 s -2 / (1+2h 4 oh 4 s -1 +oh 4 2 s -2 ) (7') However, oh 4 =2πf 4 oh 5 2 s -2 / (1+2h 5 oh 5 s -1 +oh 5 2 s -2) (8') However, ω 5 =2πf 5 It is.

[0051] In this case, the parameters that define the filter characteristics are (f 0 ,f 1 ,f 2 ,f 3 ,f 4 ,f 5 ,h 2a ,h 2b ,h 3 ,h 4 ,h 5 These values ​​were determined based on the degree of agreement with the original amplitude characteristics and the final calculation accuracy, and the combination of parameters with the best degree of approximation was (f 0 =0.45[Hz],f 1 =7.0[Hz],f 2 =0.5[Hz],f 3 =12.0[Hz],f 4 =20.0[Hz],f 5 =30.0[Hz] ,h 2a =1.0,h 2b =0.75,h 3 =0.9,h 4 =0.6,h 5 = 0.6, g = 1.262) was obtained.

[0052] Given the Laplace transform representation of a filter, a digital filter can be constructed using the z-transform. In Patent Document 2, where ΔT is the sampling interval and z is the delay operator, s -1 =(ΔT / 2)(1+ z -1 ) / (1- z -1 ) (30) and s -2 =(ΔT / 12 1 / 2 ) 2 (1+ 10z -1 +z -2 ) / (1- z -1 ) 2 (31) Using this z-transformation, digital filters expressed by the following equations (11') to (14') were obtained.

[0053] Filter 1 α 0 =8 / (ΔT) 2 +(4ω a1 +2ω a2 ) / (ΔT)+ω a1 ω a2 (11') α 1 =2ω a1 ω a2 -16 / (ΔT) 2 α 2 =8 / (ΔT) 2 -(4ω a1 +2ω a2 ) / (ΔT)+ω a1 ω a2 β 0 =4 / (ΔT) 2 +2ω a2 / (ΔT) β 1 =-8 / (ΔT) 2 β 2 =4 / (ΔT) 2 -2ω a2 / (ΔT) ω a1 =2πf a1 , ω a2 =2πf a2

[0054] Filter 2 α 0 =16 / (ΔT) 2 +17ω a3 / (ΔT) + ω a3 2 (12') α 1 =2ω a3 2 -32 / (ΔT) 2 α 2 =16 / (ΔT) 2 -17ω a3 / (ΔT) + ω a32 b 0 =4 / (ΔT) 2 +8.5h a3 / (ΔT) + ω a3 2 b 1 =2h a3 2 -8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -8.5h a3 / (ΔT) + ω a3 2 oh a3 =2πf a3

[0055] · filter3 a 0 =12 / (ΔT) 2 +12h b2 oh b / (ΔT)+ω b 2 (13') a 1 =10h b 2 -24 / (ΔT) 2 a 2 =12 / (ΔT) 2 -12h b2 oh b / (ΔT)+ω b 2 b 0 =12 / (ΔT) 2 +12h b1 oh b / (ΔT)+ω b 2 b 1 =10h b 2 -24 / (ΔT) 2 b 2 =12 / (ΔT) 2 -12h b1 oh b / (ΔT)+ωb 2 ω b =2πf b

[0056] Filter 4 α 0 =12 / (ΔT) 2 +12h c ω c / (ΔT)+ω c 2 (14') α 1 =10ω c 2 -24 / (ΔT) 2 α 2 =12 / (ΔT) 2 -12h c ω c / (ΔT)+ω c 2 β 0 =ω c 2 β 1 =10ω c 2 β 2 =ω c 2 ω c =2πf c

[0057] To obtain the output time series y(k) by applying these filters to the input time series x(k), the following equation (15') may be calculated. y(k)=[-α 1 y(k-1) -α 2 y(k-2)+β 0 x(k) +β 1 x(k-1) +β 2 x(k-2)] / α 0 (15') however, x(k) is the input time series of the filter. y(k) is the output time series of the filter, k is the number of time steps. ΔT is the sampling interval. π is the ratio of the circumference of a circle to its circumference. Let us assume that.

[0058] Also, the gain adjustment is g d Then, the input time series is x(k) and the output time series is y(k), and the following equation (16') can be calculated. y(k)=g d x(k) (16')

[0059] To summarize the above, the approximation filter calculation described in Patent Document 2 requires performing the following seven digital filter calculations in series on the acceleration record. Operation 1': f a1 =f 0 , f a2 =f 1 Then, the calculation of equations (11') and (15') Operation 2': f a3 =f 1 Then, the calculation of equations (12') and (15') Operation 3': h b1 = h 2a , h b2 = h 2b , f b =f 2 Then, the calculation of equations (13') and (15') Operation 4': h c = h 3 , f c =f 3 Then, the calculation of equations (14') and (15') Operation 5': h c = h 4 , f c =f 4 Then, the calculation of equations (14') and (15') Operation 6': h c = h 5 , f c =f 5 Then, the calculation of equations (14') and (15') Operation 7': g d = g, and calculate equation (16').

[0060] In Patent Document 2, in order to be useful for practical earthquake disaster prevention, a filter was designed and its accuracy was verified assuming that the sampling interval ΔT of the ground acceleration time series is 0.01 seconds, which is the standard for seismic intensity meters actually in operation. On the other hand, for research purposes, etc., ΔT longer than the standard 0.01 seconds may be used. In this case, ΔT and f 2 ,f 3 ,f 4 ,f 5 Depending on the combination of the above, operations 3' to 6' may become unstable and the output time series may diverge.

[0061] The divergence of the output time series occurs due to the use of the transformation (Boxer-Thaler integrator) shown in the following equation (31). s -2 =(ΔT / 12 1 / 2 ) 2 (1+ 10z -1 +z -2 ) / (1- z -1 ) 2 (31) The conversion shown in equation (31) has high calculation accuracy, but the resulting digital filter may be unstable depending on the value of ΔT.

[0062] On the other hand, if the transformation (Tustin integrator) shown in the following equation (32) is used, the resulting digital filter is stable regardless of the value of ΔT. s -2 = s -1 s -1 = (ΔT / 2)(1+ z -1 ) / (1- z -1 ) (ΔT / 2)(1+ z -1 ) / (1- z -1 ) =(ΔT / 4 1 / 2 ) 2 (1+ 2z -1 +z -2 ) / (1- z -1 ) 2 (32) However, the conversion shown in equation (32) has low calculation accuracy.

[0063] Furthermore, the transformation (Madwed integrator) shown in the following equation (33) is also possible. s -2 =(ΔT / 6 1 / 2 ) 2 (1+ 4z -1 +z -2 ) / (1- z -1 ) 2 (33) The transformation shown in equation (33) has intermediate properties between the above two transformations in terms of stability and calculation accuracy.

[0064] These three transformations of equations (31) to (33) are s -2 =γ(ΔT) 2 (1+(1 / γ-2)z -1 +z -2 ) / (1- z -1 ) 2 (34) Therefore, in this embodiment, a digital filter is configured using the conversion of equation (34) so ​​that an optimal γ can be used. Note that in this embodiment, the value of γ is not limited to 1 / 12, 1 / 6, or 1 / 4. Here, γ is determined by the approximation method of the quadratic integral of the filter. This is the parameter that specifies the

[0065] By using the conversion of equation (34), in this embodiment, instead of the digital filter expressed by equations (13') and (14') shown in Patent Document 2, a digital filter expressed by the following equations (13) and (14) is obtained. -2 Different values ​​of γ for each (in equations (13) and (14), γ b1、 Gamma b2、 Gamma c1、 Gamma c2 ) is possible.

[0066] Third filter α 0 =1+ h b2 ω b ΔT +γ b2 ω b 2 (ΔT)2 (13) a 1 =-2+ω b 2 (ΔT) 2 (1-2c) b2 ) a 2 =1- h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 b 0 =1+ h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 b 1 =-2+ω b 2 (ΔT) 2 (1-2c) b1 ) b 2 =1- h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 oh b =2πf b

[0067] · 4th filter a 0 =1+ h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 (14) a 1 =-2+ω c 2 (ΔT) 2 (1-2c) c2 ) a 2 =1- h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 β 0 =γ c1 ω c 2 (ΔT) 2 β 1 =ω c 2 (ΔT) 2 (1-2γ c1 ) β 2 =γ c1 ω c 2 (ΔT) 2 ω c =2πf c

[0068] FIG. 9 is a block diagram of a time domain filter device according to this embodiment.

[0069] In summary, the time domain filter device for filtering the time series of ground acceleration in this embodiment includes a first filter 721 for filtering in the following equation (11), a second filter 722 for filtering in the following equation (12), a third filter 723 for filtering in the following equation (13), a fourth filter 724 for filtering in the following equation (14), an output time series acquisition unit 725 for applying these filters to an input time series x(k) to acquire an output time series y(k) in the following equation (15), a gain adjustment unit 726 for adjusting the gain in the following equation (16), and a calculation unit 727 for executing Calculations 11 to 17. And, at least one of Calculations 11 to 17 is executed in the following equations (11) to (16). First filter α 0 =8 / (ΔT) 2 +(4ω a1 +2ω a2 ) / (ΔT)+ω a1 ω a2 (11) α 1 =2ω a1 ω a2 -16 / (ΔT) 2 α 2=8 / (ΔT) 2 -(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 b 0 =4 / (ΔT) 2 +2h a2 / (ΔT) b 1 =-8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -2h a2 / (ΔT) oh a1 =2πf a1 , oh a2 =2πf a2 · 2nd filter a 0 =16 / (ΔT) 2 +17h a3 / (ΔT)+ ω a3 2 (12) a 1 =2h a3 2 -32 / (ΔT) 2 a 2 =16 / (ΔT) 2 -17h a3 / (ΔT)+ ω a3 2 b 0 =4 / (ΔT) 2 +8.5h a3 / (ΔT) + ω a3 2 b 1 =2h a3 2 -8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -8.5h a3 / (ΔT) + ω a3 2 oh a3 =2πf a3 · 3rd filter a 0 =1+ h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 (13) a 1 =-2+ω b 2 (ΔT) 2 (1-2c) b2 ) a 2 =1- h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 b 0 =1+ h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 b 1 =-2+ω b 2 (ΔT) 2 (1-2c) b1 ) b 2 =1- h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 oh b =2πf b · 4th filter a 0 =1+ h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 (14) a 1 =-2+ω c 2(ΔT) 2 (1-2γ c2 ) α 2 =1- h c ω c ΔT +γ c2 ω c 2 (ΔT) 2 β 0 =γ c1 ω c 2 (ΔT) 2 β 1 =ω c 2 (ΔT) 2 (1-2γ c1 ) β 2 =γ c1 ω c 2 (ΔT) 2 ω c =2πf c Output time series acquisition section y(k)=[-α 1 y(k-1) -α 2 y(k-2)+β 0 x(k) +β 1 x(k-1) +β 2 x(k-2)] / α 0 (15) however, x(k) is the input time series of the filter. y(k) is the output time series of the filter, k is the number of time steps. ΔT is the sampling interval. π is the ratio of the circumference of a circle to its circumference. Let us assume that. Gain adjustment section y(k)=g d x(k) (16) However, g d is the gain. ·Arithmetic section Operation 11:f a1 = f 0 , fa2 = f 1 Then, the calculation of equations (11) and (15) Operation 12:f a3 = f 1 Then, the calculation of equations (12) and (15) Operation 13:h b1 = h 2a , h b2 = h 2b , f b = f 2 , γ b1 = γ 2a , γ b2 = γ 2b Then, the calculation of equations (13) and (15) Operation 14:h c = h 3 , f c = f 3 , γ c1 = γ 3a , γ c2 = γ 3b Then, the calculation of equations (14) and (15) Operation 15:h c = h 4 , f c = f 4 , γ c1 = γ 4a , γ c2 = γ 4b Then, the calculation of equations (14) and (15) Operation 16:h c = h 5 , f c = f 5 , γ c1 = γ 5a , γ c2 = γ 5b Then, the calculation of equations (14) and (15) Operation 17:g d = g, and calculate equation (16). however, f 0 is the parameter that defines the cutoff frequency of the filter, f 1 is a parameter that specifies the position of the frequency domain where the signal is attenuated at 0.5 order. f 2 , h2a , h 2b , are parameters that define the characteristics of the correction filter, f 3 , h 3 are the parameters that define the cutoff frequency and damping of the filter, f 4 , h 4 are the parameters that define the cutoff frequency and damping of the filter, f 5 , h 5 are the parameters that define the cutoff frequency and damping of the filter, g is the gain adjustment parameter, Gamma 2a、 Gamma 2b、 Gamma 3a、 Gamma 3b, Gamma 4a、 Gamma 4b、 Gamma 5a、 Gamma 5b is a parameter that specifies the approximation method of the second-order integral of the filter, It is.

[0070] Next, we derive the condition for γ necessary to obtain a stable digital filter. For a digital filter to be stable, all poles must be inside the unit circle, and the following equations (21) to (24) are the conditions for a stable digital filter. (1 / 4-1 / (2πf 2 ΔT) 2 )≦γ 2b (twenty one) (1 / 4-1 / (2πf 3 ΔT) 2 )≦γ 3b (twenty two) (1 / 4-1 / (2πf 4 ΔT) 2 )≦γ 4b (twenty three) (1 / 4-1 / (2πf 5 ΔT) 2 )≦γ 5b (twenty four)

[0071] Therefore, f b、 f cGiven a given value, γ 2b , γ 3b , γ 4b , γ 5b can be selected according to ΔT. In this case, γ 2a , γ 3a , γ 4a , γ 5a The choice of is not affected by stability requirements.

[0072] Furthermore, the time domain filter device may satisfy at least one of the following expressions (25) to (28) in parameters defining an approximation method of the second-order integral of the filter. Gamma 2a = γ 2b (twenty five) Gamma 3a = γ 3b (26) Gamma 4a = γ 4b (27) Gamma 5a = γ 5b (28)

[0073] Regarding the actual selection of γ, 2 =0.5[Hz],f 3 =12.0[Hz],f 4 =20.0[Hz],f 5 If =30.0[Hz], the sampling frequency f s = 1 / ΔT, and then we can do as follows.

[0074] (1) When five types of γ are used, f s When ≥ 80[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 12 80[Hz] > f s When ≧ 70[Hz] Gamma 2b = 1 / 12, γ 3b= 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 8 70[Hz] > f s When ≧ 60[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 6 60[Hz] > f s When ≧ 50[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 10, γ 5b = 1 / 4 50[Hz] > f s When ≧ 40[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 6, γ 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 10, γ 4b = 1 / 4, γ 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] Gamma 2b = 1 / 6, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4

[0075] (2) When four types of γ are used, f s When ≥ 80[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 12 80[Hz] > f s When ≥ 70[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 8 70[Hz] > f s When ≧ 60[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 6 60[Hz] > f s When ≧ 50[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 8, γ 5b = 1 / 4 50[Hz] > f s When ≧ 40[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 6, γ 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 8, γ 4b = 1 / 4, γ 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] Gamma 2b = 1 / 6, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4

[0076] (3) When three types of γ are used, f s When ≥ 80[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 12 80[Hz] > f s When ≧ 60[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 6 60[Hz] > f s When ≧ 40[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 6, γ 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 6, γ 4b = 1 / 4, γ 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] Gamma 2b = 1 / 6, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4

[0077] (4) When two types of γ are used, f s When ≥ 80[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 12 80[Hz] > f s When ≧ 60[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 4 60[Hz] > f s When ≧ 40[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 4, γ 5b = 1 / 4 40[Hz] > f s When ≧ 5[Hz] Gamma 2b = 1 / 12, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] Gamma 2b = 1 / 4, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4

[0078] (5) When using one type of γ to improve the calculation accuracy, Gamma 2b = 1 / 12, γ 3b = 1 / 12, γ 4b = 1 / 12, γ 5b = 1 / 12 This γ is f 2 =0.5[Hz],f 3 =12.0[Hz],f 4 =20.0[Hz],f 5 = 30.0 [Hz], it is also applicable to cases other than this.

[0079] (6) When using one type of γ and ensuring stability, Gamma 2b = 1 / 4, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4 This γ is f 2 =0.5[Hz],f 3 =12.0[Hz],f 4 =20.0[Hz],f 5 = 30.0 [Hz], it is also applicable to cases other than this.

[0080] By using equations (13) and (14), a digital filter that is always stable can be obtained according to the sampling interval ΔT. Therefore, it is possible to perform time domain filter calculations even for a ground acceleration time series in which the sampling interval ΔT has different values ​​for each step.

[0081] Sampling frequency of ground acceleration time series f sWhen f is a low value such as 5[Hz], 5 = 30 [Hz], a filter operation with a high cutoff frequency is a wasteful process. In such a case, the corresponding filter process can be bypassed. Even in such a case, the filter operation of this embodiment can be operated.

[0082] Fig. 10 shows the earthquake waveform, the waveform obtained by interpolating the earthquake waveform using a conventional method, and the waveform obtained by using the seismic intensity estimation device of this embodiment. Fig. 10(a) shows the earthquake waveform of 50 Hz, Fig. 10(b) shows the waveform obtained by interpolating the 50 Hz earthquake waveform to 100 Hz and estimating the seismic intensity using the conventional method, and Fig. 10(c) shows the waveform calculated at 50 Hz using the seismic intensity estimation device of this embodiment.

[0083] Comparing Fig. 10(b) and Fig. 10(c), they show almost the same results. However, the conventional device shown in Fig. 10(b) performs interpolation, so it is inferior in real-timeness compared to the seismic intensity estimation device of this embodiment shown in Fig. 10(c).

[0084] FIG. 11 shows a case where the measured seismic intensity estimation device of this embodiment is constructed as a system. In the figure, 101 is a measured seismic intensity estimation system, 102 is a seismometer, 103 is a communication network, and 104 is an earthquake data processing device. In this system, seismic motion observed by multiple seismometers 102 is transferred to other devices or a data center through a communication network 103, and an approximate seismic intensity value is calculated by an earthquake data processing device 104 having a measured seismic intensity estimation unit 7 at a remote location. Here, there may be multiple earthquake data processing devices performing calculation processing. Even in such a system, the advantage of the present method that an approximate seismic intensity can be obtained as quickly as possible is utilized. Since data from many observation points must often be processed simultaneously in a data center, the algorithm of the present method that does not use FFT is also useful in reducing the calculation load. [Explanation of symbols]

[0085] 1...Measured seismic intensity estimation device, 2...Processing unit, 3...Control / calculation unit, 4...Measured seismic intensity calculation unit, 5...SI value calculation unit, 6...Response value calculation unit, 7...Measured seismic intensity estimation unit, 71...Ground acceleration time series acquisition unit, 72...Time domain filter device, 721...First filter, 722...Second filter, 723...Third filter, 724...Fourth filter, 725...First filter output time series acquisition unit, 726...Gain adjustment unit, 727...Calculation unit, 73...Calculation unit, 74...Amplitude time series conversion unit, 75...Seismic intensity conversion amplitude time series conversion unit, 76...Discretization unit, 77...Duration count unit, 8...Power supply unit, 9...Measuring device, 101...Measured seismic intensity estimation system, 102...Seismometer, 103...Communication network, 104...Earthquake data processing device

Claims

1. 1. A time domain filter apparatus for filtering a time series of ground acceleration, comprising: In the following formula (11), a first filter for performing filtering is In the following formula (12), a second filter for performing filtering is In the following formula (13), a third filter for performing filtering processing; In the following formula (14), a fourth filter for performing filtering processing; In the following formula (15), an output time series acquisition unit applies these filters to an input time series x(k) to acquire an output time series y(k); A gain adjustment unit that adjusts the gain in the following formula (16); A calculation unit that executes calculations 11 to 17; Equipped 2. A time domain filter device comprising: First filter a 0 =8 / (ΔT) 2 +(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 (11) a 1 =2h a1 oh a2 -16 / (ΔT) 2 a 2 =8 / (ΔT) 2 -(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 b 0 =4 / (ΔT) 2 +2h a2 / (ΔT) b 1 =-8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -2h a2 / (ΔT) oh a1 =2πf a1 , oh a2 =2πf a2 Second filter a 0 =16 / (ΔT) 2 +17h a3 / (ΔT)+ ω a3 2 (12) a 1 =2h a3 2 -32 / (ΔT) 2 a 2 =16 / (ΔT) 2 -17h a3 / (ΔT)+ ω a3 2 b 0 =4 / (ΔT) 2 +8.5h a3 / (ΔT) + ω a3 2 b 1 =2h a3 2 -8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -8.5h a3 / (ΔT) + ω a3 2 oh a3 =2πf a3 ・Third filter a 0 =1+ h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 (13) a 1 =-2+ω b 2 (ΔT) 2 (1-2c) b2 ) a 2 =1- h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 b 0 =1+ h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 b 1 =-2+ω b 2 (ΔT) 2 (1-2c) b1 ) b 2 =1- h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 oh b =2πf b ・Fourth filter a 0 =1+ h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 (14) a 1 =-2+ω c 2 (ΔT) 2 (1-2c) c2 ) a 2 =1- h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 b 0 =c c1 oh c 2 (ΔT) 2 b 1 =ω c 2 (ΔT) 2 (1-2c) c1 ) b 2 =c c1 oh c 2 (ΔT) 2 oh c =2πf c ・Output time series acquisition section y(k)=[-α 1 y(k-1) -a 2 y(k-2)+β 0 x(k) +β 1 x(k-1) +β 2 x(k-2)] / a 0 (15) however, x(k) is the input time series of the filter. y(k) is the output time series of the filter, k is the number of time steps. ΔT is the sampling interval. π is the ratio of the circumference of a circle to its circumference. Let us assume that. ・Gain adjustment section y(k)=g d x(k) (16) However, g d is the gain. ・Calculation section Operation 11: f a1 = f 0 , f a2 = f 1 Then, the calculation of equations (11) and (15) is performed. Operation 12: f a3 = f 1 Then, the calculation of equations (12) and (15) is performed. Operation 13: h b1 = h 2a , h b2 = h 2b , f b = f 2 , γ b1 = γ 2a , γ b2 = γ 2b Then, the calculation of equations (13) and (15) is performed. Operation 14: h c = h 3 , f c = f 3 , γ c1 = γ 3a , γ c2 = γ 3b Then, the calculation of equations (14) and (15) is performed. Operation 15: h c = h 4 , f c = f 4 , γ c1 = γ 4a , γ c2 = γ 4b Then, the calculation of equations (14) and (15) is performed. Operation 16: h c = h 5 , f c = f 5 , γ c1 = γ 5a , γ c2 = γ 5b Then, the calculation of equations (14) and (15) is performed. Operation 17: g d = g, and calculate equation (16). however, f 0 is the parameter that defines the cutoff frequency of the filter, f 1 is a parameter that specifies the position of the frequency domain where the signal is attenuated at 0.5 order. f 2 , h 2a , h 2b , are parameters that define the characteristics of the correction filter, f 3 , h 3 are the parameters that define the cutoff frequency and damping of the filter, f 4 , h 4 are the parameters that define the cutoff frequency and damping of the filter, f 5 , h 5 are the parameters that define the cutoff frequency and damping of the filter, g is the gain adjustment parameter, Gamma 2a、 Gamma 2b、 Gamma 3a、 Gamma 3b, Gamma 4a、 Gamma 4b、 Gamma 5a、 Gamma 5b is the approximate method for the quadratic integral of the filter. The parameters to be specified, and The following formulas (21) to (24) are satisfied. (1 / 4-1 / (2πf 2 ΔT) 2 )≦γ 2b (21) (1 / 4-1 / (2πf 3 ΔT) 2 )≦γ 3b (22) (1 / 4-1 / (2πf 4 ΔT) 2 )≦γ 4b (23) (1 / 4-1 / (2πf 5 ΔT) 2 )≦γ 5b (24)

2. In the parameters defining the approximation method of the second-order integral of the filter, at least one of the following expressions (25) to (28) is satisfied.

2. A time domain filter device according to claim 1. c 2a = c 2b (25) c 3a = c 3b (26) c 4a = c 4b (27) c 5a = c 5b (28)

3. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (1):

3. A time domain filter device according to claim 1 or 2. Condition (1) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≥ 70[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 8 70[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 6 60[Hz] > f s When ≧ 50[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 10, c 5b = 1 / 4 50[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 6, c 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] c 2b = 1 / 12, c 3b = 1 / 10, c 4b = 1 / 4, c 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 6, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

4. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (2):

3. A time domain filter device according to claim 1 or 2. Condition (2) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≧ 70[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 8 70[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 6 60[Hz] > f s When ≧ 50[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 8, c 5b = 1 / 4 50[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 6, c 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] c 2b = 1 / 12, c 3b = 1 / 8, c 4b = 1 / 4, c 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] γ 2b = 1 / 6, γ 3b = 1 / 4, γ 4b = 1 / 4, γ 5b = 1 / 4

5. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (3):

3. A time domain filter device according to claim 1 or 2. Condition (3) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 6 60[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 6, c 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] c 2b = 1 / 12, c 3b = 1 / 6, c 4b = 1 / 4, c 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 6, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

6. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (4):

3. A time domain filter device according to claim 1 or 2. Condition (4) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 4 60[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 4, c 5b = 1 / 4 40[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 4, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

7. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (5):

3. A time domain filter device according to claim 1 or 2. Condition (5) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12

8. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (6):

3. A time domain filter device according to claim 1 or 2. Condition (6) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], c 2b = 1 / 4, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

9. The sampling interval ΔT is set to a different value for each step. A time domain filter device according to any one of claims 1 to 8.

10. a ground acceleration time series acquisition unit for acquiring a time series of ground acceleration; The time domain filter device according to claim 1 , which filters the time series of the ground acceleration acquired by the ground acceleration time series acquisition unit; a calculation unit for calculating an approximate value of a measured seismic intensity from the time series of the ground acceleration filtered by the time domain filter device; Equipped with The calculation unit is an amplitude time series converter for converting the ground acceleration time series filtered by the time domain filter device into an amplitude time series; a seismic intensity converted amplitude time series conversion unit that converts the amplitude time series converted by the amplitude time series conversion unit into a seismic intensity converted amplitude time series; a discretization unit that discretizes the seismic intensity converted amplitude time series obtained by the seismic intensity converted amplitude time series conversion unit; A duration counting unit that executes a duration count for each discretized seismic intensity conversion amplitude value obtained by the discretization unit; have A seismic intensity estimation device characterized by the above.

11. A seismic intensity approximation device according to claim 10, Multiple seismometers that measure seismic motion; A communication network; an earthquake data processing device that performs approximate seismic intensity calculation using the seismic intensity approximation device at a remote location via the communication network based on the seismic motion observed by the plurality of seismometers; Equipped A seismic intensity estimation system characterized by the above.

12. 1. A time domain filtering method for filtering a time series of ground acceleration, comprising the steps of: In the following formulas (11) to (16), at least one of operations 11 to 17 is executed.

13. A time domain filtering method comprising: First filter a 0 =8 / (ΔT) 2 +(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 (11) a 1 =2h a1 oh a2 -16 / (ΔT) 2 a 2 =8 / (ΔT) 2 -(4h a1 +2h a2 ) / (ΔT)+ω a1 oh a2 b 0 =4 / (ΔT) 2 +2h a2 / (ΔT) b 1 =-8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -2h a2 / (ΔT) oh a1 =2πf a1 , oh a2 =2πf a2 Second filter a 0 =16 / (ΔT) 2 +17h a3 / (ΔT)+ ω a3 2 (12) a 1 =2h a3 2 -32 / (ΔT) 2 a 2 =16 / (ΔT) 2 -17h a3 / (ΔT)+ ω a3 2 b 0 =4 / (ΔT) 2 +8.5h a3 / (ΔT) + ω a3 2 b 1 =2h a3 2 -8 / (ΔT) 2 b 2 =4 / (ΔT) 2 -8.5h a3 / (ΔT) + ω a3 2 oh a3 =2πf a3 ・Third filter a 0 =1+ h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 (13) a 1 =-2+ω b 2 (ΔT) 2 (1-2c) b2 ) a 2 =1- h b2 oh b ΔT +γ b2 oh b 2 (ΔT) 2 b 0 =1+ h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 b 1 =-2+ω b 2 (ΔT) 2 (1-2c) b1 ) b 2 =1- h b1 oh b ΔT +γ b1 oh b 2 (ΔT) 2 oh b =2πf b ・Fourth filter a 0 =1+ h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 (14) a 1 =-2+ω c 2 (ΔT) 2 (1-2c) c2 ) a 2 =1- h c oh c ΔT +γ c2 oh c 2 (ΔT) 2 b 0 =c c1 oh c 2 (ΔT) 2 b 1 =ω c 2 (ΔT) 2 (1-2c) c1 ) b 2 =c c1 oh c 2 (ΔT) 2 oh c =2πf c ・The output time series y(k) is the result of applying these filters to the input time series x(k). y(k)=[-α 1 y(k-1) -a 2 y(k-2)+β 0 x(k) +β 1 x(k-1) +β 2 x(k-2)] / a 0 (15) however, x(k) is the input time series of the filter. y(k) is the output time series of the filter, k is the number of time steps. ΔT is the sampling interval. π is the ratio of the circumference of a circle to its circumference. Let us assume that. ・Gain adjustment y(k)=g d x(k) (16) However, g d is the gain. Operation 11: f a1 = f 0 , f a2 = f 1 Then, the calculation of equations (11) and (15) is performed. Operation 12: f a3 = f 1 Then, the calculation of equations (12) and (15) is performed. Operation 13: h b1 = h 2a , h b2 = h 2b , f b = f 2 , γ b1 = γ 2a , γ b2 = γ 2b Then, the calculation of equations (13) and (15) is performed. Operation 14: h c = h 3 , f c = f 3 , γ c1 = γ 3a , γ c2 = γ 3b Then, the calculation of equations (14) and (15) is performed. Operation 15: h c = h 4 , f c = f 4 , γ c1 = γ 4a , γ c2 = γ 4b Then, the calculation of equations (14) and (15) is performed. Operation 16: h c = h 5 , f c = f 5 , γ c1 = γ 5a , γ c2 = γ 5b Then, the calculation of equations (14) and (15) is performed. Operation 17: g d = g, and calculate equation (16). however, f 0 is the parameter that defines the cutoff frequency of the filter, f 1 is a parameter that specifies the position of the frequency domain where the signal is attenuated at 0.5 order. f 2 , h 2a , h 2b , are parameters that define the characteristics of the correction filter, f 3 , h 3 are the parameters that define the cutoff frequency and damping of the filter, f 4 , h 4 are the parameters that define the cutoff frequency and damping of the filter, f 5 , h 5 are the parameters that define the cutoff frequency and damping of the filter, g is the gain adjustment parameter, γ 2a、 γ 2b、 γ 3a、 γ 3b, γ 4a、 γ 4b、 γ 5a、 γ 5b is a parameter that defines an approximation method for the second integral of the filter. and Gamma 2b , γ 3b , γ 4b , γ 5b satisfies the following expressions (21) to (24). (1 / 4-1 / (2πf 2 ΔT) 2 )≦γ 2b (21) (1 / 4-1 / (2πf 3 ΔT) 2 )≦γ 3b (22) (1 / 4-1 / (2πf 4 ΔT) 2 )≦γ 4b (23) (1 / 4-1 / (2πf 5 ΔT) 2 )≦γ 5b (24)

13. In the parameters that define the approximation method of the second-order integral of the filter, at least one of the following expressions (25) to (28) is satisfied.

13. The time domain filtering method of claim 12. c 2a = c 2b (25) c 3a = c 3b (26) c 4a = c 4b (27) c 5a = c 5b (28)

14. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (1): A time domain filtering method according to claim 12 or 13. Condition (1) Let the sampling frequency be f s = 1 / ΔT, and f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≧ 70[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 8 70[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 6 60[Hz] > f s When ≧ 50[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 10, c 5b = 1 / 4 50[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 6, c 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] c 2b = 1 / 12, c 3b = 1 / 10, c 4b = 1 / 4, c 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 6, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

15. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (2): A time domain filtering method according to claim 12 or 13. Condition (2) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≥ 70[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 8 70[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 6 60[Hz] > f s When ≧ 50[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 8, c 5b = 1 / 4 50[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 6, c 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] c 2b = 1 / 12, c 3b = 1 / 8, c 4b = 1 / 4, c 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 6, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

16. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (3): A time domain filtering method according to claim 12 or 13. Condition (3) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 6 60[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 6, c 5b = 1 / 4 40[Hz] > f s When ≧ 30[Hz] c 2b = 1 / 12, c 3b = 1 / 6, c 4b = 1 / 4, c 5b = 1 / 4 30[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 6, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

17. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (4): A time domain filtering method according to claim 12 or 13. Condition (4) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], f s ≧ 80 When [Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12 80[Hz] > f s When ≧ 60[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 4 60[Hz] > f s When ≧ 40[Hz] c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 4, c 5b = 1 / 4 40[Hz] > f s When ≧ 5[Hz] c 2b = 1 / 12, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4 5[Hz] > f s When ≧ 1[Hz] c 2b = 1 / 4, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

18. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (5): A time domain filtering method according to claim 12 or 13. Condition (5) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], c 2b = 1 / 12, c 3b = 1 / 12, c 4b = 1 / 12, c 5b = 1 / 12

19. The parameters that define the approximation method of the quadratic integral of the filter, the parameters that define the characteristics of the filter, and the parameters that define the cutoff frequency of the filter satisfy the following condition (6): A time domain filtering method according to claim 12 or 13. Condition (6) Let the sampling frequency be f s = 1 / ΔT, f 2 = 0.5 [Hz], f 3 = 12.0 [Hz], f 4 = 20.0 [Hz], f 5 If = 30.0 [Hz], c 2b = 1 / 4, c 3b = 1 / 4, c 4b = 1 / 4, c 5b = 1 / 4

20. The sampling interval ΔT is set to a different value for each step. A time domain filtering method according to any one of claims 12 to 19.

Images