The Formalization of Discrete Fourier Transform in HOL

Abstract

Traditionally, Discrete Fourier Transform (DFT) is performed with numerical or symbolic computation, which cannot guarantee 100% accurate analysis which may be necessary for safety-critical applications. Machine theorem proving is one of the formal methods that perform accurate analysis with completeness to some extent. This paper proposes the formalization of DFT in a higher-order logic theorem prover named HOL. We propose the formal definition of DFT and verify the fundamental properties of DFT. Two case studies are presented to illustrate usefulness and correctness of the formalized DFT, including formal verifications of Fast Fourier Transform (FFT) and cosine frequency shift.

1. Introduction

Fourier Transform (FT), through which signal is transformed from time domain to frequency one, is a fundamental and core method for signal processing. DFT is the discrete version of FT on both time domain and frequency domain which is adapted to the actual environment of computers and digital signal processing (DSP) chips. DFT has been widely used in many fields such as spectral analysis, data compressing, digital communication, signal processing, and the solution of partial differential equations [1, 2]. Therefore tremendous algorithms and packages have been implemented to perform DFT. The algorithms and packages are based on numerical or symbolic computation. The numerical computation intuitively produces approximate solutions because floating-point numbers are used to represent reals in computers. Computer algebra systems (CAS) such as MATLAB, Maple, and Mathematica offer a large collection of powerful algebraic packages and can efficiently perform symbolic computation precisely without round-off. But CAS still can output incorrect results because they have no mechanism to check correctness of results mathematically. Errors can arise in mismatching the side conditions of functions, the determination of equality of two expressions, or the definedness of expressions with respect to symbolic parameters [3, 4]. Moreover both numerical and symbolic computations are implemented in unverified huge computer programs, which are bug prone. Thus the analysis of DFT based on the above mentioned techniques cannot be relied upon for safety-critical applications of which even tiny errors can cause loss of human lives.

Formal methods have succeeded in precise analysis and verification of hardware and software systems in the past decades. Theorem proving, one of the formal methods, models the systems in logic and then reasons and verifies the model’s properties using axioms, theorems, and inference rules in the logic. The proof is performed in a theorem prover [5, 6]. The precise analysis and verification of systems employing DFT need the formal implementation of DFT in a theorem prover. Formal methods have been successfully applied in the precise analyses and verifications of hardware and software systems in the past decades. As a formal method, theorem proving can be used to model the systems based on logic and then to reason and verify the model’s properties using axioms, theorems, and inference rules in the logic. Generally, the proving tasks are performed in a theorem prover [5, 6]. For example, one needs to implement DFT formally in a theorem prover to perform the precise analyses and verifications based on DFT. In recent years, some formalization and verification [7–10] have been done on Fast Fourier Transform (FFT), which is the most popular and efficient DFT method. All these works focused on the specific formal implementations of FFT at different levels in diverse algorithms. Bjesse [7] verified FFT hardware at the netlist level. Capretta [8] formalized FFT at the power list level in the Coq theorem prover, while Gamboa [9] performed similar work in ACL2. Akbarpour and Tahar [10] formally verified FFT at different abstraction levels including the real specification level, the floating and fixed-point description level, and the RT and netlist gate levels, and they also analyzed the round-off accumulation errors. Although FFT may be the most widely used implementation of DFT, it cannot take place of DFT completely. DFT has its own theoretical advantages and is much simpler than FFT in both theory and algorithm. Therefore, compared with FFT, DFT has been all along a dominant method in theoretical analysis, especially when theorem proving is used to analyze the functions and properties of a system model. The applications of DFT in theorem proving are very important in the early stages of engineering design. Generally, DFT can be adopted to perform real-time computation on local data in the process of sampling, while FFT computation can only be done when all the data is obtained. Besides, DFT has more advantages over FFT in the choice of sampling rate, the memory consumption of data, and so forth. To the best of our knowledge, no work has been done previously to formalize DFT in any theorem prover. In the present paper, DFT is for the first time formalized in the HOL theorem prover. To some extent, the present work is to develop a theoretical tool rather than to just provide an application. Systematic and detailed formalization of DFT is presented here to support flexible formal analyses of various systems.

