Jens,
I looked at your code and
I don’t know if this makes a difference,
but I think the line
data = d_resamp(256:1792)
should be substituted by
data = d_resamp(257:1793)
since in an N=2m point FFT the DC is the m+1-st component
after fftshift (in our case it is the 1025-th component).
I still believe that ftaer correct frequency compensation the
picture should look like an unrotated QPSK, so my bet
is that there is something wrong with the frequency correction…
Achilleas