Fast iterative algorithms for birefringent filter design

Abstract

Fast iterative algorithms for designing birefringent filters with any specified spectral response are proposed. From the Jones formalism, we derive two polynomials representing the transmitted and rejected response of the filter, respectively. Once the coefficients of the filters are obtained, the orientation angle of each birefringent section and the phase shift introduced by each compensator can be determined by an iterative algorithm that gives an efficient solution to the birefringent filter design problem. Afterward, some design examples are presented to demonstrate the effectiveness of the proposed approach. In comparison with results reported in the literature, this approach provides the best performance in terms of accuracy and time complexity.

Citation

Boukharouba, A. (2021), "Fast iterative algorithms for birefringent filter design", Applied Computing and Informatics, Vol. 17 No. 2, pp. 250-263. https://doi.org/10.1016/j.aci.2018.08.004

1. Introduction

Optical filters are most commonly used in modern optical telecommunication systems [1] and also used for biomedical spectral imaging [2–4] such as Raman chemical imaging [5] and fluorescence microscopy applications [6]. Optical filters are essentially based on optical interleavers that separate an incoming spectrum into two complementary set of periodic spectra or combine them into a composite spectrum [7]. Most interleavers are based on Michelson interferometers, Mach–Zehnder interferometer (MZI) or birefringent filter principles [8–10]. The two classical designs of birefringent filters are Lyot filter [11] and Solc filter [12,13], which both types of filters use different configurations of polarizers and retardation plates to create narrow-band filters. Interference birefringent filters are optical finite impulse response (FIR) filters based on the changes induced in the state of polarization of light by birefringent materials. They are composed of a stack of retardation plates of birefringent material placed between a polarizer and an analyser. A wide range of filters can be achieved by orienting birefringent elements in an appropriate way. They play an important role in dense wavelength division multiplexing (DWDM) systems, as in gain equalization, dispersion compensation, prefiltering, and channels add/drop applications. The basic idea for the birefringent filter design is proposed by Harris et al. [14]. The desired output spectrum is first developed into the finite terms of a Fourier series and then the relative angles of both retarders and analyzer are determined according to the backward transfer method. A generalization of the procedure presented in [14] that allows the realization of impulse response having complex coefficients is detailed by Ammann and Yarborough [15]. The resulting network consists of n stages between an input and output polarizer, with each stage containing a birefringent crystal and optical compensator. Afterwards, the main problem of optical filter design becomes a problem of designing FIR filters using the Fourier series expansion of the desired frequency response. This technique is modified and improved by using a windowing technique to improve the shape of the frequency response. Classical optimization methods such as weighting least square sense and Parks-McClellan method are also used for designing digital FIR filters. To improve the performance of the classical methods, many researchers have utilized heuristic evolutionary optimization algorithms such as Genetic Algorithm (GA), Differential Evolution (DE), and Swarm Optimization (SO). For example, an optical finite impulse response (FIR) filter design methods based on crystal birefringence to produce arbitrary spectrum output are presented where a typical example of a green/magenta filter used in a liquid crystal on silicon projection display is synthesized [16,17]. An example of birefringent equalizing filter suitable for dispersion compensation in wavelength division multiplexed (WDM) communication systems is presented in [18]. A backward recursion of the transfer matrix is used to calculate the parameters of an optical filter that has an impulse response with complex coefficients. In [7] a general synthesis method for designing asymmetric flat-top birefringent interleavers is reported using a combination of digital signal processing approach and computational optimization by GAs.

The birefringent filter structure may be synthesized using a different technology based on coherent optical delay-line circuit with a two-port lattice-form configuration [19] where arbitrary filter characteristics corresponding to nth-order complex FIR digital filters can be realized by n cascaded two-port lattice-form optical delay-line circuits.

In this paper, we present iterative algorithms to design a birefringent filter whose coefficients are real or complex numbers and having an arbitrary frequency response. Some design examples are also presented to demonstrate the effectiveness of the proposed algorithms. The paper is organized as follows: Section 2 presents a theoretical analysis of the proposed algorithms for synthesising an optical FIR filter with real coefficients. Next, an extension to an arbitrary frequency response where the filter coefficients are complex numbers is detailed in Section 3. Section 4 introduces further improvements in the method given in [15] to calculate the complementary component. Design examples and simulation results are discussed in Section 5. Discussions and performance comparisons against other existing methods are presented in Section 6 and finally some conclusions are exposed in Section 7.

