Design of Fading Channel Simulator Based on IIR Filter Using Genetic Algorithm

This article designs a Rayleigh fading channel simulator using IIR filter (called Doppler filter), which is used to approximate the Jakes Doppler spectrum. We mainly focus on the design of Doppler filter and model it as an optimization problem. Non-convexity of the problem is proved in this paper and genetic algorithm (GA) is used to optimize it, which has never been used before in this area. Simulation result shows that GA converges at the solution corresponding to very accurate approximation of Jakes PSD. Finally, several statistical characteristics of the simulator are verified, including correlation functions and level-crossing rate (LCR), all of which match the theoretical predictions pretty well.


Introduction
Channel simulator is usually used in laboratory to verify system performance, which is convenient and efficient.This article designs an accurate and highly efficient channel simulator to produce fading process ( ) h t .Here ( ) h t is a com- plex Gaussian random process whose amplitude | ( ) | h t is a Rayleigh process, used to characterize the Rayleigh fading channel.
According to [1], when receiver moves relatively to the transmitter in a uniform scattering environment, Doppler shift will occur, which makes correlation functions of ( ) h t satisfy: ( ) 0 where D f is the maximum Doppler frequency.I h and Q h are the real and imaginary parts of ( ) h t respectively and they are independent.0 J is the ze- ro-order Bessel function of the first kind and its Fourier transform (2) is the power spectral density (PSD) of ( ) h t , which is called the Jakes Doppler spec- trum [1].
Channel simulator is designed to generate the fading process ( ) h t whose correlation functions and PSD satisfy (1) and ( 2) respectively.In the open literature, two methods have been proposed to design the simulator.
The first method, called "the sum of sinusoids", was first proposed by Jakes [2], where a number of sine waves were summed up to simulate the fading process.Though the method is simple and easy to be implemented, fading process generated by this method lacks of randomness and some of its high order statistic characteristics are difficult to be consistent with theoretical ones [3].
The second method is to filter noise, where complex white Gaussian noises pass through a filter to generate the fading process.This filter, which we call Doppler filter, can be implemented by both FIR [4] and IIR filter [5].Since the transition band is very steep in Jakes PSD (2), if using FIR filter to approximate it accurately, filter order needs to be high, which is expensive in terms of hardware resources.The IIR filter, however, can achieve excellent approximation with considerably low order.
Reference [5] proposed a simulator structure consisting of an IIR filter and an interpolator.It calculated coefficients of IIR filter using ellipsoid algorithm (EA) [6].To keep the filter stable, any pole gone out of the unit circle in iterations would be reflected back.This makes the iteration start from a new point, like random searching.
Reference [7] [8] also used EA but proposed a barrier function to constrain poles lying inside the unit circle.However, the convergence will be affected by the starting point, because non-convexity of the problem.To the best of my knowledge, previous articles used cascaded low order IIR filter only mentioned the non-convexity of their objective functions, however, without proving.
Based on the review above, this article realizes a channel simulator proposed in [5] using different algorithm.We focus on the design of Doppler filter inside it, which is actually an IIR filter fitting problem.Non-convexity of the problem is proved for the first time and genetic algorithm is used to solve it, which is a kind of global optimization method and performs well in solving non-convex problems.
This paper is organized as follows.Section 2 introduces the structure of the channel simulator and defines the optimization problem to design the Doppler filter.Section 3 uses genetic algorithm to solve the problem and section 4 testifies correlation functions and LCR of the fading process generated.Finally, Section 5 concludes the paper.
Notation: The following notations are used throughout this paper.Boldface lowercase letter denotes vectors.Operator | | denotes the absolute value and 2 x represents Euclidean norm of vector x .The superscript * and T de- note complex conjugation and transposition respectively.

Structure of the Simulator
Figure 1 shows the structure of the channel simulator used in this paper, which consists of three parts.Noise generator produces complex white Gaussian noise.Doppler filter filters white noise to generate correlated noise with its PSD complying with Jakes PSD.Interpolator makes Doppler rate of the fading process output satisfy the requirement of the physical channel.The Doppler rate is defined as , where s f is the sampling rate of fading process ( ) h t [5].
The ρ value at the Doppler filter usually takes 0.2.When implemented, Doppler filter remains unchanged and interpolation factor I of the interpola- tor is changed to adjust ρ at the output: . This article will focus mainly on the design of the Doppler filter next.
One important consideration in designing Doppler filter is to ensure the stability.Generally, the larger the filter order is, the more difficult to keep it stable.On account of this, a good choice is to use the cascaded low order IIR sections as the Doppler filter.Except for being easy to be stable, this structure has additional advantage that is less sensitive to quantization error compared with high order IIR filter.
Formula (3) shows the system function of the traditional biquad IIR filter, where a, b, c, d are real-valued optimization variables.Reference [9] studied the stable condition of this biquad, which was called the "Triangle Stability Region" but not intuitive.
We factorize its numerator and denominator, writing it into the form of conjugate poles and zeros (4), where zr, zi, pr and pi are real-valued optimization variables, representing the real and imaginary parts of zero and pole respectively.

