Design a linear-phase FIR filter approximating the magnitude of a given IIR filter
What you do in step 1 is simply truncate the infinite impulse response to approximate it by an FIR filter. If you use sufficiently many filter taps, the approximation becomes arbitrarily accurate. This means that the resulting FIR filter approximates the magnitude and the phase characteristic of the original IIR filter. So with this approach the phase will never become linear.

Making the impulse response symmetric to obtain phase linearity, as you do in step 2, will of course change the magnitude response.

What you should do is use the magnitude of the IIR filter as a desired response in a (linear-phase) FIR filter design routine. In that case you will get an FIR filter with an exactly linear phase and with a certain magnitude approximation error. That magnitude error can be made sufficiently small by choosing an appropriate filter order. The simplest approach is probably to use a least squares approximation, which just involves solving a system of linear equations.

Example: I use a peaking EQ filter as the IIR prototype. The coefficients are (b are the numerator coefficients, a are the denominator coefficients):


b = [1.2223e+00, 0, 7.7775e-01];
a = [1.1250e+00, 0, 8.7502e-01];


You can use the magnitude of the IIR filter's frequency response and combine it with a linear phase to obtain the desired response for the FIR filter design routine (N is the filter length). The code is Matlab/Octave syntax:


[H,w] = freqz(b,a,256);
N = 61;
D = abs(H).*exp(-1i*w*(N-1)/2);


You can use a least squares FIR filter design routine called lslevin.m, which you can find here.


h = lslevin(N,w,D,ones(length(w),1));
Hh = freqz(h,1,256);


The figure below shows the magnitudes of the two frequency responses (IIR and FIR):