2. Optical finite impulse response filter

First, we study an optical FIR filter which is composed of a stack of identical birefringent retarders with same length L placed between an input polarizer and output analyser as shown in Figure 1. The x-axis is chosen parallel to the transmission axis of the input polarizer while s and f represent respectively the slow and fast axis of the birefringent elements. The solid arrows represent the fast axes of birefringent retarders and the transmitted axis of the output analyser. φk represents the angle between the fast axis of the kth retarder and the y-axis which is the same angle between the slow axis of the retarder and the x-axis, while φp is the angle between the transmission axis of the output analyzer and the x-axis. After the input polarizer, the polarized light comes through the retarder stack where each retarder separates the input light into two components along its fast and slow axis respectively, and each component acts as the input of the next one [14,17].

The output of the optical filter must give the desired impulse response, A(t).

(1)A(t)=∑k=0nAkδ(t−ka)

where a denotes the time intervals between the impulse series.

(2)a=LΔηc

L is the length of each retarder, Δη represents the difference of their refractive indices, and c is the velocity of light in vacuum, so the phase difference caused by each retarder is expressed by:

(3)Γ=2πΔηLλ

The frequency response of the optical filter is the Fourier transform of its impulse response (1).

(4)A(ω)=∑k=0nAkexp(−jkaω)

where ω=2πf=(2πc)/λ denotes the angular frequency of the light. Hence the spectrum response of the filter is given by:

(5)A(λ)=∑k=0nAkexp(−jkΓ)

The complementary component, which is the output along the perpendicular direction of the analyzer, can be expressed as:

(6)B(λ)=∑k=0nBkexp(−jkΓ)

We assume that the optical network is lossless, which means that the energy must be conserved at all points within the network independently of the optical frequency ω.

(7)|A(ω)|2+|B(ω)|2=I02

The output, Eout, is determined by Jones matrix [16,20].

(8)Eout=∏i=1n(cosθiexp(−jΓ)sinθi−sinθiexp(−jΓ)cosθi)(cosθ0−sinθ0)

Here θi (i≠n) is the relative angle of each retarder and θn is the relative angle of the analyzer as described in [16].

(9)θ0=φ0,θ1=φ1−φ0, …, θn=φp−φn−1

The main idea of this algorithm is to transfer the matrix–vector product into tow polynomials αi(z) and βi(z).

2.1 Forward algorithm

Let’s set ci=cosθi,si=sinθi, and z=exp(−jΓ). Consequently Eout is expressed as follows:

(10)Eout=∏i=1n[cizsi−sizci][c0−s0]

For two stages Eout has the following expression:

(11)Eout=[c1zs1−s1zc1][c0−s0]=[c1c0z−s1s0−s1c0z−c1s0]

By induction, each component of the Eout is a polynomial of z for n≥1. We set α00=c0 and β00=−s0 and for (n−1) stages, Eout may be written as:

(12)Eout=[∑i=0n−1αin−1zi∑i=0n−1 βin−1zi]

For n stages Eout is given by:

(13)Eout=[cnzsn−snzcn][∑i=0n−1 αin−1zi∑i=0n−1 βin−1zi]

We can simply prove that:

(14)Eout=[∑i=0n αinzi∑i=0nβinzi]

where α0n=snβ0n−1, β0n=cnβ0n−1, αnn=cnαn−1n−1, βnn=−snαn−1n−1 and αin and βin can be calculated by the following equations for n≥2 and i=1,2...,n−1.

(15)αin=cnαi−1n−1+snβin−1

(16)βin=−snαi−1n−1+cnβin−1

From Eq. (14) and knowing that the frequency responses of the two complementary filters are defined by A(z)=∑i=0nAizi and B(z)=∑i=0nBizi respectively. Thus, αin=Ai and βin=Bi , for i=0,1,2...,n

Once the relative angles are known, the forward algorithm gives the impulse responses of the optical filters. However, the main problem of the birefringent filter synthesis is how to calculate the relative angles from the coefficients of the desired impulse response. To do this, we use a backward algorithm where the relative angles are the solution of the set of equations of the forward algorithm.

2.2 Backward algorithm