SOS z zero z zero H z z pole z pole zero zr j zi pole pr j pi
This Second-Order Section (SOS) is inherently stable, as long as all poles are strictly inside the unit circle.When implementing in a physical system, restore the conjugate form into its traditional biquad form.
Usually, multiple of SOSs need to be cascaded to approximate the objective response.Assuming we use K SOSs, and uniformly sample variable z into 1 M + points points along the unit circle: / , 0,1,..., (usually M = 500), to get the frequency response of Doppler filter after frequency sampling: wherein A is the scaling factor and x is a vector containing all optimization variables: [ , , , , ,..., , , , ] A zr zi pr pi zr zi pr pi = x (7)

Optimization Model
Ideal amplitude response of Doppler filter is the square root of Jakes PSD (2) but cannot be realized with a physical system according to Paley-Winner theorem, because at the frequency point D f f = the amplitude of Jakes PSD changes from +∞ into zero immediately.For this reason, [4] gives a modified objective response that can be achieved by a stable system.Referring to it herein, this article uses (8) as the objective amplitude response, where i is the index of frequency point i z .

/10
1 / 1 ( / ) , 0,1,..., 1 In the formula above, parameter SA denotes the stopband attenuation relative to the passband.Objective function could be defined as: where is a weight vector to emphasize the error minimization for certain frequency band [7].The complete optimization problem is:   the Doppler filter will be a minimum phase system and have the minimum group delay.Non-convexity of ( 10) confirms only if one point 0 x exists in the feasible re- gion, around which the objective function is not convex.We will find such 0 x next.
According to the first-order conditions of convexity, ( ) f x is convex if and only if the feasible region ( dom f ) is convex and ( 11) holds for all , dom f + ∆ ∈ x x x [10].
Here f ∇ denotes the gradient of ( ) f x respect to x , and ∆x is the in- cremental step of x .
For problem (10), dom f is obviously a convex set.To show the non-con- vexity of ( ) f x , we randomly generate 0 x from dom f and approximately calculate 0 ( ) f ∇ x using its finite difference in computer.Then change ∆x gradually and plot ( ) f x and its first-order approximation to see if equation ( 11) is satisfied.
An example is showed in Figure 2, in which we see that ( ) f x is not convex around 0 x .Indeed, much more than one points like 0 x could be found in si- mulation, at which ( ) f x is not convex.

Introduction of Genetic Algorithm
Genetic Algorithm (GA) simulates the evolution process of biotic population with selection, crossover and mutation operations.The population evolves after several generations to achieve a stable situation, which corresponds to the convergence of optimization problem [11].Procedures in simulating GA have been described as pseudocodes in Figure 3.Some terminologies of GA will be used in the subsequent content: A possible solution of the optimization problem is called an individual in GA.Each individual consists of a number of genes and each gene corresponds to an optimization variable.For example x in ( 7) is called an individual, and variables A, zr, zi, pr and pi in x are genes.Technical details of the pseudocodes will be given in the following.

