samples_per_bit = Fs / bit_rate;
t = (0:num_samples - 1) / Fs;
deviation_range = linspace(5, 54, 50); 
ber_accumulated = zeros(length(deviation_range), 10); 
    data = round(rand(1, num_bits));
    for k = 1:length(deviation_range)
        Fc1 = Fc0 + deviation_range(k);
        fm_signal = zeros(1, num_samples);
            index_start = (i-1) * samples_per_bit + 1;
            index_end = i * samples_per_bit;
                fm_signal(index_start:index_end) = cos(2*pi*Fc0*t(index_start:index_end));
                fm_signal(index_start:index_end) = cos(2*pi*Fc1*t(index_start:index_end));
        Eb = (sum(abs(fm_signal).^2) / num_bits) / Fs;
        noise_power = sqrt(Eb/(10^(SNR/10))/2);
        noisy_signal = fm_signal + randn(size(fm_signal)) * noise_power;
            index_start = (i-1) * samples_per_bit + 1;
            index_end = i * samples_per_bit;
            sample = noisy_signal(index_start:index_end);
            ref0 = cos(2*pi*Fc0*t(index_start:index_end));
            ref1 = cos(2*pi*Fc1*t(index_start:index_end));
            corr0 = sum(sample .* ref0);
            corr1 = sum(sample .* ref1);
            received_bit = corr1 > corr0;
            errors = errors + (data(i) ~= received_bit);
        ber_accumulated(k, iter) = errors / num_bits;
ber_average = mean(ber_accumulated, 2);
ber_theory = berawgn(SNR,'fsk',2,'coherent') * ones(size(deviation_range));
plot(deviation_range, ber_average, 'b-o');
plot(deviation_range, ber_theory, 'r--');
xlabel('Frequency deviation (Hz)');
title('Dependence of BER on frequency separation');
legend('experimental BER', 'theoretic BER');