From the forward algorithm, the value of the nth relative azimuth angle θn is given by θn=arctan(−βnn/αnn) and the coefficients of the (n−1) order filter are given by αn−1n−1=αnn/cn and β0n−1=β0n/cn. But for the other coefficients, we must solve the system of Eqs. (15) and (16) where three determinants must be formed; the denominator determinant is Δ=1 and consequently the αi−1n−1 and βin−1 have the following expressions for i=1:n−1.

(17)αi−1n−1=|αinsnβincn|=cnαin−snβin

(18)βin−1=|cnαin−snβin|=snαin+cnβin

The Algorithm 1 summarizes the whole recursive process:

3. Extension to an arbitrary frequency response

In this case the optical network consists of n retarders between a polarizer and an analyser where each retarder is composed of a birefringent crystal with equal length and optical compensator [15]. However, in [18], each retarder is a birefringent plate composed of a section of nominal length L and a section with variable length Li that acts as an optical compensator, see Figure 2. The identical lengths introduce a unitary delay whereas the variable lengths (optical compensators) introduce variable phase shifts bi between slow-axis and fast-axis components. The expression of Eout is the same as in (8) except the delay operator z=exp(−jω) becomes zexp(−jbi) where the phase shift bi is introduced by the variable length Li.

(19)Eout=∏i=1n[cizexp(−jbi−1)si−sizexp(−jbi−1)ci][c0−s0]

3.1 Forward algorithm

Following the same steps as in the Section 2.1, we can prove that:

(20)Eout=[∑i=0nαinzi∑i=0n βinzi]

where α0n=snβ0n−1, β0n=cnβ0n−1, αnn=cnexp(−jbn−1)αn−1n−1, βnn=−snexp(−jbn−1)αn−1n−1 and αin and βin can be calculated by the following equations for n≥2 and i=1,2...,n−1.

(21)αin=cnexp(−jbn−1)αi−1n−1+snβin−1

(22)βin=−snexp(−jbn−1)αi−1n−1+cnβin−1

3.2 Backward algorithm

The backward algorithm determines the relative angle of each crystal, the retardation introduced by each compensator, and the relative angle of the analyser. From the forward algorithm, the value of the nth relative azimuth angle θn is given by θn=arctan(α0n/β0n) and the coefficients of the (n−1) order filters are given by αn−1n−1=(αnn/cn)exp(−jbn−1) and β0n−1=β0n/cn. But for the other coefficients and from the Eq. (21) and (22) αi−1n−1 and βin−1 have the following expressions for n≥2 and i=1,2...,n−1.

(23)αi−1n−1=(cnαin−snβin)exp(jbn−1)

(24)βin−1=snαin+cnβin

As shown in Section 3.1, α0n must be real thus we need to add an optical compensator in the front of the analyser by choosing bp=−angle(A0) and the filter coefficients at the analyser output become αin=Aiexp(jbp) , i=0,1,2...,n. In practice, the final compensator may be removed from the filter structure because the transmittance with and without final compensator differ from each other by only a constant phase factor, and hence the final compensator can be ignored.

Moreover, α0n−1=(cnα1n−snβ1n)exp(jbn−1) must be also real and consequently, we can take α0n−1=|cnα1n−snβ1n| and bn−1=−arg(cnα1n−snβ1n) as solution especially for calculating the phase shift bi. If we count the optical compensator in the front of the analyser, we also have n compensators (bn=bp). Note that if bi=0 , a compensator is not required for this stage. The coefficients Bi are in general complex numbers; we can note that if B(ω) is a solution of (7), then exp(jμ)B(ω) is also a solution [15]. Knowing that β0n is also a real number then we can choose the coefficients of the complementary component as βin=Biexp(−jμ0), i=0,1,2...,n where μ0=angle(B0) to make β0n real. The Algorithm 2 explains the backward algorithm.

4. Complementary component calculation

Assume that A(ω) and therefore the desired Ai have been chosen. We must next find B(ω) which is polarized perpendicular to A(ω) and therefore stopped by the output analyzer.

(25)B(ω)B∗(ω)=I02−A(ω)A∗(ω)

Here I02 must be chosen greater than, or equal to, the maximum value of A(ω)A∗(ω). Having chosen I02, we can calculate B(ω) from B(ω)B∗(ω) using the method given in [15] with further improvements. Letting x=exp(−jaω), the expression of (25) becomes as follows:

(26)|B(x)|2=dn∗xn+dn−1∗xn−1+…..+d1∗x1+d0+d1x−1+…..dn−1x−(n−1)+dnx−n