The HOL system is a well-developed theorem prover from the University of Cambridge. It is equipped with the formalized mathematic theories needed for DFT such as real analysis [11], complex field [12], and matrix theory [13, 14]. The rest of the paper is organized as follows. After the formal definition of DFT in Section 2, the fundamental properties of DFT are formalized and verified in Section 3. And then in Section 4, two case studies are presented to illustrate the usefulness of the proposed formalization of DFT, including formal analysis of FFT and the cosine frequency shift. Section 5 concludes the paper.

2. Formal Definition of DFT

For a real or complex sequence with length of in time domain, DFT transforms it to a new complex sequence with the same length in frequency domain. The mathematical definition of DFT can be written as follows: where , is the frequency spectrum. For simplification, it is substituted by in the rest of the paper. The definition can be formalized in HOL as follows.

Definition 1. DFT_def,where accumulates the complex sequence from number 0 item for items. The function is the base- exponential function in which stands for the imaginary unit and stands for . The above mentioned symbols keep the same meanings in the rest of this paper.

When , the inverse DFT (IDFT) can be expressed as follows: The IDFT can be formalized as follows.

Definition 2. IDFT_def ,where stands for the DFT sequence of with length of .

3. Formal Verification of the Properties

In this section, the classical properties of DFT are formalized and verified based on Definition 1 in HOL. The formal verification of the properties can check whether or not the definition is correct and reasonable, and moreover the verified properties turn into theorems in HOL and they can be utilized directly to facilitate formal analysis of DFT based systems.

3.1. Some Useful Lemmas

To conveniently read and understand formalizations of properties of DFT, some useful lemmas verified in HOL are presented in this subsection. Consider The formalization is as follows.

Lemma 3. Equality of two complex functions’ cumulative sum      ,where functions and are functions from the natural numbers to complex numbers for convenience to DFT.

Intuitively, is an -periodic complex sequence. It holds that The formalization is as follows.

Lemma 4. Periodicity of is     .

For a complex number , IM returns the imaginary component of . One has

Lemma 5. Inverse operation of imaginary parts of complex numbers .

Suppose is a complex function; it holds that The formalization is as follows.

Lemma 6. Compound operation and accumulating operation  .

For the conjugate operator and accumulating function, there is a similar property as follows: where stands for the conjugate operator of complex numbers. The formalization is as follows.

Lemma 7. Conjugate operation and accumulating operation  ,where the function returns the conjugate of a complex number.

Accumulation of a series of complex numbers can be performed on the real and imaginary parts, respectively. Consider where the functions RE and IM return the real and imaginary components of a complex number, respectively. The formalization is as follows.

Lemma 8. Accumulating operation of complex numbers  ,  .

3.2. Implicit Periodicity

DFT is derived from the Discrete Fourier Series. Considering is an -periodic sequence, suppose is a series with -periodic extension. Then for = DFT(), is periodic. One has The formal description is as follows.

Theorem 9. DFT_PERIODIC      .

3.3. Linearity

For two -length sequences and , and are their DFT sequences. The linearity is as follows: We verified this property as the following theorem.

Theorem 10. DFT_LINEAR    +  .

The linearity of DFT shows that it can be applied in discrete linear system. The proof is based on Definition 1 and Lemma 3.

3.4. Symmetry

The time domain and frequency domain of DFT are conjugated dually and written as The formal description is as follows.

Theorem 11. DFT_SYMM         .

During verifying the theorem, we need to introduce a subgoal: in order to simplify the goal. Then apply some general operation rules in real script and Lemma 3 to make this theorem proved.

3.5. Conjugated Symmetry

The sequence after DFT is of conjugated symmetry. Consider The formal description is as follows.

Theorem 12. DFT_CONJ_SYMM     .

3.6. Frequency Shift

The phase shift of the sequence in the time domain corresponds to the circular shift of the sequence in the frequency domain. The frequency shift contains the left shift and the right shift. They can be used in formal verification of cosine frequency shift. One has The left shift and the right shift are formalized into Theorems 13 and 14, respectively.

Theorem 13. DFT_FREQUENCY_LSHIFT      

Theorem 14. DFT_FREQUENCY_RSHIFT      .

3.7. Time Shift

Because of the duality property between the time domain and the frequency one, the circular shift of in the time domain corresponds to the phase shift of in the frequency domain. Consider The formal description is as follows.

Theorem 15. DFT_TIME_SHIFT      .

3.8. Convolution

