Fast iterative algorithms for birefringent filter design

Fast iterative algorithms for designing birefringent filters with any specified spectral response are proposed. FromtheJonesformalism,wederivetwopolynomialsrepresentingthetransmittedandrejectedresponseofthe filter,respectively.Oncethecoefficientsofthefiltersareobtained,theorientationangleofeachbirefringentsectionandthephaseshiftintroducedbyeachcompensatorcanbedeterminedbyaniterativealgorithmthat givesanefficientsolutiontothebirefringentfilterdesignproblem.Afterward,somedesignexamplesarepresentedtodemonstratetheeffectivenessoftheproposedapproach.Incomparisonwithresultsreportedin theliterature,thisapproachprovidesthebestperformanceintermsofaccuracyandtimecomplexity.


Introduction
Optical filters are most commonly used in modern optical telecommunication systems [1] and also used for biomedical spectral imaging [2][3][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][9][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.

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.f 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 f 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].

Fast iterative algorithms
The output of the optical filter must give the desired impulse response, AðtÞ.
where a denotes the time intervals between the impulse series.
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: The frequency response of the optical filter is the Fourier transform of its impulse response (1).
where ω ¼ 2πf ¼ ð2πcÞ=λ denotes the angular frequency of the light.Hence the spectrum response of the filter is given by: The complementary component, which is the output along the perpendicular direction of the analyzer, can be expressed as: 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 ω.
The output, E out , is determined by Jones matrix [16,20].Here θ i (i ≠ n) is the relative angle of each retarder and θ n is the relative angle of the analyzer as described in [16].
The main idea of this algorithm is to transfer the matrix-vector product into tow polynomials α i ðzÞ and β i ðzÞ.

Forward algorithm
Let's set c i ¼ cos θ i ; s i ¼ sin θ i , and z ¼ expð−jΓÞ.Consequently E out is expressed as follows: For two stages E out has the following expression: By induction, each component of the E out is a polynomial of z for n ≥ 1.We set α 0 0 ¼ c 0 and β 0 0 ¼ −s 0 and for ðn − 1Þ stages, E out may be written as: For n stages E out is given by: We can simply prove that: where α n i and β n i can be calculated by the following equations for n ≥ 2 and i ¼ 1; 2:::; n − 1.
From Eq. ( 14) and knowing that the frequency responses of the two complementary filters are defined by AðzÞ ¼ P n i¼0 A i z i and BðzÞ ¼ P n i¼0 B i z i respectively.Thus, α n i ¼ A i and β n i ¼ B i , 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.

Backward algorithm
From the forward algorithm, the value of the nth relative azimuth angle θ n is given by θ n ¼ arctanð−β n n =α n n Þ and the coefficients of the ðn − 1Þ order filter are given by α 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 The Algorithm 1 summarizes the whole recursive process: ACI 17,2

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 L i 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 b i between slow-axis and fast-axis components.The expression of E out is the same as in ( 8) except the delay operator z ¼ expð−jωÞ becomes zexpð−jb i Þ where the phase shift b i is introduced by the variable length L i .

Forward algorithm
Following the same steps as in the Section 2.1, we can prove that: where α n n−1 and α n i and β n i can be calculated by the following equations for n ≥ 2 and i ¼ 1; 2:::; n − 1.

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ðα n 0 =β n 0 Þ and the coefficients of the ðn − 1Þ order filters are given by α But for the other coefficients and from the Eq. ( 21) and (22) α n−1 i−1 and β n−1 i have the following expressions for n ≥ 2 and i ¼ 1; 2:::

Fast iterative algorithms
As shown in Section 3.1, α n 0 must be real thus we need to add an optical compensator in the front of the analyser by choosing b p ¼ −angleðA 0 Þ and the filter coefficients at the analyser output become α n i ¼ A i expðjb p Þ , 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.
Þ must be also real and consequently, we can take α n−1 as solution especially for calculating the phase shift b i .If we count the optical compensator in the front of the analyser, we also have n compensators (b n ¼ b p ).Note that if b i ¼ 0 , a compensator is not required for this stage.The coefficients B i 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 β n 0 is also a real number then we can choose the coefficients of the complementary component as

