The present invention relates to a method and device for designing a digital filter, a program for designing a digital filter, a digital filter, a method and device for generating a numerical sequence of a desired frequency characteristic, and a program for generating a numerical sequence of a desired frequency characteristic, and in particular, to an FIR filter of a type that comprises tapped delay line made up of a plurality of delayers and which multiplies signals from respective taps severalfold and adds up multiplication results for output, and a method of designing this filter, as well as a method for generating a numerical sequence which is used to design the filter and which indicates an input frequency characteristic.
Various kinds of electronical devices provided in a variety of technical fields normally implement digital signal processing of some sort in their inside. The most important basic operations of digital signal processing include filtering processing of taking only signals within a required certain frequency band out of input signals in which respective kinds of signals and noises are mixed. Therefore, digital filters are frequently used in electronics devices of implementing digital signal processing.
IIR (Infinite Impulse Response) filters and FIR (Finite Impulse Response) filters are mostly used as digital filters. Among them, the FIR filters are advantageous as follows. Firstly the circuit is always stable since the pole of transfer function of an FIR filter is located only in the origin of the z plane. Secondly, if the filter coefficients are of a symmetrical type, it is possible to realize a completely accurate linear-phase characteristic.
In this FIR filter, the impulse response expressed in finite time length will straight be the filter coefficients. Accordingly, designing an FIR filter means to determine the filter coefficients so as to obtain a desired frequency characteristic. Several methods for calculating filter coefficients have been proposed.
For example, one of these methods determines filter coefficients by a convolution or the like using a Chebyshev approximation on the basis of the ratio of a sampling frequency to a cutoff frequency for a desired frequency characteristic. Another method determines filter characteristics by inputting a waveform for a desired frequency characteristic in the form of a numerical sequence or a function, subjecting the input numerical sequence or function to an inverse Fourier transformation (inverse FFT), and extracting real items from the result (see, for example, Patent Documents 1 and 2).
Patent Document 1: Japanese Patent Laid-Open No. 63-234617
Patent Document 2: Japanese Patent Laid-Open No. 2003-168958
However, the above conventional techniques determine an enormous number of filter coefficients of very complicated, random values. Thus, using all the filter coefficients obtained sharply increases the number of taps in a filter circuit and requires a large number of multipliers to subject the complicated, random filter coefficient values to multiplication. In other words, the above technique requires a large-scale circuit configuration. This is not practical. Thus, the filter coefficients need to be reduced to an appropriate number in a practical sense by means of windowing using a window function.
However, the windowing for a reduction in filter coefficients often discretize the filter coefficients, which significantly affect the frequency characteristic. This prevents an appropriate target frequency characteristic from being obtained. Further, the method of determining a filter coefficient subjecting the numerical sequence of an input waveform to an inverse FFT determines a frequency characteristic depending on the numerical sequence or function that expresses the input waveform. However, the determination of the numerical sequence or function is itself difficult. Thus, it is disadvantageously very difficult to offer a desired frequency characteristic regardless of whichever conventional filter designing method is used.
To obtain a desired frequency characteristic by the conventional filter designing method based on windowing, a trial and error process needs to be executed by subjecting temporarily determined filter coefficients to an FFT with the resulting frequency characteristic checked. In particular, with a technique for subjecting an input waveform to an inverse FFT, the numerical sequence or function of the input waveform must itself be determined through a trial and error process. Thus, disadvantageously, the conventional techniques require a skilled technician to put much time and effort in designing an FIR filter. This prevents an FIR filter with a desired characteristic from being easily designed.
The present invention has been made in order to solve the above problems. An object of the present invention is to allow an FIR digital filter to be easily designed with almost no trial and error processes, the FIR digital filter involving a reduced number of filter coefficients and enabling a desired frequency characteristic to be accurately provided using a small-scale circuit.
To accomplish the object, the present invention inputs a standard function to a circuit and calculates an interpolation function of a finite length from the standard function to determine an input frequency characteristic based on specifications. The present invention then subjects a numerical sequence indicative of the input frequency characteristic to an inverse Fourier transformation to obtain filter coefficients. The present invention then executes a rounding process based on the coefficient values to obtain a reduced number of filter coefficients depending on the number of process bits.
The present invention configured as described above allows an FIR digital filter having a desired frequency characteristic to be easily designed without any expertise; examples of such an FIR digital filter include a low pass filter, a high pass filter, a band pass filter, and a band elimination filter. The present invention doesn't require the windowing for reducing the number of filter coefficients. The present invention uses a numerical rounding operation to enable the number of filter coefficients (the number of taps for the digital filter) to be reduced without lowering the accuracy of the frequency characteristic. That is, the present invention enables an FIR filter with an appropriate frequency characteristic to be easily designed; the FIR filter requires a reduced number of taps and offers a pass band characteristic that allows ripple to be minimized as well as a uniform attenuation characteristic.
An embodiment of the present invention will be described with reference to the drawings.
As shown in
Then, a numerical sequence determined by the thus input interpolation function is subjected to an inverse FFT, with the resulting real items extracted (step S2). As is well known, executing an FFT process on a numerical sequence results in a frequency characteristic corresponding to the numerical sequence. Accordingly, an inverse FFT is executed on a numerical sequence indicating a frequency characteristic which is input through the interpolation function and the resulting real items are extracted, the numerical sequences required to provide the input frequency characteristic are obtained. These numerical sequences correspond to filter coefficients to be determined.
However, the numerical sequences themselves, determined, by an inverse FFT, from the interpolation function calculated in step S1, are not always configured so that their values are arranged so as to be directly used as filter coefficients. Specifically, for any types of digital filters, the numerical sequences of filter coefficients are symmetric so that their central value is largest and so that the other values decrease consistently with increasing distance from the central value, with amplitude repeated. In contrast, numerical sequences determined from the interpolation function by an inverse FFT have the smallest value in the center and the largest value at the ends of the distribution.
Thus, as shown in
The numerical sequences thus obtained can be determined to be target filter coefficients as they are. However, the present embodiment further performs a rounding operation described below to reduce the filter coefficients to a required number, and simplifying their values (step S4).
For example, if the numerical sequences resulting from the appropriate rearrangement in step S3 are y-bit data, the y-bit data is multiplied by a factor of 2X, with the resulting decimal fractions rounded to integers. Thus, x-bit (x<y) integral data is obtained and utilized as filter coefficients. Alternatively, a rounding process may be executed on y-bit data to obtain x-bit (x<y) data, which may then be multiplied by a factor of 2X to obtain integers.
Performing such a rounding operation to obtain integers enables the digital filter to be configured as shown in
The present embodiment determines the numerical sequences determined by such a rounding operation to be target filter coefficients. The above steps S3 and S4 need not necessarily be executed in this order but may be reversed.
Now, a method of calculating the input frequency characteristic in step 1 will be described in detail with reference to a specific example. Here, two methods of determining the input frequency characteristic are shown.
<First Generation Method>
The numerical sequence proposed by the present inventor and described in Japanese Patent Application No. 2003-56265 is an example of a coefficient string having a finite-base impulse response such that the impulse response has a finite value other than “0” in a local area and a value of “0” in the other areas. For example, to design a low pass filter according to specifications shown in
XF1=8/16+9/16*cos(2πt)−1/16*cos(6πt) (Equation 1)
Here, the function shown in (Equation 1) is obtained by standardization using a maximum amplitude value of 1 and a maximum frequency value of 1. The coefficients {8/16, 9/16, 0, −1/16} (“0” is the coefficient of the item cos (4πt)) of the items of (Equation 1) constitute a numerical sequence corresponding to one of the halves into which the filter coefficients {−1, 0, 9, 16, 9, 0, −1}/16 of the low pass filter described in Japanese Patent Application No. 2003-56265 are divided at the center. As described in Japanese Patent Application No. 2003-56265 in detail, the impulse response from a low pass filter having the numerical sequence {−1, 0, 9, 16, 9, 0, −1}/16 as filter coefficients is of a finite-base and passes through all sample points required to provide a smooth waveform. The numerical sequence {8, 9, 0, −1}/16 also has a finite-base impulse response.
To design a low pass filter according to such specifications as shown in
After inputting the numerical sequence of the standard function XF1, an interpolation function is determined on the basis of the standard function XF1 (step S12). The interpolation function to be determined interpolates the range between amplitude values “1” and “0” of the frequency amplitude characteristic. To determine the interpolation function, first, the ratio of the transition area to entire area of the frequency characteristic determined by the standard function XF1 (hereinafter referred to as a standard transition area ratio Rts) is determined. Here, the transition area refers to the area of the inclined part between a pass band and a stop band. With the first generation method, the transition area is considered to be the area between two representative points in the inclined part (for example, the range of the amplitude from −0.3 dB to −45 dB).
If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB in the standard function XF1 is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the value for the standardization clock Td corresponding to these amplitude values in the first half of the frequency characteristic shown in
Then, an interpolation function length Li is determined from the standard transition area ratio Rts. The term “interpolation function length Li” refers to the length (standardization clock count) of effective area of the interpolation function to be determined. The interpolation function length Li is determined from the clock width Ld of the transition area of the FIR filter to be designed and the standard transition area ratio Rts. If a low pass filter according to specifications shown in
To reduce the number of the taps in a low pass filter, it is desirable that the interpolation function length Li be an even integer larger than the calculated value. Thus, in this case, the interpolation function length Li is set to 66. An interpolation function I (LPF1) with an interpolation function length Li of 66 clocks is determined as shown in the following partition equations (Equation 2-1) and (Equation 2-2).
I (LPF1)=8/16+9/16*cos(2πt/66)−1/16*cos(6πt/66) (0/1024≦t≦65/1024) (Equation 2-1)
I (LPF1)=0 (65/1024<t≦1023/1024) (Equation 2-2)
Specifically, the interpolation function I (LPF1) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 2-1 and 2-2) from 0/1024 to 1023/1024.
Once the interpolation function I (LPF1) is thus determined, its frequency characteristic is shifted in the direction of the frequency axis (clock direction) so that the shifted interpolation function I (LPF,) joins the amplitude values “1” and “0” together (step S13). Specifically, 66 numerical sequences corresponding to the positions of the standardization clock t=0/1024 to 65/1024 determined by (Equation 2-1) are shifted to the positions of the standardization clock t=i/1024 to (i+65)/1024 (i is an integer). And all the numerical sequences at the positions of the standardization clock t=0/1024 to (i−1)/1024 are changed to “1”, And all the numerical sequences at the positions of the standardization clock t=(i+66)/1024 to 1023/1024 are changed to “0”.
Then, the frequency characteristic shown in
The interpolation function shift amount i may be set at such a value as locates the amplitude value “0.5” of the interpolation function at positions on the frequency axis corresponding to its ⅛, 2/8, and ⅜. This simplifies the filter coefficients obtained by executing an inverse FFT on the numerical sequences for the input frequency characteristic in step S2 in
<Second Generation Method>
Now, description will be given of a second method for determining the input frequency characteristic.
In
After inputting the numerical sequence of the standard function XF1, an interpolation function is determined on the basis of the standard function XF1 (step S22). The interpolation function to be determined is the one which interpolates the range between amplitude values “1” and “0” of the frequency amplitude characteristic, and which has been subjected to frequency shifting on the basis of the design specifications for the requested digital filter.
To determine the interpolation function which has been subjected to frequency shifting, first, the ratio of the transition area (based on the design specifications for the digital filter) of an interpolation function to be generated from the standard function XF1 to the transition area of the standard function XF1 (this ratio is hereinafter referred to as a requested transition area ratio Rtr). Unlike the transition area in the first generation method, the transition area in the second generation method is that area which the amplitude takes a value other than “1” and “0” in the frequency amplitude characteristic whose amplitude values are standardized between “1” and “0”. To determine the requested transition area ratio Rtr, information on two representative points (for example, the points at the amplitudes of −0.3 dB and −45 dB) in the transition area is used.
If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB in the standard function XF1 is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the value for the standardization clock Td corresponding to these amplitude values in the first half of the frequency characteristic shown in
Then, a calculation is made of the clock count Lhs from the standardization clock t=0 for the requested digital filter to the start point t=k1 of the transition area (see
Moreover, a calculation is made of the clock count Lhe from the standardization clock t=0 for the requested digital filter to the end point t=k4 of the transition area. The clock count from the start point k1 to end point k4 of the transition area is defined as Tk1−k4. Then, the clock count Lhe from the standardization clock t=0 to the transition area end point t=k4 is determined by Lhe=Lhs+Tk1−k4. Here, the clock count Tk1−k4 from the start point k1 to end point k4 of the transition area is determined by 0.5*Rtr using the clock count (=0.5) from the transition area start point (at t=0) to end point (at t=511/1024) in the standard function XF1, as well as the requested transition area ratio Rtr. As described above, Rtr=0.126963, so that Tk1−k4=0.5*0.126963=0.063482. Therefore, the clock count Lhe from the standardization clock t=0 to the transition area end point t=k4 is determined to be Lhe=0.092553+0.063482=0.156035.
As a result, the interpolation function I (LPF1) is determined as shown in the following partition equations (Equations 3-1, 3-2, and 3-3).
I (LPF1)=1 (0/1024≦t<Lhs) (Equation 3-1)
I (LPF1)=8/16+9/16*cos((2π(t−Lhs)/Rtr))−1/16*cos((6π(t−Lhs)/Rtr)) (Lhs<t≦Lhe) (Equation 3-2)
I (LPF1)=0 (Lhe<t≦1023/1024) (Equation 3-3)
Specifically, the above interpolation function I (LPF1) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 3) from 0/1024 to 1023/1024. A graph showing these 1024 numerical sequences is almost similar to that in
Then, the frequency characteristic shown in
As shown in
This enables a rounding process to sharply reduce the unwanted filter coefficients. For example, by dropping lower several bits of the filter coefficient to reduce the bits, it is possible to round all the filter coefficients with values smaller than the maximum value expressed only by the lower several bits, to “0” for discarding.
Thus, the present embodiment can reduce the filter coefficients by performing a rounding operation utilizing coefficient values. Consequently, windowing, as utilized in the prior art, is not necessarily required. As described above, the standard function input in step S1 has a finite-base impulse response. Thus, the number of filter coefficients designed on the basis of this standard function is originally smaller than that in the prior art and may be used without rounding. However, to further reduce the taps, a rounding process is preferably executed to reduce the bits.
The present embodiment is markedly different from the conventional filter designing method in the above ability to perform a rounding operation utilizing coefficient values. That is, the conventional filter designing method does not provide a sufficiently sharp distribution for the filter coefficients to be determined. Consequently, a rounding process with filter coefficient values often discards the major filter coefficients, which determine the frequency characteristic. Further, it is also difficult to obtain a frequency characteristic exhibiting a very large out-of-band attenuation amount. This prevents the required out-of-band attenuation amount from being obtained when the bits of the filter coefficient are reduced. Thus, the conventional technique cannot execute a rounding process for reducing the bits and is thus forced to reduce the filter coefficients by means of windowing. This results in a discretization error in the frequency characteristic, making it very difficult to obtain the desired frequency characteristic.
In contrast, the present embodiment can design an FIR filter without windowing, preventing a possible discretization error in the frequency characteristic. This enables the cutoff characteristic to be very significantly improved, providing an excellent filter characteristic with a rectilinear phase characteristic. In other words, an appropriate frequency characteristic can be offered which exhibits a pass band characteristic with reduced ripple and a uniform attenuation characteristic.
As shown in
In the description of this example, the standard function XF1, such as the one shown in (Equation 1), is used to design a low pass filter. However, (Equation 1) is only illustrative. For example, a standard function XF2 or XF3 expressed by (Equation 4) or (Equation 5), respectively, may be used.
XF2=1/2+cos(2πt) (Equation 4)
XF3=cos(πt)+1/8*cos(3πt)−1/8*cos(5πt) (Equation 5)
Here, the coefficients {1/2, 1} of the items of (Equation 4) constitute a numerical sequence corresponding to one of the halves into which the numerical sequence {1, 2, 1}/2 is divided at the center. Further, the coefficients {1, 1/8, −1/8} of the items of (Equation 5) constitute a numerical sequence corresponding to one of the halves into which the numerical sequence {−1, 1, 8, 8, 1, −1}/8 for a basic low pass filter L4a3 described in Japanese Patent Application No. 2003-56265 is divided at the center.
Japanese Patent Application No. 2003-56265 shows several numerical sequences for a low pass filter which are different from the numerical sequences corresponding to the coefficients of the items of (Equation 1), (Equation 4), and (Equation 5). The functions corresponding to these numerical sequences may each be used as the standard function according to the present embodiment.
As shown in
Thus, even when filter coefficients with values smaller than a predetermined threshold are discarded by a rounding process, most of the major filter coefficients remain, which determine the frequency characteristic. Consequently, the frequency characteristic is not substantially affected. This enables a rounding process to sharply reduce the unwanted filter coefficients. For example, by dropping lower several bits of the filter coefficient to reduce the bits, it is possible to round all the filter coefficients with values smaller than the maximum value expressed only by the lower several bits, to “0” for discarding.
Now, description will be given of an example in which a high pass filter according to such specifications as shown in
<First Generation Method>
To generate an input frequency characteristic according to the above first generation method, first, the standard function XF4, such as the one shown in (Equation 6), is input to the circuit. The input standard function XF4 is of a finite-base such that its impulse response has a finite value other than “0” in a local area and a value of “0” in all the other areas.
XF4=8/16−9/16*cos(2πt)+1/16*cos(6πt) (Equation 6)
Here, the function shown in (Equation 6) is obtained by standardization using a maximum amplitude value of 1 and a maximum frequency value of 1. The coefficients {8/16, −9/16, 0, 1/16} (“0” is the coefficient of the item cos (4πt)) of the items of (Equation 6) constitute a numerical sequence corresponding to one of the halves into which the filter coefficients {1, 0, −9, 16, −9, 0, 1}/16 of the high pass filter described in Japanese Patent Application No. 2003-56265 are divided at the center. As described in Japanese Patent Application No. 2003-56265 in detail, the impulse response from a high pass filter having the numerical sequence {1, 0, −9, 16, −9, 0, 1}/16 as filter coefficients is of a finite-base and passes through all sample points required to provide a smooth waveform. The numerical sequence {8, −9, 0, 1}/16 also has a finite-base impulse response.
To design a high pass filter according to such specifications as shown in
After inputting the numerical sequence of the standard function XF4, an interpolation function is determined on the basis of the standard function XF4. To determine the interpolation function, first, the standard transition area ratio Rts of the frequency characteristic determined by the standard function XF4 is determined.
If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB in the standard function XF4 is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the value for the standardization clock Tu corresponding to these amplitude values in the first half of the frequency characteristic shown in
Then, an interpolation function length Li is determined from the standard transition area ratio Rts. If a high pass filter according to the specifications shown in
To reduce the number of the taps in a high pass filter, it is desirable that the interpolation function length Li be an even integer larger than the calculated value. Thus, in this case, the interpolation function length Li is set to 66. An interpolation function I (HPF) with an interpolation function length Li of 66 clocks is determined as shown in the following partition equations (Equation 7-1) and (Equation 7-2).
I (HPF)=8/16−9/16*cos(2πt/66)+1/16*cos(6πt/66) (0/1024≦t≦65/1024) (Equation 7-1)
I (HPF)=1 (65/1024<t≦1023/1024) (Equation 7-2)
Specifically, the interpolation function I (HPF) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 7-1 and 7-2) from 0/1024 to 1023/1024.
Once the interpolation function I (HPF) is thus determined, its frequency characteristic is shifted in the direction of the frequency axis (clock direction) so that the shifted interpolation function I (HPF) joins the amplitude values “1” and “0” together. Specifically, 66 numerical sequences corresponding to the positions of the standardization clock t=0/1024 to 65/1024 determined by (Equation 7-1) are shifted to the positions of the standardization clock t=i/1024 to (i+65)/1024 (i is an integer). And all the numerical sequences at the positions of the standardization clock t=0/1024 to (i−1)/1024 are changed to “0”, and all the numerical sequences at the positions of the standardization clock t=(i+66)/1024 to 1023/1024 are changed to “1”.
Then, the frequency characteristic shown in
The interpolation function shift amount i in this case may be also set at such a value as locates the amplitude value “0.5” of the interpolation function at positions on the frequency axis corresponding to its ⅛, 2/8, and ⅜. This simplifies the filter coefficients obtained by subjecting the numerical sequences for the input frequency characteristic to an inverse FFT in step S2 in
<Second Generation Method>
To generate an input frequency characteristic according to the second generation method, first, such a standard function XF4 as shown in (Equation 6) is input to the circuit. After inputting the numerical sequence of the standard function XF4, an interpolation function containing frequency shifts is determined on the basis of the standard function XF4.
To determine an interpolation function subjected to frequency shifting, first, the requested transition area ratio Rtr of the interpolation function corresponding to the standard function XF4 is determined. To determine the requested transition area ratio Rtr, information on two representative points (for example, the points at the amplitudes of −0.3 dB and −45 dB) in the transition area is used.
If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB in the standard function XF4 is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the value for the standardization clock Tu corresponding to these amplitude values in the first half of the frequency characteristic shown in
Then, a calculation is made of the clock count Lhs from the standardization clock t=0 for the requested digital filter to the start point t=k1 of the transition area (see
Moreover, a calculation is made of the clock count Lhe from the standardization clock t=0 for the requested digital filter to the end point t=k4 of the transition area. The clock count from the start point k1 to end point k4 of the transition area is defined as Tk1−k4. Then, the clock count Lhe from the standardization clock t=0 to the transition area end point t=k4 is determined by Lhe=Lhs+Tk1−k4. Here, the clock count Tk1−k4 from the start point k1 of the transition area to the end point k4 is determined by 0.5*Rtr using the clock count (=0.5) from the transition area start point (at t=0) to end point (at t=511/1024) in the standard function XF4, as well as the requested transition area ratio Rtr. As described above, Rtr=0.126963, so that Tk1−k4=0.5*0.126963=0.063482. Therefore, the clock count Lhe from the standardization clock t=0 to the transition area end point t=k4 is determined to be Lhe=0.097715+0.063482=0.161197.
As a result, the interpolation function I (HPF) is determined as shown in the following partition equations (Equations 8-1, 8-2, and 8-3).
Specifically, the above interpolation function I (HPF) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 8) from 0/1024 to 1023/1024. A graph showing these 1024 numerical sequences is almost similar to that in
Then, the frequency characteristic shown in
As shown in
Thus, even when filter coefficients with values smaller than a predetermined threshold are discarded by a rounding process, most of the major filter coefficients remain, which determine the frequency characteristic. Consequently, the frequency characteristic is not substantially affected. Further, the frequency characteristic obtained by the filter designing method according to the present embodiment exhibits a very significant attenuation. Consequently, the desired attenuation amount can be obtained even with a slight decrease in the number of bits.
This enables a rounding process to sharply reduce the unwanted filter coefficients. For example, by dropping lower several bits of the filter coefficient to reduce the bits, it is possible to round all the filter coefficients with values smaller than the maximum value expressed only by the lower several bits, to “0” for discarding.
Thus, the present embodiment can reduce the filter coefficients by performing a rounding operation utilizing coefficient values. Consequently, windowing, as utilized in the prior art, is not necessarily required. As described above, the standard function input in step S1 has a finite-base impulse response. Thus, the number of filter coefficients designed on the basis of this standard function is originally smaller than that in the prior art. Accordingly, the standard function can be used as it is without the need for a rounding process. However, to further simplify the circuit, a rounding process is preferably executed to reduce the bits.
As shown in
In the description of this example, such a standard function XF4 as shown in (Equation 6) is used to design a high pass filter. However, the clock may be shifted by 0.5 using the standard function XF1 for a low pass filter. Further, (Equation 6) is only illustrative. For example, a standard function XF5 or XF6 expressed by (Equation 9) or (Equation 10), respectively, may be used.
XF5=−1/2+sin(2πt) (Equation 9)
XF6=cos(πt)−1/8*cos(3πt)−1/8*cos(5πt) (Equation 10)
Here, the coefficients {−1/2, 1} of the items of (Equation 9) constitute a numerical sequence corresponding to one of the halves into which the numerical sequence {−1, 2, −1}/2 is divided at the center. Further, the coefficients {1, −1/8, −1/8} of the items of (Equation 10) constitute a numerical sequence corresponding to one of the halves into which the numerical sequence {1, 1, −8, 8, −1, −1}/8 for a basic high pass filter H4a3 described in Japanese Patent Application No. 2003-56265 is divided at the center.
Japanese Patent Application No. 2003-56265 shows several numerical sequences for a high pass filter which are different from the numerical sequences corresponding to the coefficients of the items of (Equation 6), (Equation 9), and (Equation 10). The functions corresponding to these numerical sequences may each be used as the standard function according to the present embodiment.
Although not particularly shown in the drawings, not only the standard function in (Equation 6) but also such a standard function as shown in (Equation 9) or (Equation 10) allows the filter coefficients determined by the filter designing method according to the present invention to exhibit the largest value in the center of the distribution (position of the clock t=512/1024). The filter coefficients have a very sharp distribution such that their values are larger in a local area close to the center and smaller in the other areas and such that there are very large differences between the filter coefficient values close to the center and the peripheral filter coefficient values.
Thus, even when filter coefficients with values smaller than a predetermined threshold are discarded by a rounding process, most of the major filter coefficients remain, which determine the frequency characteristic. Consequently, the frequency characteristic is not substantially affected. This enables a rounding process to sharply reduce the unwanted filter coefficients. For example, by dropping lower several bits of the filter coefficient to reduce the bits, it is possible to round all the filter coefficients with values smaller than the maximum value expressed only by the lower several bits, to “0” for discarding.
Now, description will be given of an example in which a band pass filter according to such specifications as shown in
<First Generation Method>
To generate an input frequency characteristic according to a first generation method, first, for example, the standard function XF1 for a low pass filter, shown in (Equation 1), and the standard function XF4 for a high pass filter, shown in (Equation 6) (see
After inputting the standard functions XF1 and XF4, an interpolation function is determined on the basis of the standard functions XF1 and XF4. First, standard transition area ratios RtsL and RtsH of an identified frequency characteristic are determined from the standard functions XF1 and XF4.
If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB is 0.966051. The amplitude value of −45 dB is 0.005623. For example, a calculation is made of the value for the standardization clock Td corresponding to these amplitude values in the first half of the frequency characteristic shown in
Then, interpolation function lengths LiL and LiH of the low pass filter and the high pass filter are determined from the standard transition area ratios RtsL and RtsH, respectively. If a band pass filter according to specifications shown in
For the transition area of the high pass filter, the clock corresponding to 5 MHz is T5M=64 and the clock corresponding to 8.5 MHz is T8.5M=109. The clock width of the transition area to be designed is thus LdH=T8.5M−T5M=45. Therefore, the interpolation function length LiH of the high pass filter is also determined to be LiH=LdH/RtsH=69.189149.
To reduce the number of the taps in a band pass filter, it is desirable that the interpolation function lengths LiL and LiH be even integers larger than the calculated values. Thus, in this case, the interpolation function lengths LiL and LiH are both set to 70.
An interpolation function I (LPFB) for the low pass filter which has an interpolation function length LiL of 70 clocks is determined as shown in the following partition equations (Equation 11-1) and (Equation 11-2).
I (LPFB)=8/16+9/16*cos(2πt/70)−1/16*cos(6πt/70) (0/10245≦t≦69/1024) (Equation 11-1)
I (LPFB)=0 (69/1024<t≦1023/1024) (Equation 11-2)
An interpolation function I (HPFB) for the high pass filter which has an interpolation length LiH of 70 clocks is determined as shown in the following partition equations (Equation 12-1) and (Equation 12-2).
I (HPFB)=8/16−9/16*cos(2πt/70)+1/16*cos(6πt/70) (0/1024≦t≦69/1024) (Equation 12-1)
I (HPFB)=1 (69/1024<t≦1023/1024) (Equation 12-2)
Once the interpolation functions I (LPFB) and I (HPFB) for the low and high pass filters are thus determined, their frequency characteristics are shifted in the direction of the frequency axis (clock direction) so that the shifted interpolation functions I (LPFB) and I (HPFB) join the amplitude values “1” and “0” together. Specifically, 70 numerical sequences corresponding to the positions of the clock t=0/1024 to 69/1024 determined by (Equation 11-1) are shifted to the positions of the clock t=i/1024 to (i+69)/1024 (i is an integer). Further, 70 numerical sequences corresponding to the positions of the clock t=0/1024 to 69/1024 determined by (Equation 12-1) are shifted to the positions of the clock t=j/1024 to (j+69)/1024 (i>j; j is an integer). Then, all the numerical sequences at the positions of the clock t=1/1024 to (j−1)/1024 and (i+70)/1024 to 1023/1024 are changed to “0”. And all the numerical sequences at the positions of the clock t=(j+70)/1024 to (i−1)/1024 are changed to “1”.
Then, the frequency characteristic shown in
The interpolation function shift amounts i and j in this case may be set at such values as locate the amplitude value “0.5” of the interpolation function at positions on the frequency axis corresponding to its ⅛, 2/8, and ⅜. This simplifies the filter coefficients obtained by executing an inverse FFT on the numerical sequences for the input frequency characteristic in step S2 in
<Second Generation Method>
To generate an input frequency characteristic according to the second generation method, first, the standard functions XF1 and XF4, such as those shown in (Equation 1) and (Equation 6), are input to the circuit. After inputting the numerical sequences of the standard functions XF1 and XF4, an interpolation function containing frequency shifts is determined on the basis of the standard functions XF1 and XF4.
To determine an interpolation function subjected to frequency shifting, first, the requested transition area ratios RtrL and RtrH of the interpolation function corresponding to the standard functions XF1 and XF4 are determined. To determine the requested transition area ratio Rtr, information on two representative points (for example, the points at the amplitudes of −0.3 dB and −45 dB) in the transition area is used.
If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB in the standard functions XF1 and XF4 is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the value for the standardization clock corresponding to these amplitude values in the first half of the frequency characteristic shown in each of
Then, a calculation is made of the clock count LhsH from the standardization clock t=0 for the requested digital filter to the start point t=k1 of the transition area (see
Then, a calculation is made of the clock count LheH from the standardization clock t=0 for the requested digital filter to the end point t=k4 of the transition area. The clock count from the start point k1 to end point k4 of the transition area is defined as Tk1−k4. Then, the clock count LheH from the standardization clock t=0 to the transition area start point t=k4 is determined by LheH=LhsH+Tk1−k4. Here, the clock count Tk1−k4 from the start point k1 of the transition area to the end point k4 is determined by 0.5*RtrH using the clock count (=0.5) from the transition area start point (at t=0) to end point (at t=511/1024) in the standard function XF4, as well as the requested transition area ratio RtrH. As described above, RtrH=0.126963, so that Tk1−k4=0.5*0.126963=0.063482. Therefore, the clock count LheH from the standardization clock t=0 to the transition area end point t=k4 is determined to be LheH=0.097715+0.063482=0.161197.
Then, a calculation is made of the clock count LhsL from the standardization clock t=0 for the requested digital filter to the start point t=k5 of the transition area. The clock count from the transition area start point k5 to the point k2′ at −0.3 dB is defined as Tk5−k2′. The standardization clock at the point k2′ at −0.3dB is defined as Tk2′. Then, the clock count LhsL from the standardization clock t=0 to the transition area start point t=k5 is determined by LhsL=Tk2′−Tk5−k2′. Here, according to the filter specifications shown in
Then, a calculation is made of the clock count LheL from the standardization clock t=0 for the requested digital filter to the end point t=k6 of the transition area. The clock count from the start point k5 to end point k6 of the transition area is defined as Tk5−k6. Then, the clock count LheL from the standardization clock t=0 to the transition area end point t=k6 is determined by LheL=LhsL+Tk5−k6. Here, the clock count Tk5−k6 from the start point k5 of the transition area to the end point k6 is determined by 0.5*RtrL using the clock count (=0.5) from the transition area start point (at t=0) to end point (at t=511/1024) in the standard function XF1, as well as the requested transition area ratio RtrL. As described above, RtrL=0.134658, so that Tk5−k6=0.5*0.134658=0.067329. Therefore, the clock count LheL from the standardization clock t=0 to the transition area end point t=k6 is determined to be LheL=0.185473+0.067329=0.252802.
As a result, the interpolation function I (BPF) is determined as shown in the following partition equations (Equations 13-1, 13-2, 13-3, 13-4, and 13-5).
Specifically, the interpolation function I (BPF) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 13) from 0/1024 to 1023/1024. A graph showing these 1024 numerical sequences is almost similar to that in
Then, the frequency characteristic shown in
As shown in
Thus, even when filter coefficients with values smaller than a predetermined threshold are discarded by a rounding process, most of the major filter coefficients remain, which determine the frequency characteristic. Consequently, the frequency characteristic is not substantially affected. Further, the frequency characteristic obtained by the filter designing method according to the present embodiment exhibits a very significant attenuation. Consequently, the desired attenuation amount can be obtained even with a slight decrease in the number of bits.
This enables a rounding process to sharply reduce the unwanted filter coefficients. For example, by dropping lower several bits of the filter coefficient to reduce the bits, it is possible to round all the filter coefficients with values smaller than the maximum value expressed only by the lower several bits, to “0” for discarding.
Thus, the present embodiment can reduce the filter coefficients by performing a rounding operation utilizing coefficient values. Consequently, windowing, as utilized in the prior art, is not necessarily required. As described above, the standard function input in step S1 has a finite-base impulse response. Thus, the number of filter coefficients designed on the basis of this standard function is originally smaller than that in the prior art. Accordingly, the standard function can be used as it is without the need for a rounding process. However, to further simplify the circuit, a rounding process is preferably executed to reduce the bits.
As shown in
As described above, the filter designing method according to the present embodiment enables the designing of a low pass filter, a high pass filter, and a band pass filter which exhibit an appropriate ripple characteristic. Further, the windowing operation is not necessarily required, and the taps can be sharply reduced without any windowing operation. Moreover, in a rounding operation, the filter coefficients are each multiplied by a factor of 2X, with the product rounded to an integer. Multipliers used can thus be reduced. This enables the filter circuit to be implemented as an integrated circuit with a small area. Further, the determination of a filter coefficient of a desired frequency characteristic requires almost no trial and error processes, allowing an FIR filter to be easily designed.
Now, description will be given of the case in which a low pass filter according to the specifications shown in
XS1=1/2+1/2*cos(2πt) (Equation 14)
Either the above first or second generation method is also applicable to the determination of an interpolation function from the COS function XS1. Now, description will be given of a method for determining an interpolation function according to the first generation method as a representative.
To generate an interpolation function according to the first generation method, first, the standard transition area ratio Rts of the COS function XS1 is determined. If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the value for the standardization clock Td corresponding to these amplitude values in the first half of the frequency characteristic of the COS function XS1 shown in
Then, an interpolation function length Li is determined from the standard transition area ratio Rts. If a low pass filter according to specifications shown in
To reduce the number of the taps in a low pass filter, it is desirable that the interpolation function length Li be an even integer larger than the calculated value. Thus, in this case, the interpolation function length Li is set to 52. An interpolation function II (LPF1) with an interpolation function length Li of 52 clocks based on the COS function XS1 is determined as shown in the following partition equations (Equation 15-1) and (Equation 15-2).
II (LPF1)=1/2+1/2*cos(2πt/52) (0/1024≦t≦51/1024) (Equation 15-1)
II (LPF1)=0 (51/1024<t≦1023/1024) (Equation 15-2)
Specifically, the interpolation function II (LPF1) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 15-1 ad 15-2) from 0/1024 to 1023/1024.
Once the interpolation function II (LPF1) is thus determined, its frequency characteristic is shifted in the direction of the frequency axis (clock direction) so that the shifted interpolation function II (LPF1) joins the amplitude values “1” and “0” together. Specifically, 52 numerical sequences corresponding to the standardization clock t=0/1024 to 51/1024 determined by (Equation 15-1) are shifted to the positions of the clock t=i/1024 to (i+51)/1024 (i is an integer). The numerical sequences at the positions of the clock t=0/1024 to (i−1)/1024 are changed to “1”. And all the numerical sequences at the positions of the standardization clock t=(i+52)/1024 to 1023/1024 are changed to “0”.
Then, the frequency characteristic expressed by the numerical sequences thus generated is made laterally symmetric with respect to the position of the clock t=0.5. Specifically, the arrangement of all the numerical sequences except the one corresponding to the standardization clock t=0/1024, that is, the numerical sequences corresponding to the clock t=1/1024 to 511/1024, is reversed. The reversed numerical sequences are copied to the positions of the clock t=512/1024 to 1023/1024. The resulting laterally symmetric 1024 numerical sequences are determined to be the numerical sequences for the input frequency characteristic in step S1 in
Also in this case, the interpolation function shift amount i may be set at such a value as locates the amplitude value “0.5” of the interpolation function at positions on the frequency axis corresponding to its ⅛, 2/8, and ⅜. This simplifies the filter coefficients obtained by executing an inverse FFT on the numerical sequences for the input frequency characteristic in step S2 in
Such a spline function as shown below can also be used as a standard function.
XS2=1-2t2 (0/1024≦t≦511/1024) (Equation 16-1)
XS2=2 (t−1)2 (511/1024<t≦1023/1024) (Equation 16-2)
Further, a spline function XS3 expressed by (Equations 17-1 and 17-2) shown below may also be used as a standard function.
XS3=1-8t2 (0/1024<t≦255/1024) (Equation 17-1)
XS3=8 (1/2−t)2 (255/1024<t≦511/1024) (Equation 17-2)
As shown in
Thus, even when filter coefficients with values smaller than a predetermined threshold are discarded by a rounding process, most of the major filter coefficients remain, which determine the frequency characteristic. Consequently, the frequency characteristic is not substantially affected. Further, the frequency characteristic obtained by the filter designing method according to the present embodiment exhibits a very significant attenuation. Consequently, the desired attenuation amount can be obtained even with a slight decrease in the number of bits. This enables a rounding process to sharply reduce the unwanted filter coefficients. For example, by dropping lower several bits of the filter coefficient to reduce the bits, it is possible to round all the filter coefficients with values smaller than the maximum value expressed only by the lower several bits, to “0” for discarding.
Spline functions applicable to the present embodiment are not limited to the above examples. That is, using a spline function of a finite-base enables preferable results to be achieved as is the case with
A linear function can also be used as a standard function. If a linear function is used as a standard function, the standard function XL is expressed by the following partition equations (Equation 18-1) and (Equation 18-2).
XL=1−t/512 (0/1024≦t≦511/1024) (Equation 18-1)
XL=0 (511/1024<t≦1023/1024) (Equation 18-2)
Either the above first or second generation method is also applicable to the determination of an interpolation function from the linear function XL. Now, description will be given of a method for determining an interpolation function according to the first generation method as a representative.
To generate an interpolation function according to the first generation method, first, the standard transition area ratios Rts of the linear function XL is determined. If the amplitude value of the pass band is set to “1”, the amplitude value of −0.3 dB is 0.966051. The amplitude value of −45 dB is 0.005623. A calculation is made of the values for the standardization clock count Td corresponding to these amplitude values in the first half of the frequency characteristic of the linear function XL shown in
Then, an interpolation function length Li is determined from the standard transition area ratio Rts. If a low pass filter according to specifications shown in
To reduce the number of taps in a low pass filter, it is desirable that the interpolation function length Li be an even integer larger than the calculated value. Thus, in this case, the interpolation function length Li is set to 44. An interpolation function III (LPF) with an interpolation length Li of 44 clocks based on the linear function XL is determined as shown in the following partition equations (Equation 19-1) and (Equation 19-2).
III (LPF)=1−t/44 (0/1024≦t≦43/1024) (Equation 19-1)
III (LPF)=0 (43/1024<t≦1023/1024) (Equation 19-2)
Specifically, the interpolation function III (LPF) determined is 1024 numerical sequences calculated by varying the value of the clock t in (Equations 19-1 and 19-2) from 0/1024 to 1023/1024.
Once the interpolation function III (LPF) is thus determined, its frequency characteristic is shifted in the direction of the frequency axis (clock direction) so that the shifted interpolation function III (LPF) joins the amplitude values “1” and “0” together. Specifically, 44 numerical sequences determined by (Equation 19-1) and corresponding to the standardization clock t=0/1024 to 43/1024 are shifted to the positions of the clock t=i/1024 to (i+43) /1024 (i is an integer). And all the numerical sequences at the positions of the clock t=0/1024 to (i−1)/1024 are changed to “1”. And all the numerical sequences at the positions of the standardization clock t=(i+44)/1024 to 1023/1024 are changed to “0”.
Then, the frequency characteristic expressed by the numerical sequences thus generated is made laterally symmetric with respect to the position of the clock t=0.5. Specifically, the arrangement of all the numerical sequences except the one corresponding to the standardization clock t=0/1024, that is, the numerical sequences corresponding to the clock t=1/1024 to 511/1024, is reversed. The reversed numerical sequences are copied to the positions of the clock t=512/1024 to 1023/1024. The resulting laterally symmetric 1024 numerical sequences are determined to be the numerical sequences for the input frequency characteristic in step S1 in
Also in this case, the interpolation function shift amount i may be set at such a value as locates the amplitude value “0.5” of the interpolation function at positions on the frequency axis corresponding to its ⅛, 2/8, and ⅜. This simplifies the filter coefficients obtained by executing an inverse FFT on the numerical sequences for the input frequency characteristic in step S2 in
As shown in
In general, the out-of-band attenuation amount of the filter is limited by the number of bits that can be supported by hardware to be implemented. Accordingly, if there are no limitations on hardware scale, an out-of-band attenuation characteristic with more significant attenuation can be provided by increasing the number of bits x resulting from a rounding process. If a filter is designed using the standard function shown in (Equation 1), even when a filter coefficient resulting from a rounding process is made up of 16 bits, the number of taps hardly increases and the out-of-band attenuation amount of the frequency characteristic can be increased above −45 dB.
In contrast, when the COS function, the spline function, or the linear function is used as a standard function, the number of taps required increases consistently if the number of bits x of a filter coefficient resulting from a rounding process is increased. However, reducing the number of bits in a filter coefficient to some degree enables the number of required taps to be reduced to a number equivalent to that achieved with the function shown in (Equation 1). Consequently, under the conditions under which the number of bits can be reduced to some degree by a rounding operation, a filter designing method can be effectively used which uses the COS function, the spline function, or the linear function as a standard function.
With a 12-bit accuracy, which is often used in the field of digital filters, no marked difference occurs in tap count among the function shown in (Equation 1), the COS function shown in (Equation 14), and the linear function shown in (Equations 18-1 and 18-2). Thus, the filter designing method is effective regardless of whichever function is used as a standard function.
An apparatus for realizing the digital filter designing method according to the above present embodiment can be implemented using a hardware configuration, a DSP, or software. If the filter designing apparatus according to the present embodiment is to be implemented, for example, by software, it is actually composed of a CPU, or an MPU, a RAM, or a ROM in a computer. In this case, the apparatus can be implemented by operating a program stored in the RAM or ROM or a hard disk.
For example, the functions of a spreadsheet installed in a personal computer or the like can be used to execute inputting of a standard function, calculation of an interpolation function from the standard function, frequency shifting of the interpolation function, an inverse FFT operation on an input frequency characteristic generated by shifting the interpolation function, a rearrangement operation on numerical sequences, a rounding operation, and the like. In this case, the operations are actually performed by the CPU, ROM, RAM, or the like in the personal computer or the like in which the spreadsheet is installed.
The filter designing apparatus according to the present embodiment may comprise a table information storage section that stores table information on attenuation values starting with 0 dB which are associated with standardization clock values corresponding to the attenuation values on a predetermined standard function, an input device that inputs information on requested specifications as shown in FIGS. 46 to 48, and a calculation device that calculates an interpolation function using the input information on the required specifications and the above table information. With this configuration, simply by inputting the information on the required specifications to the circuit using the input device, it is possible to automatically determine an interpolation function (input frequency characteristic to be subjected to an inverse FFT) that meets the required specifications.
Alternatively, filter coefficients determined may be automatically subjected to FFT, with the results shown on a display screen as a frequency characteristic diagram. This enables the frequency characteristic of the designed filter to be visually checked. Therefore, filter designing can be more easily achieved.
To actually implement a digital filter in an electronic device or on a semiconductor IC, an FIR filter may be constructed so as to have, as filter coefficients, numerical sequences finally determined by the filter designing apparatus as described above. Specifically, as shown in
As described above in detail, the present embodiment inputs a standard function to the circuit and calculates an interpolation function from the standard function to determine an input frequency characteristic. Then, an inverse FFT is executed on a numerical sequence indicative of the input frequency characteristic to determine filter coefficients. Consequently, the coefficients for an FIR digital filter which allow a desired frequency characteristic to be provided can be easily determined without any special mathematical or electrical engineering knowledge. In particular, the same technique can be used to easily design not only a low pass filter but also a high pass filter, a band pass filter, a band elimination filter, or a comb filter.
Further, the present embodiment does not necessarily require a windowing operation for reducing the number of filter coefficients. Instead, a numerical rounding operation enables the number of filter coefficients to be reduced without lowering the accuracy of the frequency characteristic. The present embodiment can also simplify the filter coefficient values by performing an arithmetic operation on a numerical sequence determined by an inverse FFT to obtain integers. This makes it possible to sharply reduce required multipliers, filter components, to simplify the filter configuration. Furthermore, a desired frequency characteristic can be accurately provided.
The above embodiment has been described in conjunction with the usage of the standard functions XF1 to XF6, XS1 to XS3, and XL. However, standard functions that can be used for the present invention are not limited to them.
The above embodiment has been described in conjunction with the process of multiplying a numerical sequence by a factor of 2X and rounding down the resulting decimal fractions, as an example of an operation for making filter coefficients integers. However, the present invention is not limited to this. For example, the numerical sequence may be multiplied by a factor of 2X with the resulting decimal fractions rounded up or off.
As another example in which filter coefficients are made integers, a numerical sequence of filter coefficients may be multiplied by a factor of N (N is a value other than the power of 2), with the resulting decimal fractions rounded (rounded down, up, or off). To perform a rounding operation on a numerical sequence multiplied by a factor of N, the digital filter can be configured as shown in
Further, multiplying the numerical sequence by a factor of 2X enables the filter coefficients to be rounded in bits. In contrast, multiplying the numerical sequence by a factor of N enables inter-bit rounding to be executed on the filter coefficients. The rounding process in bits refers to a process of making the coefficient values integral multiples of 1/2X; if for example, the coefficient values are multiplied by a factor of 2X, with the resulting decimal fractions rounded down, then all the numerical values belonging to the range from 2X to 2x+1 are rounded to 2X. Further, the inter-bit rounding process refers to a process of making the coefficient values integral multiples of 1/N; if for example, the coefficient values are multiplied by a factor of N (for example, 2x−1<N<2X), with the resulting decimal fractions rounded down, then all the numerical values belonging to the range from N to N+1 are rounded to N. The rounding operation on the coefficient values multiplied by a factor of N makes it possible to adjust the filter coefficients to be made integers to arbitrary values other than the power of 2. This enables the precise adjustment of the number of filter coefficients (number of taps) that are used for the digital filter.
Alternatively, as an example of a rounding operation with filter coefficients made integers, all y-bit filter coefficients with a data value of smaller than 1/2X may be determined to be zero, while for y-bit filter coefficients with a data value of at least 1/2X, the data values may be multiplied by factor of 2x+X (x+X<y), with the resulting decimal fractions rounded (rounded down, up, or off).
To perform such a rounding operation as described above, the digital filter can be configured as shown in
Further, by considering all data values of smaller than 1/2X to be zero to round them down, it is possible to sharply reduce the number of required filter coefficients (taps) and to determine accurate filter coefficients composed of (x+X) bits, that is, more than x bits. This enables a more appropriate frequency characteristic to be obtained.
Further, the above embodiments are only specific examples for carrying out the present invention and are not intended to limitedly interpret the technical scope of the present invention. That is, the present invention can be carried out in various manners without departing from its sprit or major characteristics.
The present invention is useful for an FIR digital filter of a type that comprises a tapped delay line made up of a plurality of delayers and which multiplies output signals from the taps of the tapped delay line by respective filter coefficients, with the multiplication results added up for output. The present invention is also useful for a method for designing this FIR digital filter.
Number | Date | Country | Kind |
---|---|---|---|
2004-123245 | Apr 2004 | JP | national |
Number | Date | Country | |
---|---|---|---|
Parent | PCT/JP04/17667 | Nov 2004 | US |
Child | 11550345 | Oct 2006 | US |