In order to implement fractional order transfer function (FOTF) model in digital systems, transfer function discretization methods are developed to obtain discrete filters that provide acceptable amplitude or phase response approximations to FOTFs in operating frequency ranges. This paper presents an approach to improve discrete IIR filter approximations to FOTFs according to application requirements by using Particle Swarm Optimization (PSO). For this purpose, particles of PSO are initialized around the solution of a well-known analytical approximation method and then, these particles search stable IIR filter solutions in the filter design space to improve magnitude and phase response approximation performance over a desired frequency region. Specifically, PSO algorithm is used to fine-tune results of Tustin recursive discretization method according to a weighed multi-objective cost function. This cost function is expressed as the weighted sum of magnitude and phase response matching objectives with a frequency weighting. The frequency weighting is employed for prioritization of frequency regions in optimization. The proposed method also ensures the stability of optimized IIR filter approximations by enforcing particles to move into stable filter solution regions of the search space. This design approach can contribute to the digital realization of fractional order systems for practical applications.