% конкретно та формула
P = qfunc(sqrt(((3*E)/N0) * (1/(q-1))));
Pe = 2*P - P^2;

% на всякий случай программа целиком
clear;
close all;

f0 = 2400;
T = 1/600;
q = 4;
i1 = [0, 0, 1, 1];
i2 = [0, 1, 0, 1];
A = 1;
X = zeros(0,q);
Y = zeros(0,q);
for i = 1:q
    X(i) = A*(1-(2*i1(i)/(sqrt(q)-1)));
    Y(i) = A*(1-(2*i2(i)/(sqrt(q)-1)));
end

Td=T/1000; time = 0:Td:T-0.0000001;
signal = zeros(q, 1000);
basis = zeros(2, 1000);
noise = zeros(1000);
for i = 1:q
    signal(i, :) = X(i)*((2/T)^(1/2))*cos(2*pi*f0*time)+Y(i)*((2/T)^(1/2))*sin(2*pi*f0*time);
end
basis(1, :) = (2/T)^(1/2)*cos(2*pi*f0*time);
basis(2, :) = (2/T)^(1/2)*sin(2*pi*f0*time);

SNR = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
gamma = zeros(length(SNR), 1);
errors = zeros(length(SNR), 1);
Pe = zeros(length(SNR), 1);
P = zeros(length(SNR), 1);
Pr = zeros(length(SNR), 1);
E = norm(signal(1,:))^2;

for snr = 1:length(SNR)
    gamma(snr) = 10^(SNR(snr)/10);
    N0 = E/gamma(snr);
    for i = 1:100000
        noise = sqrt(N0/2)*randn(1, 1000);
        I = randi([1, 4]);
        r = signal(I, :)+noise;
        mind = 9999999999;
        n1 = sum(r .* basis(1, :))*Td;
        n2 = sum(r .* basis(2, :))*Td;
        newI = 0;
        for j = 1:q
           if ((n1 - X(j))^2 + (n2 - Y(j))^2 < mind)
              mind = (n1 - X(j))^2 + (n2 - Y(j))^2;
              newI = j;
           end
        end
        if (newI ~= I)
            errors(snr) = errors(snr) + 1;
        end
    end
    Pr(snr) = errors(snr)/100000;
    P(snr) = qfunc(sqrt(((3*E)/N0) * (1/(q-1))));
    Pe(snr) = 2*P(snr) - P(snr)^2;
end

figure(1);
semilogy(SNR, Pe);
hold on;
for i = 1:length(SNR)
    plot(SNR(i), Pr(i), 'o', 'Color', 'k');
end
grid on;