A Recipe for a Basic Trigonometry Table
Introduction
This is an article that is give a better understanding to the Discrete Fourier Transform (DFT) by showing how to build a Sine and Cosine table from scratch. Along the way a recursive method is developed as a tone generator for a pure tone complex signal with an amplitude of one. Then a simpler multiplicative one. Each with drift correction factors. By setting the initial values to zero and one degrees and letting it run to build 45 values, the entire set of values needed for a degee based trigonometric table are produced. Like my Logarithm Article [3], this one is meant to be calculable by hand, but actually won't be. More supportive methods are demonstrated and finally source code in C is provided to show the calculations actually work. References to my other articles are provided for more detailed explanations of some of the concepts.
A Point on the Complex Plane
A complex number lies on the complex plane and can be represented either a sum of a real number and an imaginary number, corresponding to a Cartesian representation, or in exponential form corresponding to polar coordinate representation. The key equation tying these two representations together is Euler's formula [1] multiplied by a length, denoted by $r$. $$ r e^{i\theta} = r \cos(\theta) + i \, r \sin(\theta) \tag {1} $$ This has a very geometric interpretation. $\theta$ is the distance along the complex unit circle measured in radians where the radial axis is positioned. The $r$ value then indicates where on the radial axis the designated point lies. By the triangle based trigonometric definitions, this point is at $(r cos(\theta), r \sin(\theta) )$ in Cartesian coordinates on the complex plane and represented as a complex number by the right hand side.
Multiplying Complex Numbers
Complex numbers are often denoted by using $z$ as a variable name. Suppose you have three of them. $$ \begin{aligned} z_1 &= r_1 e^{i\theta_1} \\ z_2 &= r_2 e^{i\theta_2} \\ z_3 &= r_3 e^{i\theta_3} \\ \end{aligned} \tag {2} $$ Suppose further that the third is the product of the first two. Following the usual rules of multiplication reveals some truths. $$ \begin{aligned} z_3 &= z_1 \cdot z_2 \\ &= r_1 e^{i\theta_1} \cdot r_2 e^{i\theta_2} \\ &= ( r_1 \cdot r_2 ) e^{i(\theta_1 + \theta_2)} \\ &= r_3 e^{i\theta_3} \\ \end{aligned} \tag {3} $$ By taking the magnitude of the last two lines: $$ r_3 = r_1 \cdot r_2 \tag {4} $$ By looking at the angles: $$ \theta_3 = \theta_1 + \theta_2 \tag {5} $$ The commonly heard statement "When you multiply two complex numbers, you add the angles and multiply the magnitudes." nicely summarizes these rules.
Complex Tone Generation
If the magnitude of a complex number is one, that means it lies on the unit circle. If the magnitudes of the values are both one then the magnitude of the product will also be one. This means the product of any two numbers on the complex unit circle will also be a number on the unit circle.
This makes generating a set of evenly spaced points along the complex unit circle simple. This just happens to be what a pure complex tone with unit magnitude is made of. Often, a convention that the angle of the initial point is designated by $\phi$, and the arc length of each step is $\omega$. The first few values readily show how it works. $$ \begin{aligned} z[0] &= e^{i\phi} \\ z[1] &= z[0] e^{i\omega} = e^{i\phi} e^{i\omega} = e^{i(\omega + \phi)} \\ z[2] &= z[1] e^{i\omega} = e^{i(\omega + \phi)} e^{i\omega} = e^{i(2 \omega + \phi)} \\ z[3] &= z[2] e^{i\omega} = e^{i(2\omega + \phi)} e^{i\omega} = e^{i(3 \omega + \phi)} \\ &: \\ z[n] &= z[n-1] e^{i\omega} = e^{i((n-1)\omega + \phi)} e^{i\omega} = e^{i(n \omega + \phi)} \\ \end{aligned} \tag {6} $$ The last line then is the rule that is the definition of a unit complex tone. $$ z[n] = e^{i(n \omega + \phi)} = \cos(n \omega + \phi) + i \, \sin(n \omega + \phi) = z_\omega^n z_\phi \tag {7} $$ To give it amplitude, multiply the signal with a scalar. That is not the concern of this article though.
Quantization and Drifting
The approach in the previous section works great theoretically, but when implemented numerically quantization errors will creep in. The more common expression is rounding errors. $$ z_\omega = e^{i\omega} = x + i \, y \tag {8} $$ When done with actual numbers, either when by doing it by hand with a finite number of digits, or on a computer with a set number of bits, the digits or bits that are neglected introduce a slight error that accumulates with repeated multiplications. This can be included in the equation as an error factor. $$ z_\omega = ( 1 + d ) e^{i\omega} \tag {9} $$ When this number is raised to a power via repeated multiplications, the error accumulates. $$ z_\omega^n = ( 1 + d )^n e^{i(n\omega)} \tag {10} $$ The parentheses in the exponent aren't necessary, but are included for conceptual clarity. Since the value of $d$ in this context is very small, an approximation can be used to estimate its value. $$ ( 1 + d )^n \approx 1 + nd \tag {11} $$ This approximation is apparent by looking at the first few lines of a binomial expansion. It is done more formally by invoking Pascal's Triangle, and the equations behind it. $$ \begin{aligned} ( 1 + d )^0 &= 1 \\ ( 1 + d )^1 &= 1 + d \\ ( 1 + d )^2 &= 1 + 2d + d^2 \\ ( 1 + d )^3 &\approx 1 + 3d + 3d^2 \\ ( 1 + d )^4 &\approx 1 + 4d + 6d^2 \\ ( 1 + d )^5 &\approx 1 + 5d + 10d^2 \\ &: \\ \end{aligned} \tag {12} $$ The higher powers of $d$ get smaller much faster than the coefficients get larger, so they can be neglected in an approximation.
The result of this accumulation of length error is that the results of the tone generation drift from having an amplitude of one.
Squaring a Unit Value
This drifting can be greatly mitigated when squaring a value. A property of complex numbers on the unit circle is that their reciprocal is also their conjugate. $$ \begin{aligned} z_\omega &= e^{i\omega} = x_\omega + i \, y_\omega \\ 1/z_\omega &= e^{-i\omega} = \frac{ x_\omega - i \, y_\omega }{x_\omega^2 + y_\omega^2} = z_\omega^* \\ \end{aligned} \tag {13} $$ For a point on the unit circle, the denominator $x^2 + y^2$ will equal one. However, suppose the values of $x$ and $y$ are off just a little bit. This might change the radius just a little bit, and the angle just a little bit. An adjustment term can be included to accommodate the radial error, and a new name for the angle $\sigma$, makes this the result: $$ \begin{aligned} z_\sigma &= ( x + i \, y ) = ( 1 + d ) e^{i\sigma} \approx ( 1 + d )z_\omega \\ 1/z_\sigma &= \frac{1}{ 1 + d } e^{-i\sigma} \approx \frac{1}{ 1 + d } z_\omega^* \\ \end{aligned} \tag {14} $$ Using regular algebra to do the division gives an approximation for the adjustment term. $$ \begin{aligned} \frac{1}{ 1 + d } &= 1 - d + d^2 - d^3 + d^4 - .... \\ &\approx 1 - d \\ \end{aligned} \tag {15} $$ This can then be applied to (14). $$ \begin{aligned} 1/z_\sigma &\approx ( 1 - d ) z_\sigma^* \approx ( 1 - d )( x - i \, y ) / ( x^2 + y^2 ) \\ z_\sigma^{**} &= ( z_\sigma^*)^* \approx ( 1/z_\sigma)^* \approx ( 1 - d )( x + i \, y ) / ( x^2 + y^2 ) \\ \end{aligned} \tag {16} $$ To make the mitigated square value, the two versions of $z_\sigma$ are multiplied together. $$ \begin{aligned} z_\omega^2 &\approx z_\sigma^2 \\ &\approx z_\sigma \cdot z_\sigma^{**} \\ &\approx z_\sigma \cdot ( 1/z_\sigma )^* \\ &\approx( 1 + d )( 1 - d )\left[ ( x^2 - y^2 ) + i \, (2xy) \right] / ( x^2 + y^2 )\\ &= ( 1 - d^2 ) \left[ ( x^2 - y^2 ) + i \, (2xy) \right] / ( x^2 + y^2 )\\ \end{aligned} \tag {17} $$ When $d$ is at the limit of the precision of the numbers being used, then $d^2$ can be neglected and this result essentially cancels out the drift error. $$ z_\omega^2 \approx \frac{( x^2 - y^2 ) + i \, (2xy) }{ x^2 + y^2 } \tag {18} $$ There is also a trick to convert the division into a multiplication. $$ \begin{aligned} g &= x^2 + y^2 - 1 \\ 1/( x^2 + y^2 ) &= 1/( 1 + g ) \\ &\approx ( 1 - g ) \\ &= 2 - (x^2 + y^2) \\ \end{aligned} \tag {19} $$ After all of that, the equation for the normalizing squaring operation can be expressed as: $$ z_\omega^2 \approx \left[ ( x^2 - y^2 ) + i \, (2xy) \right] \cdot \left[ 2 - ( x^2 + y^2 ) \right] \tag {20} $$ This will adjust the magnitude to be closer to one.
A Numerical Example
At this point, a numerical example might be useful to illustrate how it works. Starting with the Pythogorean triple of $ 3^2 + 4^2 = 5^2$, then dividing by 5, proves that the value $0.6 + i \, 0.8$ will lie exactly on the unit circle. However, suppose the value is intead close to that, say $0.61 + i \, 0.81$, what happens?
The straight up square is straightforward. $$ \begin{aligned} (0.61 + i \, 0.81)^2 &= ( 0.61^2 - 0.81^2 ) + i \, (2 \cdot 0.61 \cdot 0.81 ) \\ &= ( 0.3721 - 0.6561 ) + i \, (2 \cdot 0.4941) \\ &= -0.284 + i \, 0.9882 \\ \end{aligned} \tag {21} $$ Calculating its magnitude: $$ \| -0.284 + i \, 0.9882 \| = \sqrt{ (-0.284)^2 + (0.9882)^2 } \approx 1.0282 \tag {22} $$ Fortunately, most the calculations have already been done to calculate the mitigating factor. $$ 2 - ( x^2 + y^2 ) = 2 - ( 0.3721 + 0.6561 ) = 0.9718 \tag {23} $$ Multiplying this by the magnitude of the square shows the mitigated magnitude. $$ 1.0282 \cdot 0.9718 = 0.99920476 \tag {24} $$ This value is now very close to one. A precision error of 0.01 is much larger than those that occur at the limits of any reasonable number of digits or bits used in these calculations. The resulting error is very much in the ball park of the expected approximate value. $$ 1 - 0.01^2 = 0.9999 \tag {25} $$
Tweaking the Angle
The approximation in the previous sections is not able to adjust the angle since the desired angle is not necessarily known. The value of $\sigma$ can be found by taking the $\arg( z_\sigma^2 )$ and dividing the result by two. The $\arg$ function simply means "angle of". In programming languages, the equivalent is the atan2 function. $$ \arg( x + i \, y ) = \operatorname{atan2}( y, x ) \tag {26} $$ In contrast, taking the $\arctan$ of the value only works for two quadrants of the Cartesian plane.
The angle of the estimated square can then be adjusted to fit the desired value by multiplying it with a twist factor. $$ \begin{aligned} \delta &= \omega-\sigma \\ T &= \left(e^{i\delta}\right)^2 \\ &= e^{i2\delta} \\ &= \cos(2\delta) + i \, \sin(2\delta) \\ z_\omega^2 &\approx z_\sigma^2 \cdot T \end{aligned} \tag {27} $$ Because the value of $2\delta$ is expected to be extremely small, the Cosine and Sine values can be calculated using the series shown below with very few terms to great accuracy.
Normalizing a Value
A similar approach can be used to normalize a value instead of its square. In this case rather than multiple a value by its reciprocal conjugate, the average is taken. $$ \begin{aligned} z_n &= ( z + (1/z)^* ) / 2 \\ &\approx \left[ (1+d) e^{i\theta} + (1-d) e^{i\theta} \right] / 2 \\ &= e^{i\theta} \\ \end{aligned} \tag {28} $$ When using the Cartesian form it becomes: $$ \begin{aligned} z_n &= ( z + (1/z)^* ) / 2 \\ &= \left[ ( x + i \, y ) + ( x - i \, y )^* / ( x^2 + y^2 ) \right] / 2 \\ &= ( x + i \, y ) \left[ 1 + 1 / ( x^2 + y^2 ) \right] / 2 \\ \end{aligned} \tag {29} $$ The values from the previous section can be used to provide an easy numerical example. $$ \begin{aligned} z_n &= ( z + (1/z)^* ) / 2 \\ &= \left[ (0.61 + i \, 0.81) + [(0.61 - i \, 0.81)/ (0.61^2 + 0.81^2) ]^* \right] / 2 \\ &= \left[ (0.61 + i \, 0.81) + ( 0.593269792 - i \,0.787784478 )^* \right] / 2 \\ &= \left[ 1.203269792 + i \, 1.597784478 \right] / 2 \\ &= 0.601634896 + i \, 0.798892239 \\ \end{aligned} \tag {30} $$ Checking the magnitude. $$ \| 0.601634896 + i \, 0.798892239 \| = 1.000096674 \tag {31} $$ Which is very close to one. Then checking the angles by comparing slopes. $$ \begin{aligned} 0.81 / 0.61 &= 1.327868852 \\ 0.798892239 / 0.601634896 &= 1.327868853 \\ \end{aligned} \tag {32} $$ Which are equal to the extent of the precision used in the division calculation. This is expected to always be the same due to the angle of the reciprocal is always the negative of the angle independent of the magnitudes.
This means the angle is preserved and this operation can be repeated on a complex number that is far from the unit circle to find its unit direction without suffering any angle drift. The limit value can then be stated as: $$ \begin{aligned} z_u &= z / \| z \| \\ \| z_u \| &= 1 \\ \arg( z_u ) &= \arg( z ) \\ \end{aligned} \tag {33} $$
Simplified Normalization
Borrowing the same trick as (19), the normalizing calculation becomes: $$ \begin{aligned} z_n &= ( x + i \, y ) \left[ 3 - ( x^2 + y^2 ) \right] / 2 \\ &= z \left( 3 - zz^* \right) / 2 \\ \end{aligned} \tag {34} $$ The numerical check is even simpler for this. $$ \begin{aligned} z_n &= z \left( 1.9718 \right) / 2 = z \cdot 0.9859 \\ &= (0.61 + i \, 0.81) \cdot 0.9859 = 0.601399 + i \,0.798579 \\ \|z_n \|&= 0.999704545 \\ \end{aligned} \tag {35} $$ Using the trick to replace a division by a multiplication cost a little bit of precision. However, in practice the original value of $z$ is not likely to be 0.01 off in each coordinate, but more like 0.000001 or less, so this precision loss will likely be negligible.
This form also shows that the angle is preserved since the value is simply multiplied by a scalar.
Recursive Tone Generation
In math, with perfect precision, generating a tone from two initial values is very easy without having to know the angles. A simple recursive formula does the trick. $$ z[n] = (z[n-1])^2 / z[n-2] \tag {36} $$ If both $z[n-1]$ and $z[n-2]$ have magnitudes of one, then the next value will also have a magnitude of one, and the generated sequence will stay on the unit circle. It also allows the division to be done by multiplication without having to divide by the magnitude squared, since it is also one. $$ z[n] = (z[n-1])^2 \cdot (z[n-2])^* \tag {37} $$ This formula can be validated by using the previous tone definition equation (7). $$ \begin{aligned} z[n] &= \left( e^{i([n-1] \omega + \phi)} \right)^2 \cdot \left( e^{i([n-2] \omega + \phi)} \right)^* \\ &= \left( e^{i2([n-1] \omega + \phi)} \right) \cdot \left( e^{-i([n-2] \omega + \phi)} \right) \\ &= \left( e^{i(2n \omega - 2 \omega + 2\phi)} \right) \cdot \left( e^{i( -n \omega + 2 \omega - \phi)} \right) \\ &= e^{i(n \omega + \phi)} \\ \end{aligned} \tag {38} $$ Any two values on the unit circle can be used without knowing the angle values and the resulting sequence will be a set of points with equal spacing starting with the two values. To match the tone parameters in the definition, the initial values would be set using the parameters. $$ \begin{aligned} z[0] &= z_0 = x_0 + i\ y_0 = e^{i\phi} \\ z[1] &= z_1 = x_1 + i\ y_1 = e^{i( \omega + \phi) } \\ \end{aligned} \tag {39} $$ If the angles aren't known, the step angle $\omega$ can be recovered by dividing any two values. $$ \begin{aligned} z[n] / z[n-1] &= \frac{ (x[n] + i \, y[n] ) \cdot( x[n-1] + i \, y[n-1] )}{(x[n-1])^2 + (y[n-1])^2 } \\ &= e^{i(n\omega + \phi)} / e^{i([n-1] \omega + \phi)}\\ &= e^{i\omega} \\ \end{aligned} \tag {40} $$ Then take the $\arg$ to find the angle. The $\arg$ function ignores the magnitudes so any radial drift is not a concern in this operation. $$ \omega = \arg \left( e^{i\omega} \right) = \arg \left( z[n] / z[n-1] \right) \tag {41} $$ The value of $\omega$ defines the frequency of the tone. This is the simplest form of a time domain complex tone frequency formula. If the signal has an amplitude multiplier on it, the division operation cancels it out. This formula is sensitive to noise. For formulas that use a larger number of terms to calculate the frequency in the time domain which are then less noise sensitive can be found in my series of "Near Instantaneous ..." articles [4], [5],[6]. Although they were derived for real tone signals, they also work for complex tone signals.
Dealing with the Drift
When this is implemented numerically, the same drifting problem occurs as with the repeated multiplication approach. Instead of doing a straight up squaring, the mitigating squaring technique explained earlier should be used for the first factor, and optionally the value normalizing done on the second.
When building a sequence programmatically, at the next iteration $z[n] = x[n] + i \, y[n]$ becomes $z[n-1]$ that has its square mitigated. A value normalization can then be done at this point with hardly any extra calulations. In the squaring process, the inaccuracy introduce by not having perfect precision can be expressed in two ways that are equal using (14) and (19). $$ \begin{aligned} 1 + g &= ( 1 + d )^ 2 = 1 + 2d + d^2 \\ g &= 2d + d^2 \\ d &\approx g / 2 \\ \end{aligned} \tag {42} $$ This allows the $z[n-1]$ value to be replaced by $(1-g/2) z[n-1]$ in the sequence. Thus when the next iteration occurs and it becomes $z[n-2]$ it has already been normalized.
Finite Sequences Starting at Unity
Generated unit complex tone sequences that have $\phi=0$ always start with $z[0] = 1$. It also means that all the mitigated square values used in the calculation will also be part of the sequence and can be stored at the appropriate location and won't have to be recalculated when using the recursive method.
When its turn comes in the sequence, then the mitigated squaring process should be applied on the precalculated value to fill its slot in the sequence providing it is inside the finite range of the sequence count. Otherwise just skip to the next value which will be odd and need calculating.
In this manner, only odd terms in the sequence will need to be calculated.
Multiplicative Tone Generation
The simpler way to generate a pure unit complex tone is to start with an initial value and the desired twist factor. The next value is then simply the prior value multiplied by the twist factor. $$ \begin{aligned} T &= e^{i\omega} \\ z[0] &= e^{i\phi } \\ z[n] &= z[n-1] \cdot T \end{aligned} \tag {43} $$ In Cartesion form, this would be: $$ \begin{aligned} T &= \cos(\omega) + i \, \sin(\omega) \\ x[0] + i \, y[0] &= \cos(\phi) + i \, \sin(\phi) \\ x[n] + i \, y[n] &= x[n-1] + i \, y[n-1] \cdot T \\ x[n] &= x[n-1] \cos(\omega) - y[n-1] \sin(\omega) \\ y[n] &= x[n-1] \sin(\omega) + y[n-1] \cos(\omega) \\ \end{aligned} \tag {44} $$ This method will accumulate drift. This can be rectified by applying one of the normalization techniques each step, or every few steps depending on application requirements.
Roots of Unity
The Nth Roots of Unity are the solutions to the equation $z=1^{1/N}$. These can be defined as a unit complex tone with $\phi = 0 $ and $\omega = 2\pi/N$. $$ R[n] = e^{i \frac{2\pi}{N} n} \tag {45} $$ From this, it is easy to see the roots are a sequence of powers of the first step which can also be considered a twist factor. $$ R[n] = R[1]^n \tag {46} $$ Confirming that each value is an Nth Root of Unity simply takes raising it to the Nth power. $$ R[n]^N = \left( e^{i \frac{2\pi}{N} n} \right)^N = e^{i 2\pi n} = ( e^{i 2\pi } )^n = 1^n = 1 \tag {47} $$ The last steps can be done since $n$ is an integer. Raising 1 to a rational value will always result in a solution set which falls on a set of Roots of Unity of the denominator. $$ 1^{(a/b)} = \left[ 1^{(1/b)} \right] ^a \tag {48} $$ Since $a$ is an integer, and the sequence is unity based, the raising any $b$th root will land on a $b$th root. On the other hand, if 1 is raised to an irrational number, the set of solutions will cover the unit circle with a countable infinity number of points, none being coincident, far short of the uncountable infinity that are on it. This is part of the subject matter of a branch of Mathematics known as Real Analysis.
Usage in the Discrete Fourier Transform
The definition of a Discrete Fourier Transform (DFT) with a $1/N$ normalization is: $$ \begin{aligned} Z_k &= \frac{1}{N} \sum_{n=0}^{N-1} S_n e^{-i\frac{2\pi}{N}kn} \end{aligned} \tag {49} $$ Where $S_n$ is the $n$th value in the signal frame. In DSP literature, this is often called x[n], and the bin value $Z_k$ is referred to as X[k]. In software libraries, there is generally no normalization factor applied. The advantage of a $1/N$ factor is the DC bin ($Z_0$ or X[0]) is then the average value of the signal living up to the name derived from electrical signals of Direct Current. The other bin values will also retain about the same magnitude independent of the number of samples in the frame which equals the number of bins in the DFT.
The exponent factor in the summation can be understood in terms of the Roots of Unity. If the signal is real valued, the DFT defintion has the same form as a weighted average calculation where the signal values are the weights and the Roots of Unity are being averaged [2]. $$ e^{-i\frac{2\pi}{N}kn} = \left( e^{-i\frac{2\pi}{N}k} \right)^n = R[-k]^n = R[N-k]^n = T_k^n \tag {50} $$ The last step allows the indices to be positive. In Cartesian form: $$ e^{-i\frac{2\pi}{N}kn} = \cos \left( \frac{2\pi}{N}kn \right) + i \, \sin \left( \frac{2\pi}{N}kn \right) \tag {51} $$ If a window function is used, it gets multiplied along with the signal. A normalization factor can be incorporated into the window function values to save a multiplication each step. So, using the more common alternative nomenclature: $$ X[k] = \sum_{n=0}^{N-1} w[n] x[n] e^{-i\frac{2\pi}{N}kn} \tag {52} $$ The notable distinction between the complex tones generated in the previous sections and the exponent term in the DFT definition is the negative sign. This is also a matter of convention which places tones with $k$ cycles per frame into bin $k$ rather than $-k$. For real tone signals, this doesn't matter as they are actually a sum of a $k$ and $-k$ complex tones.
The complex tone can be extended in the negative direction to cover these cases. Since the definition uses a set of Roots of Unity for the exponent term, the periodicity of N can be used to allow for a lookup table with just one period of N values. $$ \begin{aligned} e^{-i\frac{2\pi}{N}p} &= e^{i\frac{2\pi}{N}(-p \mod N)} \\ 0 &<= -p \mod N < N \end{aligned} \tag {53} $$ Where $\mod$ is the integer modulus operator.
If z[n] is a the generated complex tone for the Nth Roots of Unity, the DFT definition can be rewritten as: $$ Z_k = \frac{1}{N} \sum_{n=0}^{N-1} S_n R[-kn \mod N] \tag {54} $$ In Cartesian form, which an implementation will use: $$ Z_k = \frac{1}{N} \sum_{n=0}^{N-1} S_n \left[ \cos\left( \frac{2\pi}{N}(kn \mod N)\right) - i \, \sin\left( \frac{2\pi}{N}(kn \mod N)\right) \right] \tag {55} $$ This simplifies the implementation as shown in the sample program in Appendix I. Each step of $n$ advance the index by $k$ and when the index reaches N or past, N is simply subtracted.
Angle Addition and Subtraction Formulas
When $r_1$ and $r_2$ are both one, (3) and (1) can be used to derive angle addition and subtraction formulas [7]. $$ \begin{aligned} e^{i(\alpha + \beta)} &= e^{i\alpha} \cdot e^{i\beta} \\ \cos(\alpha + \beta) + i \, \sin(\alpha + \beta) &= ( \cos(\alpha) + i \, \sin(\alpha) ) \cdot ( \cos(\beta) + i \, \sin(\beta)) \\ &= ( \cos(\alpha) \cos(\beta) - \sin(\alpha) \sin(\beta)) \\ &+ i \, ( \cos(\alpha) \sin(\beta) + \sin(\alpha) \cos(\beta)) \end{aligned} \tag {56} $$ By matching the real parts and the imaginary parts, the two angle addition formulas pop out. $$ \begin{aligned} \cos(\alpha + \beta) &= \cos(\alpha) \cos(\beta) - \sin(\alpha) \sin(\beta) \\ \sin(\alpha + \beta) &= \cos(\alpha) \sin(\beta) + \sin(\alpha) \cos(\beta) \\ &= \sin(\alpha) \cos(\beta) + \cos(\alpha) \sin(\beta)\\ \end{aligned} \tag {57} $$ In the last line terms are reversed as it makes it easier to remember.
The even and odd natures of the Cosine and Sine functions are needed. These come straight from an understanding of their definitions. $$ \begin{aligned} \cos(-\theta) &= \cos(\theta) \\ \sin(-\theta) &= -\sin(\theta) \\ \end{aligned} \tag {58} $$ To get the subtraction formulas substitute $\beta = -\gamma$. $$ \begin{aligned} \cos(\alpha - \gamma) &= \cos(\alpha) \cos(\gamma) + \sin(\alpha) \sin(\gamma) \\ \sin(\alpha - \gamma) &= \sin(\alpha) \cos(\gamma) - \cos(\alpha) \sin(\gamma) \\ \end{aligned} \tag {59} $$ Double angle formulas can also be derived by setting $\beta = \alpha$ in the addition formulas. $$ \begin{aligned} \cos(2\alpha) &= \cos^2(\alpha) - \sin^2(\alpha) \\ &= 2 \cos^2(\alpha) - 1\\ \sin(2\alpha) &= \sin(\alpha) \cos(\alpha) + \cos(\alpha) \sin(\alpha)\\ &= 2 \sin(\alpha) \cos(\alpha) \\ \end{aligned} \tag {60} $$ The second line in the Cosine formula comes from the fundamental relation $\cos^2(\alpha) + \sin^2(\alpha) = 1$.
Radians and Degrees
So far, all the angles have been measured using radians as units. Like with Logarithms, there is a mathematical scale (Radians) and a commonly used scale (Degrees). These are proportional with a conversion factor. Generally the distinction is more important in Calculus. The angle addition and subtraction formulas work in either scale as long as the units are consistent. In this article the symbols $\sin_D$ and $\cos_D$ will be used for the degree versions. $$ \begin{aligned} \cos_D(A) &= \cos \left( \frac{\pi}{180} A \right) \\ \sin_D(A) &= \sin \left( \frac{\pi}{180} A \right) \\ \end{aligned} \tag {61} $$
The Series that are Needed
In order to build a degree based trigonometry table, the value of $\pi$ will be needed and series solutions for Sine and Cosine values. These are done in radians. $$ \arctan(m) = m - \frac{m^3}{3} + \frac{m^5}{5} - \frac{m^7}{7} + \frac{m^9}{9} + ... \tag {62} $$ Here $m$ is the slope of the angle, calculated as $y/x$ for a Cartesian point, or $\sin(\theta)/\cos(\theta)=\tan(\theta)$. $$ \cos(\theta) = 1 - \frac{\theta^2}{2} + \frac{\theta^4}{4!} - \frac{\theta^6}{6!} + ... \tag {63} $$ $$ \sin(\theta) = \theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \frac{\theta^7}{7!} + ... \tag {64} $$ The even and odd nature of Cosine and Sine are revealed in these series formulas.
Necessary Ingredients
First the value of $\pi$ will be needed for the conversion factor from degrees to radians. In order to calculate it by hand, a special formula is used along with the series for ArcTangent.
The simplest series solution is: $$ \pi = 4 ( \arctan( 1 ) ) \tag {65} $$ The problem with this is that the series converges extremely slowly. A fancier formula that can be used is: $$ \begin{aligned} \pi &= 4 ( 4 \arctan( 1/5 ) - \arctan( 1/239 ) ) \\ &\approx 4 ( 4 \cdot 0.1973955 - 0.0041841 ) \\ &= 3.14159265 \\ \end{aligned} \tag {66} $$ This is close enough for demonstration purposes. By using extended precision as many digits as desired can be calculated. This can then be used to calculate the degrees to radians conversion factors. $$ \frac{\pi}{180} = \frac{3.14159265}{180} = 0.017453292 \tag {67} $$ This is also, by default, the value of one degree in radians.
Nature of the Table
In order to build a complete trig table, only the degree values of zero to 45 need to be calculated. The rest can be completed using some simple rules. The values from 46 to 90 can be obtained by reversing the values of Sine and Cosine and reversing the order. Then the values from 91 to 180 and be obtained using the values from 89 to 0, and changing the sign of the Cosine to negative. Then the remaining 181 to 359 flip the signs from the 1 to 179 values.
The equations for 46 to 90 are: $$ \begin{aligned} \cos_D( A ) &= \sin_D( 90 - A ) \\ \sin_D( A ) &= \cos_D( 90 - A ) \\ \end{aligned} \tag {68} $$ The equations for 91 to 180: $$ \begin{aligned} \cos_D( A ) &= -\cos_D( 180 - A ) \\ \sin_D( A ) &= \sin_D( 180 - A ) \\ \end{aligned} \tag {69} $$ The equations for 181 to 359: $$ \begin{aligned} \cos_D( A ) &= -\cos_D( A - 180 ) \\ \sin_D( A ) &= -\sin_D( A - 180 ) \\ \end{aligned} \tag {70} $$ The important point is that only the first 45 values need to be calculated. The value for zero is set by unity.
Hybrid Approach When Done by Hand
To generate a trigonometry table, several approaches have been described in this article so far. The simplest is to start with the values of $z[1] = \cos_D(1) + i \, \sin_D(1)$, and then perform repeat multiplications.
A slightly more sophisticated approach would be to start with the one degree value, then use the double angle formulas to build values for 2, 4, 8, 16, 32. Then the angle addition formulas can be used corresponding to binary digits to construct the rest of the values. Since subtraction formulas can also be used, for example 15 can be calculated from 16-1 rather than 1+2+4+8 saving two iterations and some rounding errors.
This approach can be made more accurate by using the mitigated squaring technique instead of the double angle formulas. Undoubtedly, because of the small range needed, and having the initial values calculated to an extra digit or two, will be adequate.
Lastly, recognizing it as a complex tone sequence, either the recursive or multiplicative methods described earlier can also be used.
Vicanek's Double Half Step
An alternative way of generating the complex tone sequence was derived by Martin Vicanek [8] and featured on this site by Rick Lyons [9]. This generator is very efficient and sufficiently accurate for this task if the initial values are accurate enough. Any inaccuracy translates into an elliptical shape rather than circular. This appears as drifting but actually is a bounded variation when compared to a circle. This bound can be made arbitrarily small.
The Significance of Digits
How many significant digits should a calculation use? That depends largely on the calculations and the application. With the advent of calculators and computers, this isn't as much an issue as it used to be when calculations were done with paper and pencil assisted by slide rules. What has been lost in the process is a feel for how much accuracy and precision is required in the process of calculation, but perhaps more importantly, how many are valid in the final result. An understanding of this seems to have been totally lost in the "six foot rule" applied in many place during the recent pandemic.
To get an idea, it is useful to consider a real world example. Consider a mile at 5280 feet. One digit determines a tenth of a mile, two digits about 53 feet, three digits 5 feet, four digits half a foot. For most application that is close enough, so traditionally tables were to 4 digits or five to avoid rounding errors. Five digits is just over half an inch, and six less than a sixteenth of an inch.
For the metric minded, the example is clearer and quicker. Six significant digits on a kilometer takes you down to a millimeter. As a rule of thumb, a single precision floating point in programming is about seven significant digits and double precision is up to 17.
Understanding the effects of rounding errors in a process is then key to understanding how precise the answer is.
Initial Values
All the procedures detailed above rely on having the values of Sine and Cosine for one degree. These values can be calculated directly using the series, but since the easiest way to normalize a value is by the squaring method, and the series will converge a little quicker, and Vicanek's method needs it, figuring the value for a half degree is worth it.
From the conversion value calculation (67), the value for half a degree is found by dividing it by two. $$ \theta_{1/2} = 0.017453292 / 2 = 0.008726646 \tag {71} $$ This is in the neighborhood of 1/100th. So each term in the series calculations, which have a step of $x^2$ which means more than four signficant digits are gained for each additional terms. Of course, the prior calculations also have to be carried out to the same precision.
The Cosine calculation to four terms is then $$ \begin{aligned} \cos(x) &= 1 - \frac{x^2}{2} + \frac{x^4}{4!} - \frac{x^6}{6!} + ... \\ &\approx 1 - \frac{0.008726646^2}{2} + \frac{0.008726646^4}{24} - \frac{0.008726646^6}{720} \\ &= ( 1 + 0.000000000242 ) - ( 0.000038077 + 6.1 \cdot 10^{-16} ) \\ &= 1.000000000242 - 0.000038077 = 0.999961923 \\ \end{aligned} \tag {72} $$ This is the same value a calculator returns. When doing it by hand, usually adding all the negative terms together and then doing one subtraction is easier. The fourth term was unnecessary.
The Sine calculation to three terms is then $$ \begin{aligned} \sin(x) &= x - \frac{x^3}{3!} + \frac{x^5}{5!} \\ &\approx 0.008726646 - \frac{0.008726646^3}{6} + \frac{0.008726646^5}{120}\\ &= ( 0.008726646 + 4.2 \cdot 10^{-13} ) - ( 0.000000111 ) \\ &= 0.008726535 \\ \end{aligned} \tag {73} $$ Here even the third term wasn't needed. Again, this result matches the calculator. When doing this by hand, the value of something like two degrees would also be worth calculating as it would save building up from a half while still converging very quickly.
The rest of the calculations are shown in the program example in Appendix I.
Linear Interpolation
For most cases of interpolation, simple linear interpolation will be adequate. When the angle falls between two values on the table then its value is partitioned an integer portion ($A$) and fraction ($f$) where the fraction is between zero and one. $$ \begin{aligned} \theta &= A + f \\ B &= A + 1 \\ \end{aligned} \tag {74} $$ Once A and B are known, look up the Sine and Cosine on your table and do the following calculations, depending on which you need. $$ \begin{aligned} \cos_D( \theta ) &= ( 1 - f ) \cos_D( A ) + f \cos_D( B ) \\ \sin_D( \theta ) &= ( 1 - f ) \sin_D( A ) + f \sin_D( B ) \\ \end{aligned} \tag {75} $$ The error in magnitude will be greatest in the middle, but the angle will be most correct there. With one degree spacing, in most cases it really doesn't matter, but the value can be normalized using the simple method for results that will be as accurate as the table values.
Taylor Series Interpolation
Should more accuracy be required, and your table is precise enough, a Taylor series expansion around the nearest value is simple to construct due to the way the trig functions derivatives work. First the nearest degree value is found with an offset fraction having a value between negative one half and positive one half. $$ \theta = A + f \tag {76} $$ Without delving into too much detail, the formulas become: $$ \begin{aligned} \cos_D( \theta ) &= \cos_D( A ) - f \sin_D( A ) - \frac{f^2}{2} \cos_D( A ) \\ \sin_D( \theta ) &= \sin_D( A ) + f \cos_D( A ) - \frac{f^2}{2} \sin_D( A ) \\ \end{aligned} \tag {77} $$ Of course, a third order term could be added, but it won't improve results any except in very precise tables.
Inverse Interpolation
Just like in my common logarithms article [3], a forward table can be used to make inverse evaluations. If both the Cosine and Sine are known, the next steps can be skipped. Otherwise, the one you don't have needs to be calculated from the one you do. The procedure is the same for either value. It requires taking a square root, so Newton's method will be used.
The equation to be solved is: $$ u_0 = \sqrt{ 1 - v^2 } \approx 1 - v^2/2 \tag {78} $$ From there, the solution will quickly converge using this equation: $$ u_{n+1} = ( 1 + u_n^2) / ( 2u_n ) \tag {79} $$ Where $v$ is the known value and $u$ is the unknown value.
Once both the Cosine and Sine are known, find the closest matching row in the table. The next step is to use the angle subtraction formulas to subtract the row's angle ($A$) to find the Cosine and Sine of the fractional part of the angle ($f$).
From there it simply a matter of calculating the tangent, finding the angle in radians $\delta$ using the $\arctan$ series, then converting it to degrees. $$ \begin{aligned} \tan_D(f) &= \sin_D(f) / \cos_D(f) \\ \delta &= \arctan( \tan_D(f) ) \\ &\approx \tan_D(f) - (\tan_D(f))^3/3 + (\tan_D(f))^5/5 \\ f &= \delta\cdot \frac{180}{\pi} \\ \theta_D &= A + f \\ \end{aligned} \tag {80} $$ The third term is likely not to be needed for any practical application. The tangent of an angle is still the slope of the line no matter whether the degree or radian version is being used. Same for cosine and sine.
Programming Examples
The two appendices contain source code and their output that illustrate the many of concepts and equations presented in this article, but not all. Feel free to copy and modify the code without attribution except for the tone parameter estimation routines or Vicanek's method, those deserve a shout out, the rest is standard stuff. It is written in straightforward C other than the double slash comments which most C compilers handle. Only <stdio.h> is needed, the <math.h> file is not. That is the point of this article, after all. The output was produced on a Raspberry Pi4 running Ubuntu.
Conclusion
A smorgasbord of calculation techniques related to the basic trigonometric functions have been presented in order to give the reader, particularly the novice, an idea of varioius trade-offs when solving a problem. The title goal of producing an actual trig table is done in the appendices and not done by hand, though clearly they could be. Understanding these functions is key to understanding how the DFT works.
Appendix I - Programming Example
Here is a programming example illustrating some of the methods described in this article as well as my frequency calculation formula [10] and phase and amplitude calculation [11].
#include <stdio.h>
double Cos[360];
double Sin[360];
double DegreesPerRadian;
double RadiansPerDegree;
//============================================================================
//============================================================================
// T R I G O N O M E T R Y
//============================================================================
//============================================================================
double SmallSlopeArcTan( double m, double Cutoff )
{
double a = 0.0;
double t = m;
double mm = m * m;
for( int p = 1; p < 20; p += 2 )
{
if( ( p & 2 ) == 0 ) { a += t / (double) p; }
else { a -= t / (double) p; }
if( t < Cutoff ) break;
t *= mm;
}
return a;
}
//============================================================================
double SmallAngleCos( double a, double Cutoff )
{
double c = 1.0; // Cosine
double t = 1.0; // Power term
double aa = a * a;
for( int p = 2; p < 20; p += 2 )
{
t *= aa / (double) ( p * ( p - 1 ) );
if( t < Cutoff ) break;
if( ( p & 2 ) == 0 ) { c += t; }
else { c -= t; }
}
return c;
}
//============================================================================
double SmallAngleSin( double a, double Cutoff )
{
double s = a; // Sine
double t = a; // Power term
double aa = a * a;
for( int p = 3; p < 20; p += 2 )
{
t *= aa / (double) ( p * ( p - 1 ) );
if( t < Cutoff ) break;
if( ( p & 2 ) == 0 ) { s += t; }
else { s -= t; }
}
return s;
}
//============================================================================
void NextDegreeTrigValues( int d, double cd, double sd )
{
double tx, ty, t2, nf;
tx = cd * Cos[d-1] - sd * Sin[d-1];
ty = sd * Cos[d-1] + cd * Sin[d-1];
t2 = ( tx * tx + ty * ty );
nf = ( 3.0 - t2 ) * 0.5;
Cos[d] = tx * nf;
Sin[d] = ty * nf;
}
//============================================================================
void LinearValues( double angle, double* re, double* im )
{
while( angle >= 360.0 ) angle -= 360;
while( angle < 0.0 ) angle += 360;
int A = (int) angle;
double f = angle - A;
if( f < 0.0 )
{
if( A == 0 ) { A = 359; }
else { A -= 1; }
f += 1.0;
}
int B = A + 1; if( B == 360 ) B = 0;
double g = 1.0 - f;
*re = g * Cos[A] + f * Cos[B];
*im = g * Sin[A] + f * Sin[B];
}
//============================================================================
//============================================================================
// D I S C R E T E F O U R I E R T R A N S F O R M
//============================================================================
//============================================================================
void BinValueCalculation( double Sig_re[], double Sig_im[], int k,
double Bin_re[], double Bin_im[] )
{
int d = 0;
double sum_re = 0.0;
double sum_im = 0.0;
for( int n = 0; n < 360; n++ )
{
sum_re += Sig_re[n] * Cos[d] - Sig_im[n] * -Sin[d];
sum_im += Sig_re[n] * -Sin[d] + Sig_im[n] * Cos[d];
d += k; if( d >= 360 ) d -= 360;
}
Bin_re[k] = sum_re / 360.0;
Bin_im[k] = sum_im / 360.0;
}
//============================================================================
//============================================================================
// C O M P L E X O P E R A T I O N S
//============================================================================
//============================================================================
void ComplexMultiply( double re1, double im1,
double re2, double im2,
double* reP, double* imP )
{
*reP = re1 * re2 - im1 * im2;
*imP = im1 * re2 + re1 * im2;
}
//============================================================================
void ComplexDivide( double reNum, double imNum,
double reDen, double imDen,
double* reQ, double* imQ )
{
double den2 = reDen * reDen + imDen * imDen;
*reQ = ( reNum * reDen + imNum * imDen ) / den2;
*imQ = ( imNum * reDen - reNum * imDen ) / den2;
}
//============================================================================
//============================================================================
// D A W G ' S P A R A M E T E R F O R M U L A S
//============================================================================
//============================================================================
void ComplexToneFrequencyFormula( int k, double* omega,
double Bin_re[], double Bin_im[] )
{
int M = k - 1; if( M < 0 ) M = 359;
int P = k + 1; if( P > 359 ) P = 0;
double G_reM = Bin_re[M]; double G_imM = Bin_im[M];
double G_re0 = Bin_re[k]; double G_im0 = Bin_im[k];
double G_reP = Bin_re[P]; double G_imP = Bin_im[P];
double Ave_re = ( G_reM + G_re0 + G_reP ) / 3.0;
double Ave_im = ( G_imM + G_im0 + G_imP ) / 3.0;
double K_reM = G_reM - Ave_re; double K_imM = G_imM - Ave_im;
double K_re0 = G_re0 - Ave_re; double K_im0 = G_im0 - Ave_im;
double K_reP = G_reP - Ave_re; double K_imP = G_imP - Ave_im;
double n_re = K_reM * Bin_re[M] - K_imM * Bin_im[M]
+ K_re0 * Bin_re[k] - K_im0 * Bin_im[k]
+ K_reP * Bin_re[P] - K_imP * Bin_im[P];
double n_im = K_imM * Bin_re[M] + K_reM * Bin_im[M]
+ K_im0 * Bin_re[k] + K_re0 * Bin_im[k]
+ K_imP * Bin_re[P] + K_reP * Bin_im[P];
double d_reM, d_imM;
double d_re0, d_im0;
double d_reP, d_imP;
ComplexMultiply( K_reM, K_imM, Bin_re[M], Bin_im[M], &d_reM, &d_imM );
ComplexMultiply( K_re0, K_im0, Bin_re[k], Bin_im[k], &d_re0, &d_im0 );
ComplexMultiply( K_reP, K_imP, Bin_re[P], Bin_im[P], &d_reP, &d_imP );
ComplexMultiply( Cos[1], Sin[1], d_reM, d_imM, &d_reM, &d_imM );
ComplexMultiply( Cos[1], -Sin[1], d_reP, d_imP, &d_reP, &d_imP );
double d_re = d_reM + d_re0 + d_reP;
double d_im = d_imM + d_im0 + d_imP;
double a_re, a_im;
ComplexDivide( n_re, n_im, d_re, d_im, &a_re, &a_im );
double m = a_im / a_re;
double d = SmallSlopeArcTan( m, 0.00000000000001 );
*omega = k + d * DegreesPerRadian;
}
//============================================================================
void ComplexToneVirtualBin( int Center_k, double DegreesPerStep,
double Bin_re[], double Bin_im[],
double* Virtual_re, double* Virtual_im )
{
//--- Allocate local values
double Y_re[360];
double Y_im[360];
double twist_re, twist_im;
double sin_Nhalf, cos_Nhalf;
//--- Calculate Y
double Nminus1over2 = 359.0 / 2.0;
double alpha = DegreesPerStep * RadiansPerDegree;
for( int k = Center_k - 1; k <= Center_k + 1; k++ )
{
double beta_k = k * RadiansPerDegree;
double delta_k = beta_k - alpha;
double angle_k = ( -delta_k * Nminus1over2 );
LinearValues( angle_k * DegreesPerRadian, &twist_re, &twist_im );
double half = delta_k / 2.0;
double sin_half = SmallAngleSin( half, 0.00000000000001 );
LinearValues( 360 * half * DegreesPerRadian,
&cos_Nhalf, &sin_Nhalf );
double heft = sin_Nhalf / ( 360 * sin_half );
Y_re[k] = twist_re * heft;
Y_im[k] = twist_im * heft;
}
//--- Do dot products
double ZYstar_re = 0.0;
double ZYstar_im = 0.0;
double YYstar_re = 0.0;
// double YYstar_im = 0.0;
for( int k = Center_k - 1; k <= Center_k + 1; k++ )
{
ZYstar_re += Bin_re[k] * Y_re[k] + Bin_im[k] * Y_im[k];
ZYstar_im += Bin_im[k] * Y_re[k] - Bin_re[k] * Y_im[k];
YYstar_re += Y_re[k] * Y_re[k] + Y_im[k] * Y_im[k];
// YYstar_im += Y_im[k] * Y_re[k] - Y_re[k] * Y_im[k];
// Zero by formula
}
//--- Calculate;
*Virtual_re = ZYstar_re / YYstar_re;
*Virtual_im = ZYstar_im / YYstar_re;
/*
ComplexDivide( ZYstar_re, ZYstar_im,
YYstar_re, YYstar_im,
Virtual_re, Virtual_im );
*/
}
//============================================================================
//============================================================================
// M A I N R O U T I N E
//============================================================================
//============================================================================
int main( int ArgCount, char* ArgValues[] )
{
//--- Set precision
double Cutoff = 0.00000000000000001;
//--- Calculate Pi
double m1 = 1.0 / 5.0;
double a1 = SmallSlopeArcTan( m1, Cutoff );
double m2 = 1.0 / 239.0;
double a2 = SmallSlopeArcTan( m2, Cutoff );
double pi = 4.0 * ( 4.0 * a1 - a2 );
//--- Degree to and from Radians Conversion factors
RadiansPerDegree = pi / 180.0;
DegreesPerRadian = 1.0 / RadiansPerDegree;
//--- One degree values
double cd = SmallAngleCos( RadiansPerDegree, Cutoff );
double sd = SmallAngleSin( RadiansPerDegree, Cutoff );
//--- One half degree values
double HalfDegree = RadiansPerDegree * 0.5;
double ch = SmallAngleCos( HalfDegree, Cutoff );
double sh = SmallAngleSin( HalfDegree, Cutoff );
//--- Print results
printf( "\n"
"Fundamental Values\n"
"\n"
" Pi: %17.14f \n"
" Radians Per Degree: %17.14f \n"
" Degrees Per Radian: %17.14f \n"
" Cosine( 1 Degree ): %17.14f \n"
" Sine( 1 Degree ): %17.14f \n"
" Cosine( 1/2 Degree ): %17.14f \n"
" Sine( 1/2 Degree ): %17.14f \n",
pi,
RadiansPerDegree, DegreesPerRadian,
cd, sd, ch, sh );
//--- Build first octant
Cos[0] = 1.0;
Sin[0] = 0.0;
for( int d = 1; d <= 45; d++ )
{
NextDegreeTrigValues( d, cd, sd );
}
//--- Fill second octant
for( int d = 46; d <= 90; d++ )
{
Cos[d] = Sin[90-d];
Sin[d] = Cos[90-d];
}
//--- Fill second quadrant
for( int d = 91; d <= 180; d++ )
{
Cos[d] = -Cos[180-d];
Sin[d] = Sin[180-d];
}
//--- Fill upper half
for( int d = 181; d <= 359; d++ )
{
Cos[d] = -Cos[d-180];
Sin[d] = -Sin[d-180];
}
//--- Print start of table
printf( "\n"
"First few table values\n"
"\n"
" d Cosine Sine\n"
" --- --------- ---------\n" );
for( int d = 0; d < 10; d++ )
{
printf( " %3d %11.7f %11.7f\n",
d, Cos[d], Sin[d] );
}
//--- Build a signal
double Sig_re[360];
double Sig_im[360];
double amp = 19.0;
double phase = 61.0;
double omega = 10.4;
double angle = phase;
for( int n = 0; n < 360; n++ )
{
LinearValues( angle, &Sig_re[n], &Sig_im[n] );
Sig_re[n] *= amp;
Sig_im[n] *= amp;
angle += omega;
if( angle >= 360 ) angle -= 360;
}
//--- Print start of signal
printf( "\n"
"First few signal values\n"
"\n"
" n Real Imaginary\n"
" --- --------- ---------\n" );
for( int n = 0; n < 10; n++ )
{
printf( " %3d %11.7f %11.7f\n",
n, Sig_re[n], Sig_im[n] );
}
//--- Calculate some bins
printf( "\n"
"Selected DFT bin values\n"
"\n"
" k Magnitude^2 Real Imaginary\n"
" --- ----------- ------------ -----------\n" );
double Bin_re[360];
double Bin_im[360];
for( int k = 5; k < 15; k++ )
{
BinValueCalculation( Sig_re, Sig_im, k, Bin_re, Bin_im );
printf( " %3d %11.5f %12.5f %12.5f\n", k,
Bin_re[k] * Bin_re[k] + Bin_im[k] * Bin_im[k],
Bin_re[k], Bin_im[k] );
}
//--- Calculate the parameter values from the bin values
double DegreesPerStep, Virtual_re, Virtual_im;
ComplexToneFrequencyFormula( 10, &DegreesPerStep,
Bin_re, Bin_im );
ComplexToneVirtualBin( 10, DegreesPerStep,
Bin_re, Bin_im,
&Virtual_re, &Virtual_im );
printf( "\n"
"Parameter Calculations\n"
"\n"
" Cycles per frame: %12.5f \n"
" Virtual Bin: %12.5f %12.5f\n",
DegreesPerStep, Virtual_re, Virtual_im );
}
//============================================================================
[Edit 10-05-2022] The source code was upgraded from my original posting.
Fundamental Values
Pi: 3.14159265358979
Radians Per Degree: 0.01745329251994
Degrees Per Radian: 57.29577951308234
Cosine( 1 Degree ): 0.99984769515639
Sine( 1 Degree ): 0.01745240643728
Cosine( 1/2 Degree ): 0.99996192306417
Sine( 1/2 Degree ): 0.00872653549837
First few table values
d Cosine Sine
--- --------- ---------
0 1.0000000 0.0000000
1 0.9998477 0.0174524
2 0.9993908 0.0348995
3 0.9986295 0.0523360
4 0.9975641 0.0697565
5 0.9961947 0.0871557
6 0.9945219 0.1045285
7 0.9925462 0.1218693
8 0.9902681 0.1391731
9 0.9876883 0.1564345
First few signal values
n Real Imaginary
--- --------- ---------
0 9.2113828 16.6177744
1 6.0600061 18.0069413
2 2.7098821 18.8052903
3 -0.7293490 18.9855328
4 -4.1445709 18.5417405
5 -7.4238914 17.4895922
6 -10.4587511 15.8615300
7 -13.1504010 13.7131090
8 -15.4098361 11.1139263
9 -17.1627435 8.1494689
Selected DFT bin values
k Magnitude^2 Real Imaginary
--- ----------- ------------ -----------
5 1.13536 -0.68917 0.81265
6 1.70965 -0.85437 0.98980
7 2.86265 -1.11668 1.27109
8 5.74434 -1.59750 1.78671
9 16.87969 -2.76506 3.03877
10 206.76676 -9.76995 10.55059
11 91.89680 6.57445 -6.97663
12 12.92371 2.48823 -2.59470
13 4.89471 1.54517 -1.58340
14 2.55351 1.12598 -1.13387
Parameter Calculations
Cycles per frame: 10.40000
Virtual Bin: 9.21143 16.61786
The virtual bin is as if the signal had a whole number of cycles in the frame and the bin value represents the amplitude and phase directly.
Rather than add code to do these extra operations, a calculator was used instead, and the truncated display values as well. $$ \begin{aligned} \text{Amplitude} = \sqrt{9.21143^2 + 16.61786^2} &= 19.000097727 \\ \text{Phase} = \arctan_D{16.61786/9.21143} &= 61.000000564 \\ \end{aligned} \tag {81} $$
Appendix II - Vicanek vs Dawg
This is a comparison between Vicanek's oscillator and the technique described in this article.
#include <stdio.h>
//============================================================================
void DoubleHalfStep( double x, double y, double* rx, double* ry,
double k1, double k2 )
{
double w;
w = x - k1 * y;
*ry = y + k2 * w;
*rx = w - k1 * (*ry);
}
//============================================================================
void NormalizedTwist( double x, double y, double* rx, double* ry,
double a, double b )
{
double tx, ty, t2, nf;
tx = a * x - b * y;
ty = b * x + a * y;
t2 = ( tx * tx + ty * ty );
nf = ( 3.0 - t2 ) * 0.5;
*rx = tx * nf;
*ry = ty * nf;
}
//============================================================================
int main( int ArgCount, char* ArgValues[] )
{
//--- Allocate common variables
double x, y, rx, ry;
//--- Set known initial values from Appendix I output
double cd = 0.99984769515639;
double sd = 0.01745240643728;
double ch = 0.99996192306417;
double sh = 0.00872653549837;
//--- Run Vicanek Double Half Step
double k1 = sh / ch; // Tangent( theta/2 )
double k2 = sd; // Sine( theta )
x = 1.0;
y = 0.0;
printf( "\n"
"Vicanek's Quadrature Oscillator\n"
"\n"
" n Real Imaginary Magnitude^2 - 1\n"
" --- --------- --------- ---------------\n" );
for( int n = 0; n <= 45; n++ )
{
printf( " %2d %11.7f %11.7f %9.1e\n",
n, x, y, x*x + y*y - 1.0 );
DoubleHalfStep( x, y, &rx, &ry, k1, k2 );
x = rx;
y = ry;
}
printf( "\n" );
//--- Run Dawg Normalized Twist
printf( "\n"
"Dawg's Complex Tone Generator\n"
"\n"
" n Real Imaginary Magnitude^2 - 1\n"
" --- --------- --------- ---------------\n" );
double a = cd; // Cosine( theta )
double b = sd; // Sine( theta )
x = 1.0;
y = 0.0;
for( int n = 0; n <= 45; n++ )
{
printf( " %2d %11.7f %11.7f %9.1e\n",
n, x, y, x*x + y*y - 1.0 );
NormalizedTwist( x, y, &rx, &ry, a, b );
x = rx;
y = ry;
}
printf( "\n" );
}
//============================================================================
Here is the output, which contains the values needed to fill a trig table.
Vicanek's Quadrature Oscillator
n Real Imaginary Magnitude^2 - 1
--- --------- --------- ---------------
0 1.0000000 0.0000000 0.0e+00
1 0.9998477 0.0174524 2.2e-16
2 0.9993908 0.0348995 4.4e-16
3 0.9986295 0.0523360 6.7e-16
4 0.9975641 0.0697565 1.3e-15
5 0.9961947 0.0871557 1.8e-15
6 0.9945219 0.1045285 2.7e-15
7 0.9925462 0.1218693 3.6e-15
8 0.9902681 0.1391731 4.4e-15
9 0.9876883 0.1564345 5.8e-15
10 0.9848078 0.1736482 6.9e-15
11 0.9816272 0.1908090 8.4e-15
12 0.9781476 0.2079117 1.0e-14
13 0.9743701 0.2249511 1.2e-14
14 0.9702957 0.2419219 1.4e-14
15 0.9659258 0.2588190 1.6e-14
16 0.9612617 0.2756374 1.8e-14
17 0.9563048 0.2923717 2.1e-14
18 0.9510565 0.3090170 2.3e-14
19 0.9455186 0.3255682 2.6e-14
20 0.9396926 0.3420201 2.9e-14
21 0.9335804 0.3583679 3.2e-14
22 0.9271839 0.3746066 3.5e-14
23 0.9205049 0.3907311 3.8e-14
24 0.9135455 0.4067366 4.1e-14
25 0.9063078 0.4226183 4.4e-14
26 0.8987940 0.4383711 4.8e-14
27 0.8910065 0.4539905 5.1e-14
28 0.8829476 0.4694716 5.5e-14
29 0.8746197 0.4848096 5.9e-14
30 0.8660254 0.5000000 6.2e-14
31 0.8571673 0.5150381 6.6e-14
32 0.8480481 0.5299193 7.0e-14
33 0.8386706 0.5446390 7.4e-14
34 0.8290376 0.5591929 7.8e-14
35 0.8191520 0.5735764 8.2e-14
36 0.8090170 0.5877853 8.6e-14
37 0.7986355 0.6018150 9.0e-14
38 0.7880108 0.6156615 9.4e-14
39 0.7771460 0.6293204 9.8e-14
40 0.7660444 0.6427876 1.0e-13
41 0.7547096 0.6560590 1.1e-13
42 0.7431448 0.6691306 1.1e-13
43 0.7313537 0.6819984 1.2e-13
44 0.7193398 0.6946584 1.2e-13
45 0.7071068 0.7071068 1.2e-13
Dawg's Complex Tone Generator
n Real Imaginary Magnitude^2 - 1
--- --------- --------- ---------------
0 1.0000000 0.0000000 0.0e+00
1 0.9998477 0.0174524 2.2e-16
2 0.9993908 0.0348995 2.2e-16
3 0.9986295 0.0523360 -1.1e-16
4 0.9975641 0.0697565 2.2e-16
5 0.9961947 0.0871557 -1.1e-16
6 0.9945219 0.1045285 0.0e+00
7 0.9925462 0.1218693 2.2e-16
8 0.9902681 0.1391731 2.2e-16
9 0.9876883 0.1564345 0.0e+00
10 0.9848078 0.1736482 2.2e-16
11 0.9816272 0.1908090 0.0e+00
12 0.9781476 0.2079117 2.2e-16
13 0.9743701 0.2249511 0.0e+00
14 0.9702957 0.2419219 2.2e-16
15 0.9659258 0.2588190 0.0e+00
16 0.9612617 0.2756374 4.4e-16
17 0.9563048 0.2923717 2.2e-16
18 0.9510565 0.3090170 0.0e+00
19 0.9455186 0.3255682 0.0e+00
20 0.9396926 0.3420201 0.0e+00
21 0.9335804 0.3583679 0.0e+00
22 0.9271839 0.3746066 2.2e-16
23 0.9205049 0.3907311 0.0e+00
24 0.9135455 0.4067366 2.2e-16
25 0.9063078 0.4226183 -1.1e-16
26 0.8987940 0.4383711 0.0e+00
27 0.8910065 0.4539905 2.2e-16
28 0.8829476 0.4694716 0.0e+00
29 0.8746197 0.4848096 -2.2e-16
30 0.8660254 0.5000000 0.0e+00
31 0.8571673 0.5150381 0.0e+00
32 0.8480481 0.5299193 2.2e-16
33 0.8386706 0.5446390 2.2e-16
34 0.8290376 0.5591929 0.0e+00
35 0.8191520 0.5735764 2.2e-16
36 0.8090170 0.5877853 0.0e+00
37 0.7986355 0.6018150 0.0e+00
38 0.7880108 0.6156615 0.0e+00
39 0.7771460 0.6293204 2.2e-16
40 0.7660444 0.6427876 0.0e+00
41 0.7547096 0.6560590 2.2e-16
42 0.7431448 0.6691306 2.2e-16
43 0.7313537 0.6819984 -1.1e-16
44 0.7193398 0.6946584 0.0e+00
45 0.7071068 0.7071068 0.0e+00
Notice that Vicanek's drifts just a little bit due to the truncated initial values. This may look severe, but if you make a kilometer a millimeter, it is still less than a millimeter there, aka $10^{-12}$.
References
[1] Cedron Dawg, "The Exponential Nature of the Complex Unit Circle", 2015
[2] Cedron Dawg, "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity", 2015
[3] Cedron Dawg, "A Recipe for a Common Logarithm Table", 2017
[4] Cedron Dawg, "Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1)", 2017
[5] Cedron Dawg, "Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 2)", 2017
[6] Cedron Dawg, "Exact Near Instantaneous Frequency Formulas Best at Zero Crossings", 2017
[7] Cedron Dawg, "Angle Addition Formulas from Euler's Formula", 2019
[8] Martin Vicanek, "A New Recursive Quadrature Oscillator", 2015
[9] Rick Lyons, "A New Contender in the Quadrature Oscillator Race", 2022
[10] Cedron Dawg, "Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT", 2017
[11] Cedron Dawg, "Phase and Amplitude Calculation for a Pure Complex Tone in a DFT using Multiple Bins", 2018
- Comments
- Write a Comment Select to add a comment
To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers:
