%% Mìøení spektrální masky pomocí dat ze spektrálního analyzátoru
% Skript slouží k meøení a vyhodnocování hodnota frekvenèního zdvihu a
% modulaèního výkonu signálu pomocí dat ze spektrálního analyzátoru.
% Volitelné parametry pro mìøení signálu lze nastavit v sekci pro inicializaci promìných
%% Inicializace promenych
clear all, clc
Tunered_freq = 88.2e6;   %   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 = 15*60;
Frequency_Step = fsamp/Nsamp;
addressMXA = '169.254.200.160'; % IP adresa pøipojeného pøístroje

%% Pøíprava matic pro ukládání dat:
DemodData =nan(Nsamp,round(Measure_Time*Frequency_Step*1.1));
fn_max = nan(1,round(Measure_Time*Frequency_Step*1.1));
IQsamples = nan(Nsamp,round(Measure_Time*Frequency_Step*1.1));
Pm = nan(1,round(Measure_Time*Frequency_Step*1.1));
TimeMarker = nan(ceil(Measure_Time/60)+1);
%% Inicializace objektu
disp('Inicializace objektù');

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

%% IF filter
disp('Vytváø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.

%  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 TimeMarker_index=-N:N
    b(TimeMarker_index+N+1)=3.011*(fco/SR)*exp(-pi*(a*TimeMarker_index/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);


%% Navázání spojení s analyzátorem
disp('Pøipojování')
signalAnalyzerObject = tcpip(addressMXA, 5025);
signalAnalyzerObject.InputBufferSize = 5e6;
signalAnalyzerObject.Timeout = 10;
fopen(signalAnalyzerObject);

fprintf(signalAnalyzerObject, '*RST'); % Reset nastavení pøístroje
instrumentInfo = query(signalAnalyzerObject, '*IDN?'); % Vyžádání identifikaci pøístroje
disp(['Instrument identification information: ' instrumentInfo]);

% Nastavení pøijmu IQ dat
fprintf(signalAnalyzerObject,[':SENSe:FREQuency:CENTer ' num2str(Tunered_freq)]); % Nastavení naladìné frekvence
fprintf(signalAnalyzerObject,['TRAC:IQ:STAT ON']);
fprintf(signalAnalyzerObject,['TRAC:IQ:SET NORM,10 MHz,1000000 Hz,IMM,POS,0,50000']); % Nastavení parametrù mìøení: IF filtr 10 MHz, fsamp = 1 MHz Nsamp = 50000.
fprintf(signalAnalyzerObject,['FORMat REAL,32']); % Formát pøijatých dat
fprintf(signalAnalyzerObject,['TRACe:IQ:DATA:FORMat IQBLock']); % Poøadí pøijatých IQ dat

%% Pøíprava mìøícího cyklu

Cutoff =15;
P19kHz = ((19e3/kFM)^2)/2;
Vec_fn  = 0:1e3:100e3;
SumCon = zeros(1,length(Vec_fn));
TimeMarker_index = 0;
stepCount = 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 [dB]','FontSize', 14,'FontWeight','Bold')
xlim([1 length(Pm)]);

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

% ziskavani IQ dat
disp('Získávání IQ složek')
while    Measure_Time>etime(tStop,tStrat)
    
    stepCount = 1+stepCount; %pocet cyklu
    fprintf(signalAnalyzerObject,['TRAC:IQ:DATA?']);
    data_IQBLock = binblockread(signalAnalyzerObject);
    
    IQ_val_bin = reshape(de2bi(data_IQBLock)',32,2*Nsamp)';
    Vec = repmat((-23:1:-1),length(IQ_val_bin),1);
    Man = 1+sum((IQ_val_bin(:,1:23).*2.^Vec),2);
    Exponent = bi2de(IQ_val_bin(:,24:31))-127;
    Sgn = (-1).^IQ_val_bin(:,32);
    IQ_val = Sgn.* Man.* 2.^Exponent;
    I_val =  IQ_val(1:length(IQ_val)/2,:);
    Q_val =  IQ_val(length(IQ_val)/2+1:end,:);
    IQsamples = I_val + 1i*Q_val;
    
    
    FiltredData = filter(b,1,IQsamples);
    DemodData = pFMDemod(FiltredData);
    DemodData([Cutoff:end]) = DemodData([Cutoff:end])-mean(DemodData([Cutoff:end]));
    
    Pp = (1/(Nsamp-Cutoff))*sum(abs(DemodData([Cutoff:end])).^2); 
    Pm (:,stepCount) = 10*log10(Pp/P19kHz); % Výpoèet modulaèního výkonu
    
    Umax = max(DemodData([Cutoff:end]));
    fn_max (:,stepCount) = kFM.*Umax; % Maximální zjištìný frekvenèní zdvih
    
    fn((stepCount-1)*length(DemodData([Cutoff:end]))+1:(stepCount)*length(DemodData([Cutoff:end]))) =abs(DemodData([Cutoff:end]))*kFM;
    fn_con = ceil(fn((stepCount-1)*length(DemodData([Cutoff:end]))+1:(stepCount)*length(DemodData([Cutoff:end]))));
    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))= [];


% uvloneni objektu
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)
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(4)
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')
