This isn't the primary approach I would take to isolate phase equalization to a band of frequencies but to show the process to follow for the OP's approach that won't lead to the time domain aliasing that occurred I offer the following simple replacement of linear phase for regions outside of interest consistent with the average phase slope given by the delay of the filter:
Consider a phase response in the frequency domain as follows from an example 5000 bin DFT, together with a best fin linear line consistent with the average delay of the signal:
Below is the IFFT prior to any selection
imp = fft.ifft(np.exp(1j*phase))
plt.plot(fft.fftshift(imp))
We note that bin 0 represents "DC", bins 1 through 2500 the positive frequencies, and bins 2501 through 4999 represent the negative frequencies. (With bin 4999 the negative of bin 1, bin 4998 the negative of bin 2, etc).
If we consider this phase response to be the compensation we want, but only isolated to selective frequencies, we have to consider if the signal that we are compensating is real (after compensation, in which case the phase will be conjugate symmetric after compensation) in which case we need to include both positive and negative frequencies, or if the signal is complex in which case we need to decide if it is the positive or negative frequency we would like to compensate.
Showing both cases:
First for a real result after equalization:
f1 = 800
f2 = 1500
N = 5000
mask = np.zeros(N)
mask[f1:f2+1] = 1
mask[N-f1:N+1 -f2] = 1
result = np.where(mask, phase, plin )
With the following result for the impulse response:
And again for the case of selecting only the positive frequencies:
mask = np.zeros(N)
mask[f1:f2] = 1
result2 = np.where(mask, phase, plin )
I do not know the specifics of what the OP had done, but the time domain aliasing shown is consistent with zeroing every other bin in the frequency domain. For example for the phase plot shown in the beginning if we zeroed every other sample a similar time domain aliasing would occur:
phase2 = np.copy(phase)
phase2[::2]= 0
imp4 = fft.ifft(np.exp(1j*phase2))
plt.figure()
plt.plot(fft.fftshift(imp4))
