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');