filter design - Convolution with an impulse response does not give the same result as the frequency response

I am trying to design a linear-phase FIR filter with the following frequency and phase responses:

Given its characteristics (peak filter at 1kHz with 9dB gain and Q=7), I got the following IIR biquad coefficients:

$$b_0 = 1.0889938082113808$$ $$b_1 = -1.8806408829603927$$ $$b_2 = 0.80787505079350063$$ $$a_0 = 1$$ $$a_1 = -1.8806408829603927$$ $$a_2 = 0.89686885900488144]$$

I then computed the impulse response of that biquad filter and truncated it to 1024 taps. I got the following filter:

To make the phase of my FIR filter linear, I applied an FFT on its impulse response, I put the imaginary part to 0 and came back to time domain by doing an inverse FFT.

The final taps I got can be found there.

Now, when I plot the frequency/phase responses of this filter using Numpy/Scipy or Matlab, I get exactly the same plot as the one of my expected filter (see figure 1). But, when I apply this filter on a pink noise signal by convolution and measure the transfer function, I get this instead:

I am pretty lost here, I would assume that convolving my noise signal with my filter and measuring the transfer function would give me the same curve as the one I computed theoritically but it seems not to be the case (that is, in my mind, the definition of what convolution does).

Some additional information:

1. I have verified that the problem does not come from my implementation of the convolution operation, by testing with many different implementations; if you still believe the problem come from there, you can try yourself with the taps I provided.


2. I trust the software I used to measure the transfer function depicted in the last plot, and anyway I could hear that the transfer function is not the one I expect.


3. 
   

   To plot the transfer functions from Python, I used the following code:

   
   

   import matplotlib.pyplot as plt
   
   import numpy as np
   from scipy import signal
   
   def plot_transfer_function(b, a=1):
       """Plot the frequency and phase response of a digital filter."""
       w, h = signal.freqz(b, a)
   
       fig = plt.figure()
       plt.title('Digital filter frequency response')
   
   
       # Frequency response
       ax1 = fig.add_subplot(211)
       plt.semilogx((w * SAMPLING_RATE) / (2 * np.pi), 20 * np.log10(abs(h)), 'b')
       # Make the plot pretty
       plt.ylabel('Amplitude [dB]', color='b')
       plt.xlabel('Frequency [Hz]')
       plt.axis([20, 20000, -12, 12])
       plt.grid()  # show the grid
   
       # Phase response
   
       ax2 = fig.add_subplot(212)
   
       angles = np.unwrap(np.angle(h))
       plt.semilogx((w * SAMPLING_RATE) / (2 * np.pi), np.degrees(angles), '#aa0000')
       # Make the plot pretty
       plt.ylabel('Angle (degrees)', color='#aa0000')
       plt.axis([20, 20000, -180, 180])
       plt.grid()
   
       fig.axes[0].set_axis_off()  # disable main figure axis
   
   
       plt.show()
   
   plot_transfer_function(TAPS)  # "TAPS" is an array containing the 1024 taps
   

So, what is wrong with my reasoning here? Why do I not measure the same frequency/phase response when calculating it and when measuring it after a convolution?

---

Additional question: for my FIR filter to have linear phase, I thought it should be symmetric (or anti-symmetric). Well, when I plot my 1024 taps, they appear to be pretty symmetric, excepted the first which does not have its equivalent at the end of the impulse response. Why is it so?