diff --git a/applications/lofar2/model/pfb_os/dsp.py b/applications/lofar2/model/pfb_os/dsp.py
index 19db600633d9292fa2b11dc12dde7dd515fb449d..c30de5e7f01d30f7cb67ac560d3e5c274d1e9ab8 100644
--- a/applications/lofar2/model/pfb_os/dsp.py
+++ b/applications/lofar2/model/pfb_os/dsp.py
@@ -121,33 +121,56 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
     return h, f, HF
 
 
-def fourier_interpolate(HF, Ncoefs):
+def fourier_interpolate(HFfilter, Ncoefs):
     """Use Fourier interpolation to create final FIR filter coefs.
 
     HF contains filter transfer function for N points, in order 0 to fs. The
     interpolation inserts Ncoefs - N zeros and then performs IFFT to get the
     interpolated impulse response.
+
+    Use phase correction depenent on interpolation factor Q to have fractional
+    time shift of hInterpolated, to make it symmetrical. Similar as done in
+    pfs_coeff_final.m and pfir_coeff.m. Use upper = conj(lower), because that
+    is easier than using upper from HFfilter.
     """
-    # . insert Ncoefs - N zeros between positive and negative frequencies
-    N = len(HF)
+    N = len(HFfilter)
     K = N // 2
+
+    # Interpolation factor Q can be fractional, for example:
+    # - Ncoefs = 1024 * 16
+    #   . firls: N = 1024 + 1  --> Q = 15.98439
+    #   . remez: N = 1024 --> Q = 16
+    Q = Ncoefs / N
+
+    # Apply phase correction (= time shift) to make the impulse response
+    # exactly symmetric
+    f = np.arange(N) / Ncoefs
+    p = np.exp(-1j * np.pi * (Q - 1) * f)
+    HF = HFfilter * p
+
     HFextended = np.zeros(Ncoefs, dtype=np.complex_)
     if is_even(N):
-        # Copy DC and K positive frequencies (including fs/2) in lower part
-        HFextended[0:K + 1] = HF[0:K + 1]
-        # Copy K - 1 negative frequencies in upper part
-        HFextended[-(K - 1):] = HF[K + 1:]
+        # Copy DC and K - 1 positive frequencies in lower part (K values)
+        HFextended[0:K] = HF[0:K]  # starts with DC at [0]
+        # Create the upper part of the spectrum from the K - 1 positive
+        # frequencies and fs/2 in lower part (K values)
+        HFflip = np.conjugate(np.flip(HF[0:K + 1]))  # contains DC at [0]
+        HFextended[-K:] = HFflip[0:K]  # omit DC at [K]
+        # K + K = N values, because N is even and K = N // 2
     else:
-        # Copy DC and K positive frequencies in lower part
-        HFextended[0:K + 1] = HF[0:K + 1]
-        # Copy K negative frequencies in upper part
-        HFextended[-K:] = HF[K + 1:]
-    hinterpolated = np.fft.ifft(HFextended)
-    if np.allclose(hinterpolated.imag, np.zeros(Ncoefs), rtol=c_rtol, atol=c_atol):
-        print('hinterpolated.imag ~= 0')
+        # Copy DC and K positive frequencies in lower part (K + 1 values)
+        HFextended[0:K + 1] = HF[0:K + 1]  # starts with DC at [0]
+        # Create the upper part of the spectrum from the K positive
+        # frequencies in the lower part (K values)
+        HFflip = np.conjugate(np.flip(HF[0:K + 1]))  # contains DC at [0]
+        HFextended[-K:] = HFflip[0:K]  # omit DC at [K]
+        # K + 1 + K = N values, because N is odd and K = N // 2
+    hInterpolated = np.fft.ifft(HFextended)
+    if np.allclose(hInterpolated.imag, np.zeros(Ncoefs), rtol=c_rtol, atol=c_atol):
+        print('hInterpolated.imag ~= 0')
     else:
-        print('WARNING: hinterpolated.imag != 0')
-    return hinterpolated.real
+        print('WARNING: hInterpolated.imag != 0 (max(abs) = %f)' % np.max(np.abs(hInterpolated.imag)))
+    return hInterpolated.real
 
 
 ###############################################################################