Note that |B(x)|2 has got the same roots as xn|B(x)|2 which is a polynomial of degree 2n.

(27)xn|B(x)|2=∑k=02nDkxk

where Dk=−∑m=0kAmAn+m−k∗ for k=0:n−1,Dn=I02−∑m=0n |Am|2, and Dn+k=Dn−k∗ for k=1:n. The roots of the polynomial xn|B(x)|2 are simply calculated using the eigenvalues of its companion matrix. It is obvious that if rk is a root of (26), then (1/rk∗) is also a root. One of these two roots is associated with B(x) and the other with B∗(x). if |rk|<1 then |1/rk∗|>1, thus we can retain the roots whose amplitudes less than unity as roots of B(x) as straightforward solution. Once the roots of B(x) are found, we use the same steps as in [15] to get the coefficients of the polynomial B(ω),Bk. if we have p(x)=∏k=1n (x−rk)=∑k=0n pkxk then B(x)=qp(x)=∑k=0n qpkxk where q=D(n)/∑k=0n pk2. Consequently, the coefficients of the polynomial B(ω) are given by Bk=qpk.

Notice that the all above algorithms are developed for I02=1. However, if I0 has an arbitrary value, the output of the optical filter calculated using the orientation angle found by the forward algorithm must be multiplied by I0 to have the desired impulse response of the optical filter. Multiplying by I0, which represents an amplification/attenuation, has no effect on the shape of the filter response. For all the following examples the chosen value of I0 is 2 (I02=2). Consequently, the output of the filter must be multiplied by 2 to have the desired impulse response.

5. Design examples and simulation results

In order to demonstrate the effectiveness of the proposed algorithms for optical filter design, some design examples are presented and discussed.

5.1 Flat-top birefringent interleaver filter

we study the case of an asymmetric flat-top birefringent interleaver synthesized using Parks-McClellan optimal equiripple FIR filter design algorithm. The resulting filter H(z) is a seventh-order filter, and its coefficients Ak are calculated in [7]. However, Bk are calculated using the algorithm stated in Section 4, and the relative angles are found by the backward algorithm of the Section 2 where the filter coefficients are real. Once the orientation angles are found, we can calculate the coefficients A^k of the optical filter using the forward algorithm and compare them with the desired ones. Table 1 shows the filter parameters obtained by the proposed algorithms.

We can notice that the calculated coefficients of the optical filter are exactly equal to the desired ones and consequently the desired filter and the obtained one have the same spectrum.

5.2 Multi-channel selector

The second example is a multi-channel selector which is an optical frequency filter designed to select signals at certain frequencies from eight frequency-division multiplexing (FDM) signals [19]. In this case, the optical filter is synthesized to select three signals with frequencies of f1,f3 and f7 from eight FDM signals. Thus, the transmittance values at the three frequency points of (−7/16)fo,(−3/16)fo,(5/16)f0 must be 1, and the transmission values at the other frequency points must be 0 where the number of expansion coefficients is set at 16. The inverse discrete Fourier transform is used to obtain the filter coefficients. As the coefficients are complex numbers, we must use the algorithms of the Section 3 to calculate the parameters of the filter. Table 2 shows the obtained opto-geometrical parameters of the designed filter. The phase shifts bk are non-zero, because the expansion coefficients are complex numbers. Once the relative angles and phases shifts are obtained using the backward algorithm, we can calculate the output of the optical filter using the forward algorithm where the calculated coefficients A^k match exactly the desired coefficients Ak, which confirms the accuracy of the proposed approach.

5.3 Dispersion compensation

Birefringent equalizing filters are interesting examples of optical filters whose coefficients of the impulse response are complex numbers. They are suitable for dispersion compensation in wavelength division multiplexed (WDM) communication systems. The filter coefficients Ak are calculated in [18]. The coefficients Bk of the complementary component are calculated using the procedure described in Section 4 whereas the orientation angles, the phase shift introduced by each stage and the coefficients A^k of the obtained impulse response are calculated using the algorithms described in Section 3. Table 3 illustrates the coefficients Bk; relative angles θk, the phase shifts bk and the coefficients A^k of the obtained impulse response.

We can also notice that the calculated coefficients A^k match exactly the desired coefficients Ak.

5.4 Particle Swarm Optimization (PSO)

