%% Mìøení spektrální masky pomocí RTL-SDR pøijímaèe
% Skript slouží k meøení a vyhodnocování hodnota frekvenèního zdvihu a
% modulaèního výkonu signálu pomocí pøipojeného RTL-SDR zaøízení.
% Volitelné parametry pro mìøení signálu lze nastavit v sekci pro inicializaci promìných
%% Inicializace promìných
clear all, close all, clc
disp('Inicializace promìných');
% Nastavení parametry pro mìøení:
Freq_Offset = 180e3;
Tunered_freq = 86e6-Freq_Offset;       % Naladìná frekvence v Hz
fsamp = 1e6;           % Vzorkovací frekvence Dovolené rozsahy: 225–300 kHz a 900–3200 kHz
Nsamp = 50000;        % Poèet odebraných vzorkù. Maximální uvedená hodnota je 375000
Measure_Time =1*60; % Délka mìøení v sekundách
Freq_step = fsamp/Nsamp; %Vzdálenost jednotlivých bodù frekvenèního spektra navzorkovaného signálu
RadioAddress = '0'; % Urèení portu pro RTL-SDR zaøízení
Datatype = 'single';  % Typ pøijatých dat
ppm =42;              % Hodnota sloužící k frekvenèní korekci nepøesnosi RTL-SDR zaøízení
Gain =0;             % Nastavené zesílení

%% Pøíprava matic pro ukládání dat:
fn_max = nan(1,round(Measure_Time*Freq_step*1.1));
Pm = nan(1,round(Measure_Time*Freq_step*1.1));
TimeMarker = nan(ceil(Measure_Time/60));
fn = nan(ceil(1.1*Nsamp*Measure_Time*Freq_step),1);
lost = nan(1,round(Measure_Time*Freq_step*1.1));
latency = nan(1,round(Measure_Time*Freq_step*1.1));
%% Inicializace objektù
disp('Inicializace objektù');
% Inicializace objektu comm.SDRRTLReceiver sloužícího k pøíjmo dat z RTL
obj_rtlsdr = comm.SDRRTLReceiver(...
    RadioAddress,...
    'CenterFrequency',Tunered_freq,...
    'EnableTunerAGC' , false, ...
    'SampleRate',fsamp,...
    'SamplesPerFrame',Nsamp,...
    'OutputDataType',Datatype,...
    'FrequencyCorrection',ppm,...
    'EnableBurstMode',true,...
    'NumFramesInBurst',10,...
    'TunerGain',Gain);

%  Inicializace objektu comm.FMDemodulator sloužícího k demodulaci dat
kFM = 75e3; % kFM
pFMDemod = comm.FMDemodulator('SampleRate',fsamp,...
    'FrequencyDeviation',kFM);


%% Vytvoøení IF filtru

% Tato èást byla získána z webových stránek spoleènosti Mathwork. Konkrétnì
% z: https://www.mathworks.com/matlabcentral/fileexchange/12606-1d-gaussian-lowpass-filter
% Autor uvedl popis funkce v anglickém jazyce, proto je tento popis
% ponechán. Byla zmìnìna hodnota konstanty "N" z N=ceil(0.398*SR/fco)
% na N=1.25*ceil(0.398*SR/fco); pro zvýšení potlaèení v nepropustném pásmu
% filtru.

disp('Vytváøení IF filtru');
%  Return coefficients of Gaussian lowpass filter.
% SR=sampling rate, fco=cutoff (-3dB) freq, both in Hz.
% Coeffs for FIR filter of length L (L always odd) are computed.
% This symmetric FIR filter of length L=2N+1 has delay N/SR seconds.
% Examples of use
%    Compute Gaussian filter frequency response for SR=1000, fco=50 Hz:
%    freqz(gaussfiltcoef(1000,50),1,256,1000);
%    Filter signal X sampled at 5kHz with Gaussian filter with fco=500:
%    y=filter(gaussfiltcoef(5000,500),1,X);
% SR, fco are not sanity-checked.  WCR 2006-10-11.
b=0;
fco = 250e3/2;
SR = fsamp;
a=3.011*fco;
% N=ceil(0.398*SR/fco);   %filter half-width, excluding midpoint
%Width N corresponds to at least +-3 sigma which captures at least 99.75%
%of area under Normal density function. sigma=1/(a*sqrt(2pi)).
N=1.25*ceil(0.398*SR/fco); 
L=2*N+1;                %full length of FIR filter
for Vec_fn=-N:N
    b(Vec_fn+N+1)=3.011*(fco/SR)*exp(-pi*(a*Vec_fn/SR)^2);
end;
%b(k) coeffs computed above will add to almost exactly unity, but not
%quite exact due to finite sampling and truncation at +- 3 sigma.
%Next line adjusts to make coeffs b(k) sum to exactly unity.
b=b/sum(b);
%% Pøíprava mìøícího cyklu
% Kontrola pøipojení RTL-SDR pøijímaèe
disp('Provádím kontrolu pøipojení RTL');
if isempty(sdrinfo(obj_rtlsdr.RadioAddress)) % kontrola pøipojení RTL
    error('Komunikace s RTL selhala');
end