Complementary component calculation
Assume that AðωÞ and therefore the desired A i have been chosen.We must next find BðωÞ which is polarized perpendicular to AðωÞ and therefore stopped by the output analyzer.
Here I 2 0 must be chosen greater than, or equal to, the maximum value of AðωÞA*ðωÞ.Having chosen I 2 0 , 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: Note that BðxÞj 2 has got the same roots as x n BðxÞj 2 which is a polynomial of degree 2n. x where The roots of the polynomial x n jBðxÞj 2 are simply calculated using the eigenvalues of its companion matrix.It is obvious that if r k is a root of (26), then (1=r * k ) is also a root.One of these two roots is associated with BðxÞ and the other with B*ðxÞ.if jr k j < 1 then j1=r * k j > 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ðωÞ; B k .if we have pðxÞ ¼ Consequently, the coefficients of the polynomial BðωÞ are given by B k ¼ qp k .
Notice that the all above algorithms are developed for I 2 0 ¼ 1.However, if I 0 has an arbitrary value, the output of the optical filter calculated using the orientation angle found by the forward algorithm must be multiplied by I 0 to have the desired impulse response of the optical filter.Multiplying by I 0 , which represents an amplification/attenuation, has no effect on the shape of the filter response.For all the following examples the chosen value of I 0 is ffiffi ffi 2 p (I 2 0 ¼ 2).Consequently, the output of the filter must be multiplied by ffiffi ffi 2 p to have the desired impulse response.

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.

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 A k are calculated in [7].However, B k 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 b 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.

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 f 1 ; f 3 and f 7 from eight FDM signals.Thus, the transmittance values at the three frequency points of ð−7=16Þf o ; ð−3=16Þf o ; ð5=16Þf 0 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 b k 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 b A k match exactly the desired coefficients A k , which confirms the accuracy of the proposed approach.

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 A k are calculated in [18].The coefficients B k 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 b A k of the obtained impulse response are calculated using the algorithms described in Section 3. Table 3 illustrates the coefficients B k ; relative angles θ k , the phase shifts b k and the coefficients b A k of the obtained impulse response.We can also notice that the calculated coefficients b A k match exactly the desired coefficients A k .

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 α n k calculated by the forward algorithms and the desired output for both above cases.
The parameters to be determined are the relative angle θ k of each crystal and the retardation introduced by each compensator b k .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 ndimensional vector, x i ¼ fx i1 ; x i2 ; . . .::; x in g, and velocity v i ¼ fv i1 ; v i2 ; . . .::; v in g, where i ¼ 1; 2; . . .; p and p denotes the size of the swarm.In each generation t þ 1, particle i adjusts its velocity v tþ1 ik and position x tþ1 ik for each dimension d by referring to the personal best position pbest t ik and the swarm's best position gbest t ik as follows [22]: where w is the inertia weight, c 1 and c 2 are positive constants, called acceleration coefficients, r 1 and r 2 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 fx ik g 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 gbest m ik 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 gbest m ik 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 is the output of the optical filter calculated by the first forward algorithm and A ¼ fA k g 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.

ACI 17,2
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.

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−1 k−1 and β k−1 0 ) calculated as a ratio between two numbers in each stage and the two last coefficients (α 0 0 and β 0 0 ) are also expressed as a simple ratio between two numbers as illustrated in the first backward algorithm.Moreover the coefficients B k 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, B k , 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.

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.

Figure 1 .
Figure 1.Basic configuration of stack of birefringent retarders with same length.

Figure 2 .
Figure 2. Generalized structure of an optical FIR filter.

4 i 3 .
Table Coefficients B k , relative angles θ k , the phase shifts b k and the obtained coefficients Âk , for k ¼ 0 : 20.