The opto-geometrical parameters of the birefringent filter can also be calculated using the heuristic evolutionary optimization algorithms. The cost function to be minimized is expressed according to the optical filter output αkn calculated by the forward algorithms and the desired output for both above cases.

(28)f(θk, bk)=1n∑k=0n(αkn−Ak)2

The parameters to be determined are the relative angle θk of each crystal and the retardation introduced by each compensator bk. In this work, we can take as example the PSO algorithm, which is the famous one of the heuristic evolutionary optimization algorithms [21]. Suppose that the search space is a n-dimensional space, then the ith particle can be represented by a n-dimensional vector, xi={xi1,xi2,…..,xin}, and velocity vi={vi1,vi2,…..,vin}, where i=1,2,…,p and p denotes the size of the swarm. In each generation t+1, particle i adjusts its velocity vikt+1 and position xikt+1 for each dimension d by referring to the personal best position pbestikt and the swarm’s best position gbestikt as follows [22]:

(29)vikt+1=wvikt+c1r1t(pbestikt−xikt)+c2r2t(gbestikt−xikt)

(30)xikt+1=xikt+vikt

where w is the inertia weight, c1 and c2 are positive constants, called acceleration coefficients, r1 and r2 are two random numbers in the range [0,1]. The algorithm searches for the best solution through an iterative process by minimizing the above cost function. The positions {xik} are the parameters of the optical filter such as the relative angles for filters with real coefficients and the relative angles with their phase shifts for filters with complex coefficients. If the algorithm converges after m iterations, the global best gbestikm is the position that has produced the smallest cost function value of all positions occupied by the swarm through last iteration m. Consequently, we can retain the global best gbestikm where k=0,1,…,n as the optimal filter coefficients. We consider the simplest example of the above three examples, flat-top birefringent interleaver filter, which is a filter with real coefficients that needs to determine only the relative angles θk. Thus, the cost function becomes f(θk)=(1/n)∑k=0n (αkn−Ak)2 where αkn is the output of the optical filter calculated by the first forward algorithm and A={Ak} is the desired output of the first example as indicated in the first column of the Table 1. The PSO parameters are as follows: population size, npop=100; inertia weight, w=0.99; acceleration coefficients, c1=1.5,c2=2.0; lower and upper bound of variables, varmin=−1.5,varmax=1.5; lower and upper bound of velocities velmin=−0.3,velmax=0.3. Table 4 shows the obtained relative angles and their corresponding filter coefficient using PSO with different number of iterations.

Compared to the proposed algorithms, the PSO with 50,000 iterations and for the simplest example it could not even find the exact coefficients of the desired response. As a result, the heuristic evolutionary optimization algorithms are not suitable for birefringent FIR filter design.

6. Discussion

As shown above in the Section 5, simulation results and comparisons with a state of the art methods show that the proposed algorithm is faster, easier and more accurate to calculate the optical structure of the birefringent filters. In [18,7] the procedure of the parameter calculation is complicated, not clear and needs a great number of basic arithmetic operations (addition, subtraction, multiplication and division) compared to our algorithm. For example the phase shift in [18] is expressed as a ratio between two quantities and needs four multiplications, two additions and one division, whereas in our algorithm only two multiplications and one addition are required as expressed above in the second backward algorithm. In addition, the algorithm propose by [7] need two recurrent equations to calculate the coefficients of the (k−1)th stage from ones of the kth stage. However, in the proposed algorithm we have two coefficients (αk−1k−1 and β0k−1) calculated as a ratio between two numbers in each stage and the two last coefficients (α00 and β00) are also expressed as a simple ratio between two numbers as illustrated in the first backward algorithm. Moreover the coefficients Bk are calculated using FIR filter design techniques in [7], which needs a heavy computation and lacks accuracy compared to the proposed algorithm. However, we have developed a concise and accurate algorithm to calculate the coefficients, Bk, of the complementary component. The forward algorithm is used to simulate the optical filter by calculating the filter coefficients from the relative angles and the phase shifts obtained using the backward algorithm and compare them with the desired ones. It may be also used with the classical optimisation methods and the heuristic evolutionary optimization algorithms where the cost function to be minimized is a function of the filter parameters, which is easily calculated using the forward algorithm. The proposed approach is also compared with the PSO algorithm, which is the famous one of the heuristic evolutionary optimization algorithms. As illustrated in Table 4, the PSO with 50,000 iterations and for the simplest example it could not even find the exact coefficients of the desired response. Moreover, the evolutionary optimization methods are realized by multiple iterations of updating parameters until convergence. Consequently, they are also not competent in designing optical FIR filters with complex coefficients where the phase shifts are included due to their heavy computation and unguaranteed convergence. On the other hand, the proposed algorithm does well in designing optical FIR filters of any desired spectral shape and with any order using a simple iterative calculation with a minimum number of arithmetic operations.