The convolution can be obtained as the product’s IDFT of two sequences’ DFT. The property facilitates a remarkable simplification. For two -length sequences and , their circular convolution is defined as follows: The formal definition is as follows.

Definition 16. DFT_CONVOLUTION  .
Let be DFT of ; the circular convolution theorem holds as follows: The formal theorem is as follows.

Theorem 17. DFT_CONVOLUTION_TIME      .

According to the theorem, the convolution can be calculated from the IDFT of . Benefiting from Fast FT algorithms, the computing complexity can be reduced greatly.

4. Applications

To illustrate the correctness and usefulness of our work, we formally verify FFT and cosine frequency shift as examples in this section. DFT always has been implemented by employing FFT, which is one of the most important algorithms. The formal verification of FFT is the basis to verify applications which employed FFT. Frequency shift is a significant property which should be taken into account at DFT applications. The two examples are chosen for their representativeness.

4.1. FFT

Without question, DFT is extremely useful in many fields but it often takes too long time to compute DFT directly by the definition, especially for long date sequences. Therefore, many FFT algorithms had been developed to perform DFT quickly. Directly computing DFT of an -length sequence takes arithmetical operations; by comparison, FFT algorithms take arithmetical operations to produce the same result. According to Formula (1) of DFT, the rotating operator , which is a periodic complex exponential sequence, enables the complexity reduction.

For the periodic complex-power sequence , it holds conjugated symmetry and reducibility as follows: The formalization is as follows.

Lemma 18. WN_CONJ_SYMM ,

The formalization is as follows.

Lemma 19. WN_REDUCE     .

Because of the conjugated symmetry of the rotating operator, some multiplying terms can be merged and the amount of the multiplications reduces by a half roughly; because of the periodicity and the reducibility of the rotating operator, a DFT of a long length sequence can be factorized into multiple DFTs of short sequences, where the shorter the sequence, the lower the computing complexity. Based on the ideas above mentioned, many FFT algorithms had been developed. In this section, we formally verify the well-known radix-2 Cooley-Tukey algorithm.

From the reducibility of , taking -length DFT and , it holds that The formalization is as follows.

Lemma 20. WN_REDUCE_HALF  .

For computing DFT of -length sequence , , the sequence is divided into two length subsequences in terms of the parity of , as follows:

Suppose that ,  ,   ; FFT can be described as follows:

FFT is consisting of the inversion processing and the butterfly shaped computation. And according to the above expressions, the FFT of is the concatenation of Formula (23). The formal description of FFT is shown in Algorithm 1.

FFT:  (k N m:num)(f:num->real) (f1:num->real) (f2:num->real).   (N <> 0) ∧ (N = 2 m) ∧ (0 <=  k ∧ k < N DIV 2) ∧ ((f1(r)  =  f(2 r)) ∧ (f2(r)  =  f(2 r  +  1))) ∧   (DFT f1 ((N DIV 2):num) k = csum (0,N DIV 2) (λr. f1 r exp (i (-2 pi/(N/2) r k)))) ∧   (DFT f2 ((N DIV 2):num) k = csum (0,N DIV 2) (λr. f2 r exp (i (-2 pi/(N/2) r k)))) ∧   (DFT f N k = csum (0,N) (λn. f n exp (i (-2 pi/N n k)))) ∧   (DFT f N (k + (N DIV 2)) = csum (0,N) (λn. f n exp (i (-2 pi/N n (k + (N/2))))))     ⇒     ((DFT f N k = DFT f1 ((N DIV 2):num) k + exp (i (-2 pi/N k))     DFT f2 ((N DIV 2):num) k) ∧     (DFT f N (k + (N DIV 2)) = DFT f1 ((N DIV 2):num) k + exp (i (-2 pi/N (N/2)))       exp (i (-2 pi/N k)) DFT f2 ((N DIV 2):num) k))     ⇒       ((DFT f N k = DFT f1 ((N DIV 2):num) k + exp (i (-2 pi/N k))       DFT f2 ((N DIV 2):num) k) ∧       (DFT f N (k + (N DIV 2)) = DFT f1 ((N DIV 2):num) k - exp (i (-2 pi/N k))        DFT f2 ((N DIV 2):num) k))

Algorithm 1 

The formal definition of FFT.

