JPH04177461A - High speed technical calculation method - Google Patents

Abstract

PURPOSE:To perform a technical calculation at a high speed without using any expensive numerical arithmetic processor and at the same time to reduce the necessary capacity of a memory byperforming a numerical operation through an integer arithmetic and in a block floating point format. CONSTITUTION:A fixed floating point arithmetic using an integer arithmetic is expressed in a block floating point number that showed the data in an exponent part or a coefficient part common to all data and a fixed decimal point of each data. Then a block floating point arithmetic applying an integer arithmetic is carried out. It is supposed that a decimal, point is set at the 13th bit of an integer and then that 8192=2<13> of an integer is equivalent to 1.01 of a fixed decimal point. In this fixed decimal point format, the value covering -4.0 (minimum) through about 3.9999 (maximum) can be expressed. The expressible number most approximate to 0 is about 0.000122 and a number approximate to 1.0 has the decimal accuracy equivalent to about 4 digits. The diagram shows a program used for calculation of a complex number end applying a fixed decimal point in a C language.

Description

DETAILED DESCRIPTION OF THE INVENTION [Field of Industrial Application], 1: (7) The invention relates to speeding up technical calculations by a computer.

[Conventional technology]

FIGS. 5(a) and 5(b) are complex number arithmetic programs in the C language published on pages 209 to 211 of "Scientific and Technical Computation in C" by Shin Koike, published by CQ Publishing.