Solving the Problem Using GA
GA parameters: Supposing that Doppler filter uses K = 7 cascaded SOSs, the referenced GA parameters are listed in.Table 1.Generally, the larger Nind is, the greater probability that the global optimal solution may be contained in the population, so Nind should be configured big enough to keep a good diversity of the genes.
Initializing population: This step generates the initial population and each member must be feasible, meaning that individuals must satisfy constraints in (10).Firstly, generate genes A, zr, zi, pr and pi from uniform distribution in the range [−1, 1] as an individual.
Selection: Selection is to select Nparent individuals from each subpopulation to produce offspring.Individuals stronger will have more opportunities to be selected and mate, while those with poor adaptability will be eliminated gradually.
All individuals selected will be divided into two parts, where the first part will do crossover operation, and another part will mutate.In crossover, 2 Nxov ⋅ individuals need to be selected as parents because two parental individuals will produce one child.As for mutation, Nmut parents are needed because one single parent mutates to produce one child.
In this paper, a tournament selecting method [12] has been used.S (tournament size, e.g. 4) players are randomly chosen from each subpopulation in the tournament.A player is an individual and the one that has the best fitness value will be selected.Repeat this process until Nparent parental individuals are selected in each subpopulation.
Crossover: Crossover simulates the mating process of natural population.In the operation, parental individuals are paired first and each pair consists of a father and a mother.Then randomly select half of genes from the father, and the complementary half from the mother.Finally combine these two parts of genes together to produce a child.
, ; 0 , i Gen i i Gen i Gen i father rnd child mother rnd The subscript 1, , 4 + is the index of gene in the individual and rnd is a random integer generated from 0 and 1. Totally Nxov children in each sub-population are need to be produced in this step.
Mutation: In a biotical population, gene mutation happens when passing genes from parents to children, which will produce new genes and increase the genetic diversity of the population.Gaussian mutation [13] is one kind of classical mutation methods in GA, where a random variable, generated from the Gaussian distribution (16) with zero mean and standard variance σ , will be superimposed on the original gene as a mutation variation.At the beginning of the GA, each gene is dispersedly distributed in the entire feasible region.The variation magnitude can be relatively big at this moment to offer more possibilities for each gene.However, as the evolution progresses, each gene begins to converge and variation magnitude should be reduced accordingly.Because after evolution, gene values that were far away from their values, now have been eliminated already by natural selection.If the variation amplitude generated is too large, the gene newly mutated will still be eliminated with great probability in the subsequent evolution, which indicates that such mutation has no meaning.
This article chooses the standard variance of x to be σ of the Gaussian distribution.As the iteration progresses, x will gradually converge, σ track- ing this tendency becomes smaller.Thus the variation magnitude will also decrease, which avoids invalid searching.
Combination: After mutation, all offspring have been produced and they will be combined to generate a new population, which consists of three parts.The first one is the elite individuals retained from the population of the last generation.The second one is the offspring created by crossover operation and the last part is the children generated by mutation, whose genes are mutated versions of their parents' genes.
Population newly generated will be delivered to calculate fitness again and the generation number Gen will be updated.So the evolution continues and optimization variables converge to the global or near global optimum solution gradually.
Migration: Migration occurs every MigG generations past promoting the communication between subpopulations.Elites of each subpopulation are gathered up first in this step, then they are randomly divided into Nsub parts, passed to each subpopulation to replace their worst fit individuals.Subpopulation and migration mechanism can efficiently improve the probability of getting the global optimal in GA.

Simulation Results
Simulating GA described above in Matlab, we get the Doppler filter.Figure 4 shows the amplitude response of the designed Doppler filter.We see that it matches closely with the targeted response in passband, the transition band is very sharp and the response in stopband fluctuates slightly around 60 dB SA = − .There are two criterions to measure the statistical accuracy of fading channel simulator.The first one is the correlation function.As an example, Figure 5 shows the correlation characteristics of the generated fading process and we see that both auto-correlation and cross-correlation functions simulated match well with theoretical curves.
Another statistic characteristic of the fading channel is the LCR, which indicates the speed of the envelope of ( ) h t passing through level Z in the down- ward direction.For Rayleigh fading channel, LCR satisfies (17), where / rms Z Z λ = is the signal level being crossed normalized with the rms of the fading process [1].   ) and another one is when λ is small ( / 10 rms Z Z < ).However, [5] explained that these two situations and both caused by sparse sampling of the fading process, unrelated to the particular simulation techniques used.
retical ones, expect for two situations where the LCR values are underestimated.However, according to [5], these situations are common problems and unrelated to the particular simulation techniques used.We see that mismatches disappear as interpolation factor I increases.

Conclusion
This paper designs the channel simulator producing correlated Gaussian noises as the Rayleigh fading waveforms.An optimization model has been built to design the Doppler filter, and its non-convexity is proved with the aid of computer simulation.Genetic algorithm is used to solve the problem and it converges to the solution corresponding to very accurate approximation of Jakes PSD.Also statistical characteristics of the generated fading process match well with theoretical predictions, which confirms the validity of the simulator we designed.

Figure 1 .
Figure 1.Structure of the channel simulator.
th k SOS located inside the unit circle to keep Doppler filter stable. Constraint lying inside the unit circle too, so that
feasibility, add the individual into the initial population if its poles and zeros are all within the unit circle, otherwise discard.Repeat this process until Nind feasible individuals are created to form an initial subpopulation.Calculating fitness: Fitness is used to quantitatively describe each individual's adaption to the environment.This article chooses the value of objective function as fitness and the value smaller represents better adaptation.this step, compute the objective function value of individuals in each subpopulation first, then sort all individuals in an ascend order according to their fitness values.

Figure 6
Figure6shows the simulated LCRs, which are in good agreement with theo-

Figure 6 .
Figure 6.Level-crossing rate with interpolation factor 10,100,1000 I = respectively.The simulated results match well with theoretical ones, expect for two situations.The first one is when the channel changes extremely fast ( 0.02 ρ =) and another one is when λ is small (

Table 1 .
Parameters of GA.