clear
clc

% Загрузка исходного файла из ChipScope
fid = fopen('50.prn', 'r');
fgetl(fid);  % Пропускаем первую строку
y = fscanf(fid, '%g', [3 Inf]);
fclose(fid);
y = y(3, :) - mean(y(3, :));

% =========== Параметры программы ================
Thr = 0.25;            % Порог обнаружения сигнала
dF = 10000;            % Полоса вокруг спектрального отсчёта, относящаяся к нему

fd = 125e6;            % Частота дискретизации входного сигнала

% =========== Параметры программы ================

N = size(y, 2);

S = abs(fftshift(fft(y .* parzenwin(N)')));
S = S / max(S);
S = S(fix(N/2+1):N);

T=N/fd;
NB = ceil(dF*T)+5;     % Количество соседних спектральных составляющих

ind1 = find(S > Thr);  % Всё, что выше порога, считаем сигналом

a = zeros(size(S));    % Добавляем сюда соседние отсчёты - NB штук с каждой стороны
a(ind1) = ones(size(ind1));
b = zeros(size(S));
for i=-NB:NB
    b = b | circshift(a, [0 i]);
end

ind1 = find(b == 1);   % Отсчёты, соответствующие сигналу
ind2 = find(b == 0);   % Отсчёты, соответствующие шуму

if size(ind1, 2) > 6*NB
    fprintf('Результат сомнителен\n');
end

sinad = 10*log10(sum(S(ind1).^2) / sum(S(ind2).^2)) + 2.7;

N = size(y, 2);
f = ((fix(N/2):N-1)-N/2)/(N/2)*fd/2/1e6;
plot(f, 20*log10(S), 'r')
hold on
plot(f(ind2), 20*log10(S(ind2)), 'b-')
grid on
hold off
xlabel('f, MHz');
ylabel('S, dB');
set(gcf, 'color', 'white');

fprintf('SINAD: %g dB\n', sinad);