7. Conclusion

In this paper, iterative algorithms for designing optical FIR filters with any specified spectral response have been presented. They have been tested using different examples and it is observed that they provide exact results in many applications such as asymmetric flat-top birefringent filter, multi-channel selector, and dispersion compensation in wavelength division multiplexed (WDM) communication systems. However, for PSO based birefringent filters, the algorithm must be run many times with a large number of iterations to obtain good results. Moreover, the evolutionary optimization algorithms are extremely sensitive to starting points and the objective function is multimodal and highly non-lineare, which make them very expensive in terms of execution time. In the proposed algorithms, such complicated problem is reduced to find only the roots of polynomial of degree 2n and the exact solution is determined using iterative algorithms with only a few number of operations.

Finally, knowing that liquid crystal tunable filters are used in optical telecommunication systems and they are also used in multispectral and hyperspectral imaging systems because of their high image quality and rapid tuning over a broad spectral range. Consequently and as a future work, we will try to replace the variable sections of the filter structure with liquid crystal cells whose birefringence can be controlled and tuned with a small voltage. In this way and keeping the same filter structure, we propose to synthesize liquid crystal tunable filters by tuning only the birefringence of the liquid crystal using the same iterative algorithms.

Figures

Figure 1

Basic configuration of stack of birefringent retarders with same length.

Figure 2

Generalized structure of an optical FIR filter.

Table 1

Filter parameters obtained by the proposed algorithms, for k=0:7.

Ak	Bk	θk(rad)	A^k
0.0197	0.0003	0.0169	0.0197
−0.0587	0.0019	−0.0469	−0.0587
0.0043	−0.0027	−0.0067	0.0043
0.5347	−0.0183	0.4293	0.5347
0.5347	0.0321	0.4293	0.5347
0.0043	0.0572	−0.0067	0.0043
−0.0587	−0.2362	−0.0469	−0.0587
0.0197	1.1662	−1.5539	0.0197

Table 2

Obtained opto-geometrical parameters of the designed filter, for k=0:15.

Ak	Bk	θk (rad)	bk (rad)	A^k
0.1875	−0.0081 + 0.0034i	−0.1396	2.7489	0.1875
0.0577 − 0.0239i	0.0039 − 0.0039i	0.0504	2.7489	0.0577 − 0.0239i
−0.0442 + 0.0442i	−0.0009 + 0.0022i	0.0494	2.7489	−0.0442 + 0.0442i
0.0239 − 0.0577i	−0.0000 + 0.0357i	0.0481	−0.3927	0.0239 − 0.0577i
−0.1875i	0.0085 + 0.0205i	0.1524	−0.3927	−0.1875i
−0.0239 − 0.0577i	−0.0130 − 0.0130i	0.0519	2.7489	−0.0239 − 0.0577i
0.0442 + 0.0442i	0.0134 + 0.0056i	0.0519	2.7489	0.0442 + 0.0442i
−0.0577 − 0.0239i	0.0759 + 0.0000i	0.0516	−0.3927	−0.0577 − 0.0239i
−0.1875	0.0337 − 0.0140i	0.1539	−0.3927	−0.1875
−0.0577 + 0.0239i	−0.0229 + 0.0229i	0.0493	2.7489	−0.0577 + 0.0239i
0.0442 − 0.0442i	0.0108 − 0.0261i	0.0504	2.7489	0.0442 − 0.0442i
−0.0239 + 0.0577i	−0.0000 − 0.1193i	0.0512	−0.3927	−0.0239 + 0.0577i
0.1875i	−0.0196 − 0.0474i	0.1437	−0.3927	0.1875i
0.0239 + 0.0577i	0.0335 + 0.0335i	0.0434	2.7489	0.0239 + 0.0577i
−0.0442 − 0.0442i	−0.0400 − 0.0165i	0.0452	−0.3927	−0.0442 − 0.0442i
0.0577 + 0.0239i	1.3342	1.5240	0	0.0577 + 0.0239i