The FFT can be formally verified as an implementation of DFT. Firstly, use DEF_def and WN_REDUCE_HALF to rewrite the goal. We also need to introduce a subgoal “” to simplify the current goal. Then we can use some tactics in hol4 and some appropriate operation rules in real script and complex script to extend and transform the goal as well as WN_PERIODIC, WN_CONJ_SYMM, WN_REDUCE, and other properties mentioned above. The procedure of validation is shown in Algorithm 2.

REWRITE_TAC [DFT_def, WN_REDUCE_HALF] THEN REPEAT STRIP_TAC THENL   //the initial goal was decomposed into 2 subgoals   [ASM_REWRITE_TAC []   //the 1st  subgoal proved!,   POP_ASSUM MP_TAC THEN   Q.SUBGOAL_THEN “exp (i (-2 pi/N (N/2)))  = -1” ASSUME_TAC   //the 1st  auxiliary subgoal was added!   THENL     [REWRITE_TAC [EXP_IMAGINARY] THEN     Q.SUBGOAL_THEN “-(2) (pi:real)/N (N/2):real = -pi”  ASSUME_TAC     THENL       [REWRITE_TAC[real_div] THEN ONCE_REWRITE_ TAC       [REAL_MUL_ASSOC] THEN REWRITE_TAC [GSYM REAL_NEG_LMUL] THEN       ONCE_REWRITE_TAC [REAL_MUL_SYM] THEN       REWRITE_TAC [GSYM  REAL_MUL_ASSOC] THEN       Q.SUBGOAL_THEN “inv (N:real) N:real = 1:real”  ASSUME_TAC       //the 2nd  auxiliary subgoal was added!       THENL         [REWRITE_TAC [GSYM REAL_1] THEN         MATCH_MP_TAC REAL_MUL_LINV THEN REWRITE_TAC [REAL_0] THEN         METIS_TAC [],         ASM_REWRITE_TAC [] THEN         REWRITE_TAC [GSYM real_1,  REAL_MUL_RID] THEN REWRITE_TAC       [REAL_MUL_ASSOC] THEN Q.SUBGOAL_THEN “inv (2:real) 2:real = 1:real”       ASSUME_TAC         //the 3rd  auxiliary subgoal was added!         THENL           [Q.SUBGOAL_THEN “inv 2:real <> 0:real” ASSUME_TAC           //the  4th  auxiliary subgoal was added!           THENL             [REWRITE_TAC [REAL_INV_1OVER]             THEN SRW_TAC[] [],             REWRITE_TAC [GSYM REAL_1]             THEN MATCH_MP_TAC REAL_MUL_LINV             THEN REWRITE_TAC [REAL_0]             THEN SRW_TAC[] []],           ASM_REWRITE_TAC []           THEN SRW_TAC [] [GSYM real_1,  REAL_MUL_LID]]],         ASM_REWRITE_TAC []         THEN REWRITE_TAC[COS_NEG, SIN_NEG] THEN         REWRITE_TAC [COS_PI, SIN_PI] THEN         REWRITE_TAC [COMPLEX_1] THEN         REWRITE_TAC [complex_of_real, complex_neg] THEN REWRITE_TAC [RE, IM,       FST,  SND]       ]       //all the auxiliary subgoals was proved!     ONCE_ASM_REWRITE_TAC [] THEN     ONCE_REWRITE_TAC [COMPLEX_MUL_COMM] THEN     REWRITE_TAC [GSYM COMPLEX_NEG_LMUL] THEN     REWRITE_TAC [COMPLEX_MUL_LID] THEN     REWRITE_TAC [GSYM COMPLEX_NEG_RMUL] THEN     METIS_TAC [GSYM complex_sub]]     //the second subgoal proved!   ]

Algorithm 2 

The validation of FFT.

4.2. Cosine Frequency Shift

Frequency shift is an important property of DFT, and it can derive many other similar properties like cosine frequency shift. This subsection presents formal verification of cosine frequency shift. A typical cosine frequency shift can be formulated as follows:

To verify the formula above using the frequency shift property of DFT, should be transformed to exponential formula with Euler’s formula as follows.

Lemma 21. EULER_COS:

Its formalization is shown in Algorithm 3.

