How to remove Phase Offset from PSK signal?

Question

Suppose you have captured IQ data from a PSK modulated signal with intermittent blanks between signal transmission as shown in the figure below.

I want to be able to automate and run an algorithm on a multitude of data captures like this to determine the phase offset and then plot the Phase corrected like the figure below in Matlab. What's the best way to do this?

clc;clear;close all;

rSymbolDuration = 1e-3; 
rSampleRate     = 1e6;
rPhaseOffset    = 10;

a = repelem([1 0 -1 0],round(rSampleRate*rSymbolDuration));
b = repelem([0 1 0 -1],round(rSampleRate*rSymbolDuration));

caWaveform = [];

while length(caWaveform)<100e3
    caWaveform = [caWaveform ,...
        [awgn(repmat(complex(a,b),1,randi(3,1,1)),35),...
        awgn(repmat(complex(ones(rSampleRate*2*rSymbolDuration,1)',...
        ones(rSampleRate*2*rSymbolDuration,1)'),...
        1,randi(3,1,1)),35)+(-1-1i)]];
end

xData = (1:length(caWaveform))./rSampleRate;
caWaveformR = caWaveform .* exp(1j.*2.*pi.*(-rPhaseOffset.*xData));

figure;plot(xData,real(caWaveformR));hold on;plot(xData,imag(caWaveformR));grid on;xlabel('time(secs)');ylabel('\pi');title('IQ vs Time')
figure;plot(xData,angle(caWaveformR));grid on;xlabel('time(secs)');ylabel('Phase(rads)');title('Phase vs Time')
figure;plot(real(caWaveformR),imag(caWaveformR),'.');grid on;xlabel('I');ylabel('Q');title('IQ Scatter Plot')
figure;plot(xData,angle(caWaveform));grid on;xlabel('time(secs)');ylabel('Phase(rads)');title('Phase vs Time')

Answer 1

The blanks are easy to detect in this case, given simultaneous I and Q values below threshold, with a threshold based on the overall signal standard deviation of the signal and the noise with no signal present. Have a running determination of both metrics and set a threshold half way in between. When either I or Q are above threshold, use the resulting I+jQ toward the running metric for signal standard deviation and compute the angle using atan2(Q, I). When the signal is below the threshold, use the resulting I+jQ toward the running metric for noise standard deviation (in both cases this is a sum of squares divided by total number of samples in the buffer used to track the metrics).