Table 3

Coefficients Bk, relative angles θk, the phase shifts bk and the obtained coefficients Âk, for k=0:20.

Ak	Bk	θk (rad)	bk (rad)	A^k
2.696110−3+4.276510−4 i	−7.378810−6+8.126510−21 i	−0.0027	−1.5995	2.696110−3+4.276510−4 i
−1.653810−3+1.279210−2 i	−2.003510−6+2.181310−20 i	0.0128	1.5686	−1.653810−3+1.279210−2 i
−4.499510−2−5.718710−3 i	8.034710−5+3.062310−19 i	0.0449	1.5731	−4.499510−2−5.718710−3 i
1.561710−2−1.206910−1 i	1.633610−5+2.631310−18 i	0.1196	1.5708	1.561710−2−1.206910−1 i
2.475810−1+3.203410−2 i	−2.775110−4+1.411410−17 i	0.2378	1.5692	2.475810−1+3.203410−2 i
−4.760610−2+3.723310−1 i	−3.705310−5+4.955210−17 i	0.3239	1.5693	−4.760610−2+3.723310−1 i
−3.561410−1−4.502710−2 i	1.107610−4+1.105110−16 i	0.2429	−1.5755	−3.561410−1−4.502710−2 i
1.159310−2−9.116510−2 i	−3.630210−5+2.078710−16 i	0.0270	1.5739	1.1593(10−2−9.116510−2 i
−2.487610−1−3.095210−2 i	2.103910−4+5.678410−16 i	0.2468	1.5740	−2.487610−1−3.095210−2 i
2.902910−2−2.301510−1 i	−1.666610−5+1.608510−15 i	0.1082	−1.5751	2.902910−2−2.301510−1 i
−1.766210−1−2.168410−2 i	1.808910−4+3.053110−15 i	0.2125	1.5751	−1.766210−1−2.168410−2 i
2.902910−2−2.301510−1 i	1.073810−5+3.557210−15 i	0.1082	−1.5740	2.902910−2−2.301510−1 i
−2.487610−1−3.095210−2 i	3.595310−5+3.451510−15 i	0.2468	−1.5739	−2.487610−1−3.095210−2 i
1.159310−2−9.116510−2 i	2.942710−5+5.171610−15 i	0.0270	1.5755	1.159310−2−9.116510−2 i
−3.561410−1−4.502710−2 i	−2.953810−4+8.336410−15 i	0.2429	−1.5693	−3.561410−1−4.502710−2 i
−4.760610−2+3.723310−1 i	−1.127310−5+8.880410−15 i	0.3239	−1.5692	−4.760610−2+3.723310−1 i
2.475810−1+3.203410−2 i	9.695610−5+3.668510−15 i	0.2378	−1.5708	2.475810−1+3.203410−2 i
1.561710−2−1.206910−1 i	−3.925410−7−1.716910−15 i	0.1196	−1.5731	1.561710−2−1.206910−1 i
−4.499510−2−5.718710−3 i	−2.146910−5−6.727210−16 i	0.0449	−1.5686	−4.499510−2−5.718710−3 i
−1.653810−3+1.279210−2 i	−5.958110−7+2.536710−15 i	0.0128	1.5995	−1.653810−3+1.279210−2 i
2.696110−3+4.276510−4 i	1.0099	1.5681	0.1573	2.696110−3+4.276510−4 i

Table 4

filter coefficient using PSO with different number of iterations, for k=0:7.

Ak	Iterations = 10,000		Iterations = 30,000		Iterations = 50,000
	θk (rad)	A^k	θk (rad)	A^k	θk (rad)	A^k
0.0197	0.4315	0.0197	1.0926	0.0197	1.4735	0.0197
−0.0587	−1.1422	−0.0588	1.1563	−0.0587	0.0669	−0.0587
0.0043	0.6651	0.0039	−0.5953	0.0042	0.3078	0.0042
0.5347	0.7994	0.5341	−0.8517	0.5345	−1.0761	0.5346
0.5347	0.8004	0.5341	−0.8517	0.5345	−1.0761	0.5346
0.0043	0.6637	0.0040	−0.5953	0.0042	0.3079	0.0042
−0.0587	−1.1424	−0.0588	1.1563	−0.0587	0.0670	−0.0587
0.0197	−1.1370	0.0198	−0.4783	0.0197	−0.0973	0.0197