val EULER_COS = store_thm(“EULER_COS”, “ n:real. (cos n,0) = (1)/(2) (exp (i n) + exp (i (-n)))”, GEN_TAC THEN REWRITE_TAC [EXP_IMAGINARY] THEN REWRITE_TAC [complex_add, RE, IM, FST, SND] THEN REWRITE_TAC [COS_NEG, SIN_NEG] THEN SRW_TAC [REAL_DOUBLE] THEN REWRITE_TAC [complex_scalar_lmul, RE, IM, FST, SND] THEN SRW_TAC [REAL_MUL_ASSOC])

Algorithm 3 

The formalization of EULER_COS.

Based on the definition of DFT and the lemma EULER_COS, Formula (25) can be formally described and verified as in Algorithm 4.

val DFT_MOD_COS = store_thm(“DFT_MOD_COS”, “  (k N m:num)(f:num->real). ((N):real <> (0) ∧ (0 <= k ∧ k < N)) ∧ (DFT f n (k  -  m) = csum(0,N) (∖(n:num). (f n) exp(i ((-2) pi/N n (k - m))))) ∧ (DFT f n (k  +  m) = csum(0,N) (∖(n:num). (f n) exp(i ((-2) pi/N n (k + m))))) ⇒   ((csum(0,N) (∖(n:num). (f n (cos((2) pi/N n (m)),0) exp (i (-2 pi/N n k))))) =  (1)/(2) (DFT f n (k-m)) + (1)/(2) (DFT f n (k  +  m)))” REPEAT STRIP_TAC THEN ASM_REWRITE_TAC [] THEN REWRITE_TAC[EULER_COS] THEN REWRITE_TAC [COMPLEX_ADD_SCALAR_LMUL,  COMPLEX_ADD_RDISTRIB, CSUM_ADD]  THEN Q.SUBGOAL_THEN “  (a b:real) c:complex d:complex. a (b c) d = b (a c d)” ASSUME_TAC //the 1st  auxiliary subgoal was added! THENL   [REPEAT STRIP_TAC THEN REWRITE_TAC[COMPLEX_LMUL_SCALAR_LMUL]  THEN   ONCE_REWRITE_TAC [COMPLEX_SCALAR_MUL_COMM] THEN   GEN_REWRITE_TACLAND_CONV [GSYM COMPLEX_SCALAR_MUL_COMM] THEN   GEN_REWRITE_TAC LAND_CONV [COMPLEX_SCALAR_LMUL] THEN ONCE_REWRITE_TAC [COMPLEX_SCALAR_MUL_COMM] THEN   SRW_TAC [] [GSYM COMPLEX_SCALAR_RMUL]   //the 1st auxiliary subgoal proved!,   ASM_REWRITE_TAC [] THEN REWRITE_TAC [CSUM_RMUL] THEN   REWRITE_TAC [GSYM COMPLEX_ADD_SCALAR_LMUL] THEN   AP_TERM_TAC THEN REWRITE_TAC [real_div, REAL_NEG_LMUL] THEN   REWRITE_TAC [GSYM real_div] THEN   BINOP_TAC   //The initial goal was decomposed into 2 subgoals   THENL     [MATCH_MP_TAC DFT_FREQUENCY_RSHIFT THEN     METIS_TAC []     //The 1st subgoal proved,     MATCH_MP_TAC DFT_FREQUENCY_LSHIFT THEN     METIS_TAC []     The 2nd subgoal proved     ]   ]

Algorithm 4 

The formalization of cosine frequency shift.

In the description of cosine frequency shift as shown in Algorithm 4, is directly expanded by the definition of DFT to facilitate the proving process.

5. Conclusion

DFT is widely used in DSP and linear time invariant (LTI) systems which could be in safety-critical fields where formal methods are expected to be employed to assure accurate and complete verification. This paper presented the formalization of DFT in HOL, including formal definition of DFT, formal verification of properties of DFT, and two case studies demonstrating usefulness and correctness of our work. Our work can be equipped as a library of HOL to facilitate formal verification of systems employing DFT.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors thank Professor Shengzhen Jin for his helpful suggestions. This work was supported by the International Cooperation Program on Science and Technology (2010DFB10930 and 2011DFG13000), the National Natural Science Foundation of China (61170304, 61104035, 61373034, 61303014, and 61472468), the Project of Construction of Innovative Teams and Teacher Career Development for Universities and Colleges Under Beijing Municipality, and the Scientific Research Base Development Program of the Beijing Municipal Commission of Education (TJSHG201310028014).