Проблемы с демодуляцией QAM-16 сигнала

В момент демодуляции создается последовательность из 4 значений, хотя коэффициент модуляции и демодуляции равен 16. В чем может быть проблема и как это исправить?

clc
close all
clear all
M = 16;
L=10000;
T=10;
W=43; %SNR

%% Передатчик
bits1 = randi([0 M-1], 1, L);
bits2 = randi([0 3], 1, L/T);

mod1 = qammod(bits1, M,'UnitAveragePower', true);
% scatterplot(mod1)

mod2 = qammod(bits2, 4, 'UnitAveragePower',true);
% scatterplot(mod2)

for i = 1:length(bits2)
    sym(i,:)=[mod1((i-1)*T+1:i*T) mod2(i)];
end
f=0;
for k=1:L/T
   for j=1:T+1
       f=f+1;
       symbols(f)=sym(k,j);
   end
end
% scatterplot(symbols)
%% Канал передачи
% symbols = awgn(symbols,W,"measured");

pnoise = comm.PhaseNoise('Level',[-70 -104 -110],'FrequencyOffset',[1e4 1e5 2e5], 'SampleRate', 28e6);
symbols2 = pnoise([zeros(1e5,1);symbols']);
% symbols2 = pnoise(symbols);

fc = 1e6; % Carrier frequency in Hz
fs = 28e6; % Sample rate in Hz.
phNzLevel = [-70 -104 -110]; % in dBc/Hz
phNzFreqOff = [1e4 1e5 2e5]; % in Hz
Nspf = 6e6; % Number of Samples per frame
freqSpan = 400e3; % in Hz, for spectrum computation

sinewave = dsp.SineWave( ...
    Amplitude=1, ...
    Frequency=fc, ...
    SampleRate=fs, ...
    SamplesPerFrame=Nspf, ...
    ComplexOutput=true);
pnoise = comm.PhaseNoise( ...
    Level=phNzLevel, ...
    FrequencyOffset=phNzFreqOff, ...
    SampleRate=fs);
sascopeRBW100 = spectrumAnalyzer( ...
    SampleRate=fs, ...
    Method="welch", ...
    FrequencySpan="Span and center frequency", ...
    CenterFrequency=fc, ...
    Span=freqSpan, ...
    RBWSource="Property", ...
    RBW=100, ...
    SpectrumType="Power density", ...
    SpectralAverages=10, ...
    SpectrumUnits="dBW", ...
    YLimits=[-150 10], ...
    Title="Resolution Bandwidth 100 Hz", ...
    ChannelNames={'signal','signal with phase noise'}, ...
    Position=[79 147 605 374]);

x = sinewave();
y = pnoise(x);

sascopeRBW100(x,y)

symbols2 = symbols2(1e5+1:end);
%% Приёмник
prim(1,:)=symbols2(1:T);
for i = 1:length(bits2)-1
    prim(i+1,:)=symbols2(i*T+i+1:i*T+i+T);
end

for b = 1:length(mod2)
    qam4(b)=symbols2((T+1)*b);
end

h=0;
for u=1:L/T
   for l=1:T
       h=h+1;
       priem(h)=prim(u,l);
   end
end

for g = 1:L/T
    phase_error(g) = angle(qam4(g) / mod2(g));
    compensated4(g) = qam4(g) .* exp(-1i * phase_error(g));
end
for v = 1:L/T
    phase_errorM((v-1)*T+1:v*T)=phase_error(v);
end
for f = 1:L
    compensatedM(f) = priem(f) .* exp(-1i*phase_errorM(f));
end


demod=qamdemod(compensatedM, M, 'bin','OutputType','bit');
% scatterplot(qam4)
% scatterplot(symbols2)
% scatterplot(compensatedM)

[number,ratio]=biterr(bits1,demod);
evm = lteEVM(demod,bits1);

figure('Position',[200 200 1080 540])
subplot(1,2,1)
scatter(real(symbols2),imag(symbols2),300,".")

subplot(1,2,2)
scatter(real(compensatedM),imag(compensatedM),300,".")