[SOLVED] M2M4M6 moments based snr estimation

clc 
clear 
% Generate data. 
M = 16; % Alphabet size 
Pd = 2000; % Length of data 
k = log2(M); % Number of bits per symbol
x = randi([0,1],Pd,1); % Random binary bit stream
hBitToInt = comm.BitToInteger(k);
xsym = step(hBitToInt,x); % Convert the bits in x into k-bit symbols
qamsig = modulate(modem.qammod(M),xsym); % QAM signal 
snr_theory = 0:2:20; % Theoretical value of SNR 
QQ = 10.^(snr_theory/10); % SNR in decimal 
P1=1/4;
P2=3/4;
R(1)=1;
R(2)=3.15;
c2 = 1;
c4 = R(1).^4.*P1+R(2).^4.*P2;
c6 = R(1).^6.*P1+R(2).^6.*P2;
c8 = R(1).^8.*P1+R(2).^8.*P2;
c10 = R(1).^10.*P1+R(2).^10.*P2; 
for a=1:length(snr_theory) 
    SNR=num2str(snr_theory(a)); 
    for n=1:2000 
        rxsig = awgn(qamsig,snr_theory(a)); % Add Gaussion noise
        M2 = mean(rxsig.^2); % Second order moment 
        M4 = mean(rxsig.^4); % Fourth order moment 
        M6 = mean(rxsig.^6); % Sixth order moment
        m=(M2.*M4)./M6;
        v=[((m.*c6)-c4) ((9.*m.*c4)-c4-4) ((18.*m)-6) ((6.*m)-2)];
        r=roots(v);
        r(2);
        snr_est_M2M4M6(n) = 10.*log10(abs(r(2)));
        M2M4M6_snr_est(a) = mean(snr_est_M2M4M6); % estimation value 
        snr_est12(n) = 10.^(snr_est_M2M4M6(n)/10); % SNR in decimal 
        snr_est2(a) = mean(snr_est12);  
        NMSE_M2M4M6(a) = mean(((QQ(a)-snr_est12)/QQ(a)).^2); % Normalized MSE 
        STD_M2M4M6(a) = var(snr_est_M2M4M6); %standard deviation 
        NB_M2M4M6(a) = mean((QQ(a)-snr_est12)./QQ(a)); % normalized bias 
    end 
end 
figure(1) 
%plot estimation value 
plot(snr_theory,snr_theory); 
hold on 
plot(snr_theory,M2M4M6_snr_est,'-o'); 
hold off 
figure(2) 
% plot estimation NMSE 
semilogy(snr_theory,NMSE_M2M4M6,'-*') 
figure(3) 
% plot estimation STD 
plot(snr_theory,STD_M2M4M6,'-*') 
figure(4) 
% plot estimation NB 
plot(snr_theory,NB_M2M4M6,'-*') 
grid