Understanding Phase-offset in 4-QAM/QPSK

function [audio_recorder] = receiver(fc)
    fs = 22050;
    % https://se.mathworks.com/help/matlab/ref/audiorecorder.html
    % audiorecorder creates and returns an audiorecorder object with these example properties:
    % SampleRate: 8000
    % BitsPerSample: 8 | can either be 8, 16 or 24
    % NumChannels: 1 | To obtain the ID of a device, use audiodevinfo -> https://se.mathworks.com/help/matlab/ref/audiodevinfo.html
    audio_recorder = audiorecorder(fs, 24, 1);

    audio_recorder.UserData.fs            = fs;
    audio_recorder.UserData.fc            = fc;
    audio_recorder.UserData.stopBit       = 0;

    % Attach callback function
    % How often the function should be called in seconds
    time_value = 0.5;

    % Attach a function that should be called every timer_value seconds.
    set(audio_recorder, 'TimerPeriod', time_value, 'TimerFcn', @audioTimerFcn)

    % Add user data for callback function
    audio_recorder.UserData.receive_complete    = 0;    % this is a flag that the while loop in the GUI will check
    audio_recorder.UserData.pack                = [];   % allocate for data package
    audio_recorder.UserData.pwr_spect           = [];   % allocate for PSD
    audio_recorder.UserData.const               = [];   % allocate for constellation
    audio_recorder.UserData.eyed                = [];   % allocate for eye diagram

    % start recording
    record(audio_recorder);

end

% CALLBACK FUNCTION
% This function will be called every [time_value] seconds, where time_value
% is specified above. Note that, as stated in the project MEMO, all the
% fields: pwr_spect, eyed, const and pack need to be assigned if you want
% to get outputs in the GUI.
function audioTimerFcn(recObj, event, handles)

    [~, tsamp, const, bps, fsfd, Rs, Ts, preamble, rrc, stopBit, NoBits, ~] = common();
    fc = recObj.UserData.fc;
    fs = recObj.UserData.fs;
    
    disp('Callback triggered')
    disp(' ')


    % implement an RRC pulse shape
    pulse_train = rtrcpuls(rrc.roll_off_factor, Ts, fs, rrc.span);
    
    % disp(length(pulse_train));
    % captured_signal = getaudiodata(recObj)';
    % Generate a random delay between 10 and 500
    % Check Ex8.m file on Canvas
    % delay = zeros(1, randi([10, 20], 1, 1));
    % delay_data = awgn(delay, 10, 'measured');

    cs = getaudiodata(recObj)';

    % received signal with delay
    % cs = [delay_data, info];

    % Remove the stop bit
    cs = cs(recObj.UserData.stopBit + 1 : end);

    % Normalize Rx signal
    cs = cs./max(abs(cs));

    % Demodulation, Put the baseband signal on carrier signal
    demod_rx = sqrt(2) * cs.*exp(1j * 2 * pi * fc * (0:length(cs) - 1) * tsamp);       

    % Timer Synchronization
    upsampled_preamble = upsample(preamble, fsfd);

    % shape the known preamble like the rx preamble
    preamb_train = conv(pulse_train, upsampled_preamble);

    % correlate the processed preamble and demodulated signal to find the similarity
    corr = conv(demod_rx, fliplr(preamb_train));         

    % noise threshold
    noiseThreshold = 10;                                                         

    % find location of max correlation
    [maxCorr, Tmax] = max(corr);                                               

    % find delay (it will take length preamble)
    Tx_hat = Tmax - length(preamb_train);                                       
    
    % pass the signal only if
    %  1. the correlation is not with noise
    %  2. the delay is positive
    %  3. the whole signal is captured
    a = noiseThreshold > abs(maxCorr);
    b = Tmax < length(preamb_train);
    d = length(demod_rx) - Tmax;
    e = fsfd * (NoBits / 2);
    c = d < e;

    % disp([a, b, c]);

    if(a || b || c)
        fprintf('No signal found \n')
        return;
    end
    
    fprintf("Delay occured is %s samples\n", num2str(Tx_hat));

    % Take note of the stop bit
    recObj.UserData.stopBit = length(cs);

    % ===================================================================
    %       Matched filter calculations
    % ===================================================================

    % Create a matched filter impulse response from the pulse_train
    matched_filter = fliplr(conj(pulse_train));

    % Now we run collected signal through the filter
    matched_filter_output = filter(matched_filter, 1, demod_rx);

    % Remove the transient response in the FIR bandpass filter
    % Transient response is the time the system takes from
    % an equilibrium position to another steady state
    matched_filter_output = matched_filter_output(length(matched_filter) : end);

    % ===================================================================
    %       Phase Synchronization
    % ===================================================================   
    % Match the received preamble symbols with what we have so far
    rx_preamble = matched_filter_output(Tx_hat + 1 : fsfd : Tx_hat + length(upsampled_preamble));

    % Normalize the preamble symbols to 1 unit
    rx_preamble = rx_preamble./mean(abs(rx_preamble));

    % At this stage, we can already extract the phase difference
    % between the known preambles and the received ones.
    phase = angle(rx_preamble) - angle(preamble);
    phase = mod(phase, 2 * pi);

    if(var(phase) > 5)
        error = 5;
        phase(phase < (2 * pi + error) & phase > (2 * pi - error)) = phase(phase < (2 * pi + error) & phase > (2 * pi - error)) - 2 * pi;
    end

    avgPhase = mean(phase);

    % Get the real signal received
    matched_filter_output = matched_filter_output(Tx_hat+  length(upsampled_preamble) + 1 : Tx_hat + length(upsampled_preamble) + (NoBits / 2) * fsfd);

    % Get the correct phase used at the Tx side
    matched_filter_output = matched_filter_output * exp(-1j * avgPhase);

    fprintf("The signal is phase shifted by %s deg CCW\n", num2str(avgPhase * 180/pi));

    % ===================================================================
    %      Extract symbols from the signal
    % =================================================================== 
    % Sample the matched filter signal at fs
    rx = matched_filter_output(1:fsfd:end);

    % Normalize the vector
    rx = rx./mean(abs(rx));


    % ===================================================================
    %      Turn the symbols to bits
    % ===================================================================
    % Minimum Eucledian distance detector
    % Relate the detection to Detection region

    % Compute the distance to each received symbol from where it should be
    metric = abs(repmat(rx.', 1, 4) - repmat(const, length(rx), 1)).^2;   

    % find the closest for each received symbol
    [tmp, msg_hat] = min(metric, [], 2);     

    % Get the index of the symbol in the constellation
    msg_hat = msg_hat' - 1;             

    % Convert symbols to bits
    msg_hat = de2bi(msg_hat, 2, 'left-msb')';                                   
    data_hat = msg_hat(:)';
    
    bits = data_hat(1 : NoBits);

    x = rx(1 : (NoBits / 2));
    pulse_train = matched_filter_output;
   
end

function [const, fs, Tsamp, bitperSymb, alpha, fsfd, Rs, Ts, preamb, span] = common()

    const = [(1 + 1i), (1 - 1i), (-1 -1i), (-1 + 1i)]/sqrt(2);
    fs = 22050;
    Tsamp = 1 / fs;
    M = length(const);                                                          
    bitperSymb = log2(M);                                                       
    alpha = 0.4;                                                                
    fsfd = 120;                                                                 
    Rs = fs/fsfd;                                                               
    Ts = 1/Rs;    
    preamb = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];    
    span = 6;

end