% Zpìtné posunutí frekvenèního spektra signálu
Vec_freq =180e3;
Vec_t = (1:Nsamp);
Vector_ft = Vec_freq'*Vec_t;
FC_matrix = (exp(((1i*2*pi).*Vector_ft)./fsamp))';


Vec_fn  = 0:1e3:100e3;
SumCon = zeros(1,length(Vec_fn));
Cutoff =15;
TimeMarker_index = 0;

% Pøíprava parametrù pro vykreslování hodnot
figure (1)
Handle =  plot(Pm);
grid on
grid minor
xlabel('n[-]','FontSize', 14,'FontWeight','Bold')
ylabel('Pm [dBr]','FontSize', 14,'FontWeight','Bold')
xlim([1 length(Pm)]);

disp('Pøipravuji sbìr IQ dat');
IQsamples = step(obj_rtlsdr);
DemodData = pFMDemod(IQsamples);


stepCount = 0;
k_Pm = ((19e3/kFM)^2)/2;

% Nastaveni kontroly casu
tStrat = clock;
tStop = clock;

% ziskavani IQ dat

while    Measure_Time>etime(tStop,tStrat) % cyklus získávání dat
    stepCount = 1+stepCount; %pocet cyklu
    [IQsamples, ~, lost(:,stepCount), latency(:,stepCount)] = step(obj_rtlsdr);% naètení dat z RTL
    IQsamples_sweep = (IQsamples).*FC_matrix;
    FiltredData = filter(b,1,IQsamples_sweep); 
    DemodData = pFMDemod(FiltredData);
    DemodDataCut = DemodData([Cutoff:end])-mean(DemodData([Cutoff:end]));
    
    Pp = (1/(Nsamp-Cutoff))*sum(abs(DemodDataCut).^2);
    Pm (:,stepCount) = 10*log10(Pp/k_Pm);% Výpoèet modulaèního výkonu
    
    Umax = max(DemodDataCut);
    fn_max (:,stepCount) = kFM.*Umax; % Maximální zjištìný frekvenèní zdvih
    
    fn((stepCount-1)*length(DemodDataCut)+1:(stepCount)*length(DemodDataCut)) =abs(DemodDataCut)*kFM;
    fn_con = ceil(fn((stepCount-1)*length(DemodDataCut)+1:(stepCount)*length(DemodDataCut)));
    
    %Zvìtšení matice pro hodnoty fn pøi pøekroèení maxima
    if max(fn_con) >max(Vec_fn)
        Vec_fn = [Vec_fn max(Vec_fn)+1e3:1e3:ceil(max(fn_con/1e3))*1e3];
        SumCon = [SumCon  (zeros(length(length(SumCon)+1:length(Vec_fn)),1))'];
        SumNorm = SumCon/max(SumCon);
    end
    
    Sumfn_con =sum((fn_con>Vec_fn));
    SumCon = SumCon+Sumfn_con; % Akumulace hodnot fn
    SumNorm = SumCon/max(SumCon); % Normování hodnot
    
    
    
    try
        set(Handle,'YData',Pm)
        pause(0.0005)
    catch
        break;
    end
    
    tStop = clock;
    etime(tStop,tStrat)
    % Zaznamenání indexu vzorku po úsecích 60 sekund
    if  etime(tStop,tStrat)>=60*TimeMarker_index
        TimeMarker_index = TimeMarker_index +1;
        TimeMarker(TimeMarker_index) =  stepCount;
        
    end
end
% Odrtranìní nevyužitých prvkù matic
fn(isnan(fn))= [];
TimeMarker(isnan(TimeMarker))= [];
fn_max(isnan(fn_max))= [];
Pm(isnan(Pm))= [];
lost(isnan(lost))= [];
latency(isnan(latency))= [];

% Uvolnìní objektù
release(obj_rtlsdr)
release(pFMDemod)

%% Vyhodnocování výsledkù mìøení
disp('Vyhodnocování výsledkù mìøení')
figure(2)
semilogy(SumNorm);
xlim([1 length(SumCon)]);
grid on
grid minor
xlabel('f [kHz]','FontSize', 14,'FontWeight','Bold')
ylabel('Množství hodnot fn>f','FontSize', 14,'FontWeight','Bold')
figure(3)
plot(lost)
xlabel('n[-] [min]','FontSize', 14,'FontWeight','Bold')
ylabel('Ztracených vzorkù','FontSize', 14,'FontWeight','Bold')
figure(4)
stem(fn_max);
ylabel('Maximální frekvence [Hz]','FontSize', 14,'FontWeight','Bold')
xlabel('n[-]','FontSize', 14,'FontWeight','Bold')
xlim([1 length(fn_max)]);
grid on
grid minor
Pm_60s =0
for c = 1:length(TimeMarker)-1
    Pm_60s(c) = mean(Pm(TimeMarker(c):TimeMarker(c+1)));
end
figure(5)
plot(Pm_60s)
grid on
grid minor
xlabel('t [min]','FontSize', 14,'FontWeight','Bold')
ylabel('Pm [dB]','FontSize', 14,'FontWeight','Bold')

disp('Vyhodnocování ukonèeno')