(1) defines the format complex that represents complex numbers, and converts the real part of complex numbers into double-precision floating point (hereinafter referred to as do
(abbreviated as uble), the imaginary part is double
It is raining at 1m of e. (2) is the complex type addition function cadd defined in (1), and (3
) is a complex type subtraction function csul, (4)
is a function cmul for multiplication of complex type, and (
5) is a complex type division function cdiv. (6) is a function cabs that calculates the absolute value of a complex type number, and (7) is a function to- that calculates a complex type number where each component is two double type numbers.
It is complex.

Figure 6 is an FFT program using complex number operations that was published on page 218 of the same book, ``Scientific and Technical Calculations Using rc''. (10) is the declaration part of the subroutine name cfft, (11) is the beginning of the pit reversal operation part, (12
) indicates the beginning of the butterfly ridge.

FIGS. 7(a) and 7(b) are practical complex FFT programs based on FIG. 6, with the addition of the split wing section after FFT calculation, which was omitted in the program of FIG. (10)～
(12) is the same as in Figure 6, and (13) is the FFT
It shows the division part that divides the elements of the array Z after the calculation by n. FIG. 8 is a flowchart of the pit reversal section, FIG. 9 is a flowchart of the butterfly ridge section, and FIG. 10 is a flowchart of the split wing section.

Next, the operation will be explained.

The function cadd in (2) of FIG. 5 is a complex variable a.

The receiver takes b as an argument and returns the complex number Z, which is the sum, Z, re (real part of Z) = a, re (real part of a) +
b, re (real part of b) Z, im (imaginary part of Z) = a
, im (imaginary part of a) + b, im (imaginary part of b) are calculated.

The function csub in (3) in Figure 5 similarly returns the difference of complex numbers, and X, re = a, re - b, re, X, im
=a, im-b, im is being calculated.

The function cmul in (4) in Figure 5 similarly returns the product of complex numbers: X, re=a, reX b, re -a, imX b,
imX, im=a, reX b, 1m -a, imX
b, re are being calculated.

The function cdiv in (5) of FIG. 5 similarly returns the quotient of a complex number, and calculates the following. Also, when the denominator is O, an error is displayed.

The function cabs in (6) in Figure 5 doubles the absolute value of a.
The function returned by le calculates l a l = (a, re) ``10 (a, im).

The function tocolIlplex in (7) in Figure 5 is dou
A function that takes the variables X and y of ble as arguments, converts them to a complex number Z, and returns Z, re=x, Z, im=
Executing y.

FIGS. 6 and 7 are examples of FFT programs using complex number operations.

FIG. 8 is a flowchart of the pit reversal portion of the above program, and explains the part (11) in FIG.

FIG. 9 is a flowchart of the butterfly calculation section of the program, and explains the part (12) in FIG. 7.

FIG. 10 is a flowchart of the split wing section (13) in FIG. 7 (omitted in FIG. 6).

In the FFT program cfft shown in FIG. 7, an array Z of complex number data, an integer 1 . An integer variable f that selects FFT and RFT (if f = 1, RFT, f--1
If f-1, then FFT is executed) is given as an argument and executed, and if f-1, data obtained by inverse Fourier transform of the data in array Z is obtained in array Z. If f=-1, data obtained by Fourier transforming the data in array Z is obtained in array Z.

The FFT program is a program that performs a discrete Fourier transform at high speed, and is a program that calculates the frequency spectrum G (n) of the discrete number data g (k) that changes over time using the following relationship. However, N is the number of data and is usually expressed as a power of Z (2, 4, 8, ..., 128,
256,512,1024. ), and WN is 6-
It is a number expressed by JM.

FFT programs are mainly used to determine the frequency components of electrical and mechanical signals, and in the past were only used on large computers, but in recent years they have been used on personal computers and measuring instruments. It is also used above.

If a complex number FFT is described using the complex number operations (2) to (7) in FIG. 5, a short program as shown in FIG. 7 will be obtained.

When the program shown in FIG. 7 was executed on a personal computer PC9801RX2 using Borland's C compiler Turbo-C, it took about 13 seconds when the number of data N = 1024 and without a numerical calculation processor.

If a numerical calculation processor 180287-10 is attached, the same calculation can be executed in about 3 seconds. However, 180287-10
Currently, it is expensive at about 40,000 to 60,000 yen.

[Problem to be solved by the invention]

Conventional technical calculation methods have been performed using double-precision floating point arithmetic supported by the C language, as described above. In double-precision floating-point operations, each data is expressed using a 64-bit exponent system, so a high precision of 15 digits can be obtained with decimal numbers, but on the other hand, even basic operations such as addition, subtraction, multiplication, and division can be performed using 100 digits of integer operations. It takes ~1000 times as much time and requires 8 bytes of memory for 1 data. For this reason, conventional technical calculation programs take a long time and require a large amount of memory. Numerical arithmetic processors exist to speed up floating point arithmetic operations, but they are quite expensive and are therefore not often implemented in general-purpose personal computers.

This invention was made to solve the above problems, and is a high-speed technology that can obtain results in a short time without using an expensive numerical processor and requires less memory. The purpose is to obtain a calculation program.

[Means to solve the problem]

The high-speed technical calculation method according to the present invention represents data in fixed point or block floating point, and performs calculations in fixed point arithmetic using integer arithmetic or block floating point.

[Effect]

The high-speed technical computation method of this invention uses fixed-point arithmetic using integer arithmetic instead of double-precision floating-point arithmetic, which has higher precision than necessary, requires a large amount of memory, and takes a long time, or The exponent part or coefficient part common to the data and blocks expressed as fixed point numbers for each piece of data are expressed in floating point numbers, and by performing block floating point operations using integer operations, it is possible to speed up technical calculations and to store data. Reducing the amount of memory required.

〔Example〕

An embodiment of this invention will be described below.

In this embodiment, fixed-point numbers are expressed using integers. The representation of integers differs depending on the computer, but for example, the C for personal computer PC: 9801 series
In the compiler Turbo-Cver 2.0 (Borland), integers are 16-bit signed integer types (hereafter i
(abbreviated as nt) is the standard. A 16-bit signed integer can represent 21'' numbers, with a minimum of -32768~
A maximum of 32767 can be expressed.

In this embodiment, it is assumed that there is a decimal point at the 13th bit of this integer, and the integer 8192=213 is set to correspond to the fixed decimal point 1.01. This fixed-point format can represent a minimum of -4.0 to a maximum of about 3.9999. The closest number to O that can be represented is approximately 0,000
122, and a number around 1.0 has an accuracy of about 4 digits in decimal.

This fixed-point operation can be performed in exactly the same way as an integer in the case of addition and subtraction, but in the case of multiplication, 1. OXl, 0=1.0
It is necessary to adjust so that Furthermore, the product of 16-bit integers reaches a maximum of about 1 billion, which cannot be expressed by 16-bit signed integers. Therefore, in multiplication, it is necessary to convert each integer into a 32-bit signed integer type (hereinafter abbreviated as long int).

A 32-bit signed integer can represent numbers from a minimum of about -200 million to a maximum of about 2.1 billion. 81 after multiplication by long int
Multiplication is possible by dividing by 92 and converting it back to int type. For example, fixed-point number a = 2.0. The operation of assigning the product of b = 0.5 to fixed point C is expressed in C language as follows: tnt a, b, c; a=(fnt) (2,0m8192.0) ;b=(
fnt)(0,518192,0);c=(int)(
(long int)as(long 1nt)b
/8192); However, 1nt) is a cast operator that instructs type conversion to int. Also (long
1nt) is a cast operator that instructs type conversion to long int. Fixed-point division is also lon
Need to calculate with g int 9 For example a = 1.0
．． To assign b divided by 8 by b = 2.0 to C, int a, b, c; a = (int) (1, O * 8192.O); b = (in
t) (2,0 zero 8192.O);c=(int)((l
long 1nt) a book 8192/(long 1n
t) b) Can be written in H and C languages.

In this way, fixed-point operations can be performed using a combination of integer operations.

Furthermore, it is also possible to represent a complex number with two fixed point numbers. In other words, the real part of a complex number can be represented by one fixed-point number, and the imaginary part by another fixed-point number. If this is defined as bcotnp type in C language, it will be as follows.

typedef 5truct(int x; int
y;)bcompHbcomp Zl, Z2. Z3; In this case, the real part of the complex number Z1 is Zl, x, and the imaginary part is zi,
Corresponds to y. Z1 Tomi 0.7 + J0.7. Z2=0.
To assign the sum of 5+J O,8 to Z3, declare it using the above definition formula, then Zl, x=(int) (0,7 zero 8192.0)
;zl,y=(int)(0,7 8192.0);Z
2. x=(int)(0,5 zero 8192,0)iZ2.
y=(int) (0,8 lines 8192.0): Z3. x
=Z1. x+Z2. x; Z3. y=Z1. y+Z2. You can write it as y;

As you can see, it would be difficult to read and write the conversion from a floating point number to a complex number using a fixed point number or the calculation of a complex number as it is, so
It is preferable to do this as a C language function.

FIGS. 1(a) and 1(b) are programs for complex number calculations using fixed point numbers in C language. The first declaration section is exactly the same as described above.

The function bcadd in (2) is an addition function, and the complex part of the function cadd in the conventional embodiment is
It has been changed to mp.

The function bcsub in (3) is a subtraction function, and this is also a function in which the complex part of the function csub in the conventional embodiment is changed to bcomp.

The function bcmul in (4) is a function that calculates a product, and the complex part of the function cmul in the conventional embodiment is
omp, and then calculate the multiplication of the real and imaginary parts. On the left side, convert the integer type from int to long int, perform the multiplication, and return it to a fixed point, so divide by 8192 and change the type back to int. It is rewritten as follows.

The function bcdiv in (5) is a function that calculates the quotient, and coo+plex of the function cdiv of the conventional embodiment is
p and performs integer type conversion. Since it is a fixed-point quotient, you can multiply the dividend by 8192 and then divide by the divisor, but the numerator is the sum of the products of two fixed-point numbers, for example (2, O + J0.0) / (2, 0+
J.O.

0), the real part of the numerator can be written as an integer (163
84) 2+0=268435456, and 8
When multiplied by 192, it becomes approximately 2.2X10” long i
It cannot be expressed using nt. The maximum value of long 1nt is approximately 2X10", which is 8 times the 268435456 in the previous example, so when calculating the quotient, for example, the numerator can be up to 8 times as long as the absolute value of the dividend and divisor is 2 or less. Therefore, the numerator is multiplied by 8 and the denominator is 8192÷8.
Divide by -1024, then perform the division, then lon
g Changed from 1nt back to int type. Therefore, the absolute value of the divisor is an integer and fi = 32, and the fixed point is 0.0.
If it is smaller than 04, an error will occur. At this time, an error message is displayed.

To further speed up the division, it is possible to stop multiplying the numerator and divide the denominator by 8192, but in this case, if the absolute value of the divisor is an integer and is less than 0.01 in fi#90 fixed point, an error will occur.

To increase the precision of division, for example, set the absolute value of the dividend and divisor to be less than or equal to 1, and calculate 231÷(819
2) Multiply the numerator by 32 with a margin of ”=32, and the denominator by 2.
All you have to do is divide by 56 and then do the division.

In this case, the absolute value of the divisor is an integer, 4r256=16. An error occurs when the value is less than 0.002 in a fixed point number.

In FIG. 1, (6) is a function bcabs for calculating the absolute value of a complex number. C01D of cabs in conventional embodiment
-plex is replaced with bcomp.

In Figure 2, (7) is a double whose absolute value is 2 or less.
This is a function tabcoa+p that calculates the number of bcomp whose real and imaginary parts are the two numbers of e, and is almost the same as tocomplex in the conventional embodiment, but is multiplied by 8192 to convert it to a fixed-point format.

(8) is a subroutine bcnorm that divides the n elements of the array Z[) of complex numbers bcompO by 2, and is used in the FFT program described later.

(9) is a subroutine bcnorm2 that adjusts n elements of the array Z() of bcomp so that the maximum absolute value of each component (real part, imaginary part) is 1 or less.

In order to make the maximum value of M consideration less than or equal to 1, all elements of Z() are divided by the variable con8t, and the variable amp corresponding to the common exponent part of the array Z[] is multiplied by canst. This subroutine supports the block floating point format consisting of Z and amp.

(10) is a complex FFT using the complex operations mentioned above.
It is a program. The algorithm is almost the same as the conventional example, but in the conventional example, the division by the number of data N was performed all at once at the end of the program, but in this program, the data is divided into a subroutine each time a set of butterfly product operations is performed. By dividing by 2 using bcnorm,
Executing division by N. In the FFT program, each time a set of butterfly products is performed, the maximum value of the data is approximately 2.
Since it has the property of multiplying, by dividing by 2 each time to prevent the fixed point from overflowing midway,
A division by N is executed.

Flowcharts explaining different parts of the algorithm of this complex FFT program are shown in FIGS. 11 and 12. In Figure 11, (13) is the added bcno
This shows the execution position of the rm subroutine. Instead of proceeding to the division by N program indicated by ■ in the conventional program of FIG. 9, in FIG. 11 the subroutine is ended and the program returns to the main program.

FIG. 12 is a flowchart showing the algorithm of the bcnorm subroutine, in which each element of array Z is divided by two. In this way, addition, subtraction, multiplication, division, and complex number FFT can be realized using complex number bcmp in fixed-point decimal format.

Figure 4 shows the main program in C language created to test the above complex number FFT program.
Set the real parts of the four data to draw a sine curve,
The imaginary parts are all set to 0, and a complex number FFT program is executed.

In other words, we are looking for the frequency spectrum when inputting data that appears to be oscillating sinusoidally on the time axis.

The subroutine prJime is a program that displays the time, and is used to calculate the time consumed by the 20FFT program.

Results of running this main program and measuring the time 1
FFT of 024 points was completed in about 2 seconds.

Approximately 6 times faster speed than the conventional example was achieved. It is a little faster than one using a numerical processor, but the accuracy is about 2 to 3 digits. For example, 8-bit A/
This accuracy seems to be sufficient for processing data input by a D converter or the like. Also, the results are shown on the screen (4
This level of accuracy is sufficient even for an application program that displays a graphic (00 dots x 640 dots).

Although the program was written in the C language in the above embodiment, a similar program can be created in other structured high-level languages such as PASCAL.

Further, in the above embodiment, a program for executing FFT was written, but any technical calculation program that processes a plurality of similar data can be speeded up using the same method.

Also, in the above example, the FFT is calculated using complex number calculations, but it is also possible to separate the data of the real part and imaginary part of the complex number, arrange them in separate arrays, and calculate them using normal real number calculations. It is.

It is also possible to perform various technical calculations using the complex number representation and complex number operations of the data used in the above example, but depending on the type of technical calculation, the precision of the fixed point may be 10.
Since it is approximately 4 digits in base, care must be taken as a large error may occur during calculation.

〔Effect of the invention〕

As described above, according to the present invention, numerical calculations are performed in block floating point format using integer calculations, making it possible to perform technical calculations at high speed without using an expensive numerical calculation processor. The memory capacity can also be reduced.

[Brief explanation of drawings]

FIG. 1 is a diagram showing a function method for adding, subtracting, multiplying, and dividing complex numbers and determining the absolute value according to an embodiment of the present invention, FIG. 2 is a diagram showing a method for processing complex numbers according to an embodiment of the present invention, and FIG. FIG. 4 is a diagram showing a complex FFT according to an embodiment of the present invention, FIG. 5 is a diagram of a conventional complex FFT function, and FIG. 6 is a diagram showing a conventional complex FFT function. F
FT diagram, Figure 7 is a conventional practical complex FFT diagram, Figure 8
FIG. 9 is a flowchart of the pit reversal part of the complex number FFT program common to the conventional and one embodiment of the present invention.
10 is a flowchart of the division section of a conventional complex number FFT program. FIG. 11 is a flowchart of the butterfly seed section of a complex number FFT program according to an embodiment of the present invention. FIG. 12 shows the subroutine BCNO of a complex number FFT program according to an embodiment of the present invention.
It is a flowchart of RH. (1) is a complex number structure declaration, (2) is a complex number addition function, (3) is a complex number subtraction function, (4) is a complex number multiplication function, (5) is a complex number division function, (6) is a function that calculates the absolute value of a complex number, (7) is a function that creates a complex number from two real numbers, (8) is a subroutine that divides n elements of an array of complex numbers by 2, and (9) is an array of complex numbers. A subroutine that makes the maximum value of each component of the elements smaller than 1, (1
0) is a complex number FFT subroutine, (11) is a pit reversal section, (12) is a butterfly seed section, and (13) is a division section. In addition, in the figures, the same reference numerals indicate the same or corresponding parts.

Claims (1)

[Claims] A high-speed technical computation method characterized by using blocked floating-point arithmetic or fixed-point arithmetic using integer arithmetic in technical computation that processes multiple pieces of data on a structured general-purpose high-level language.