%% Mìøení spektrální masky pomocí dat ze spektrálního analyzátoru
% Skript slouží k meøení a vyhodnocování spektrální masky 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 promìných
clear all, close all, clc
disp('Inicializace promìných');
Tunered_freq = 88.2e6;    % Naladìná frekvence v Hz
Nsamp = 1900; % Poèet odebraných vzorkù. Maximální uvedená hodnota je 375000.
fsamp = Nsamp*500;     % Vzorkovací frekvence. Dovolené rozsahy: 225–300 kHz a 900–3200 kHz
Measure_Time =5*60; % Délka mìøení v sekundách
Frequency_Step = fsamp/Nsamp; %Vzdálenost jednotlivých bodù frekvenèního spektra navzorkovaného signálu
Sweep_step =10;% Vykreslení každého x-tého bodu ve spektru
addressMXA = '169.254.200.160'; % IP adresa pøipojeného pøístroje

%% Pøíprava matic pro ukládání dat:
data_IQBLock = nan(2*4*Nsamp,round(Measure_Time*(Nsamp/fsamp)*1.1));
IQsamples = nan(Nsamp,round(Measure_Time*(Nsamp/fsamp)*1.1));
%% Vytvorení masky
% Tato èást slouží k vytvoøení vektoru, který reprezentuje spektrální
% masku. Konstanta "Peak" urèuje pevnou úroveò masky pokud jsou mìøena data
% z generátoru.
disp('Vytvoøení masky')
AbsPosition = true; % "true" pro pevnou úroveò masky, "false" pro sledování maxima spektra
Peak = -20;
MaskSpan = 305e3;

MaskPoint = [-152.5e3/Frequency_Step -40;
    -124e3/Frequency_Step -30;
    -107.5e3/Frequency_Step -15;
    -74e3/Frequency_Step 0;
    74e3/Frequency_Step 0;
    107.5e3/Frequency_Step -15;
    124e3/Frequency_Step -30;
    152.5e3/Frequency_Step -40];

M1 = linspace(MaskPoint(1,2),MaskPoint(2,2),(MaskPoint(2,1)-MaskPoint(1,1)+1));
M2 =   linspace(MaskPoint(2,2),MaskPoint(3,2),(MaskPoint(3,1)-MaskPoint(2,1)+1));
M3 =   linspace(MaskPoint(3,2),MaskPoint(4,2),(MaskPoint(4,1)-MaskPoint(3,1)+1));
M4 =   linspace(MaskPoint(4,2),MaskPoint(5,2),(MaskPoint(5,1)-MaskPoint(4,1)+1));
M5 =   linspace(MaskPoint(5,2),MaskPoint(6,2),(MaskPoint(6,1)-MaskPoint(5,1)+1));
M6 =   linspace(MaskPoint(6,2),MaskPoint(7,2),(MaskPoint(7,1)-MaskPoint(6,1)+1));
M7 =   linspace(MaskPoint(7,2),MaskPoint(8,2),(MaskPoint(8,1)-MaskPoint(7,1)+1));
MaskVector = [M1(1:end-1) M2(1:end-1) M3(1:end-1) M4(1:end-1) M5(1:end-1) M6(1:end-1) M7(1:end)];
clear M1 M2 M3 M4 M5 M6 M7;

%% Vytvoøení RBW 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í RBW 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 = 5.3e3;
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 k=-N:N
    b(k+N+1)=3.011*(fco/SR)*exp(-pi*(a*k/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,950000 Hz,IMM,POS,0,1900']); % Nastavení parametrù mìøení: IF filtr 10 MHz, fsamp = 950 kHz Nsamp = 1900.
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

%Vytvoøení matice pro frekvenèní rozmítání
Vec_freq = (-152.5:0.5:152.5).*1e3;
Vec_t = (1:Nsamp);
Vector_ft = Vec_freq'*Vec_t;
Sweep_Matrix = (exp(((1i*2*pi).*Vector_ft)./fsamp))';

% Pøíprava matic pro zaznamenání vypoèteného spektra
Pmax = -inf(size(Sweep_Matrix,2),1);
Pmh = -inf(size(Sweep_Matrix,2),1);

% Pøíprava parametrù pro vykreslování hodnot
Freq_Axis = MaskSpan.*(-(MaskSpan/Frequency_Step)/2:(MaskSpan/Frequency_Step)/2)/(MaskSpan/Frequency_Step);
figure(2)
hold on
Handle1 = plot(Freq_Axis,Pmax,'-b');
Handle2 = plot(Freq_Axis,MaskVector','-g');
grid on
% Nastavení limitù a popisu osy x a y
ylim([-100 00])
xlim([-152.5e3 152.5e3])
xlabel('Frequence [Hz]','FontSize', 14,'FontWeight','Bold')
ylabel(['PSD [dBm/',sprintf('%.0f',Frequency_Step),' Hz]'],'FontSize', 14,'FontWeight','Bold')

% Vynulování pomocných ukazatelù
stepCount = 0;
stepCount2 = 0;

%% Mìøící cyklus

% Zaznamenání aktuálního èasu pro urèení délky mìøení
tStrat = clock;
tStop = clock;

disp('Zahájení mìøení')
while    Measure_Time>etime(tStop,tStrat) % cyklus získávání dat
    
    stepCount = 1+stepCount; %pocet cyklu
    stepCount2 = 1+stepCount2;
    %Pøíjem dat
    fprintf(signalAnalyzerObject,['TRAC:IQ:DATA?']);
    data_IQBLock(:,stepCount) = binblockread(signalAnalyzerObject);
    
    Vh = (stepCount2:Sweep_step:size(Sweep_Matrix,2)); % Index vybraných posunutých signálù
    if stepCount2 >= Sweep_step+1
        stepCount2 = 1;% Pøi vykreslení celého spektra, se obnový ukazatel pro opìtovné vykreslování spektra
    end
    
    % Dekódování pøijatých dat
    IQ_val_bin = reshape(de2bi(data_IQBLock(:,stepCount))',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(:,stepCount) = I_val + 1i*Q_val;
    
    IQsamples_sweep = (IQsamples(:,stepCount)).*Sweep_Matrix(:,Vh);% Simulace rozmítání
    
    FilterIQ = filter(b,1,IQsamples_sweep.*tukeywin(Nsamp,0.1));% RBW filtrace
    
    % Positive peak detektor
    FilterIQ= sqrt(real(FilterIQ).^2+imag(FilterIQ).^2);
    FilterIQmax = zeros(size(FilterIQ,1),size(FilterIQ,2));
    for k = 1: Nsamp/10
        FilterIQmax(10*(k-1)+1,:) = max(FilterIQ([10*(k-1)+1:10*k],:));
    end
    
    FilterIQ = 10*filter(b,1,FilterIQmax);% VBW filtrace
    FilterIQ = (FilterIQ.^2)./50; % Pøevedení na výkonové hodnoty
    
    % Vybrání zobrazené hodnoty spektra ze všech frekvenènì posunutých
    % signálù a pøevedení na dBm hodnoty
    Pmax(Vh) = (10*log10(max(FilterIQ))+30)';
    
    % Funkce Max Hold
    Pmh(Pmax>Pmh) = Pmax(Pmax>Pmh);
    Pmax = Pmh;
    
    % Stanovení pozice masky 
    if AbsPosition
        Position = Peak;
    else
        Position = max(Pmax);
    end
    
    % Vykreslení vypoètených hodnot spektra a masky
    try
        set(Handle1,'YData',Pmax)
        set(Handle2,'YData',MaskVector'+Position) % Pro pevnou hodnotu masky vložit promìnou "Peak".Pro sledování maxima vložit hodnotu "max(Pmax)"
        
    catch
        % Pøi zjištìní problému (napøíklad pøi zavøení zobrazovaného okna) je mìøení ukonèeno
        break;
    end
    
    % Zmìna barvy masky pøi jejím pøekroèení
    if any(Pmax > MaskVector'+Position)
        set(Handle2,'Color','r')
    else
        set(Handle2,'Color','g')
    end
    pause(0.0005)
    
    tStop = clock;
    % Zobrazení prùbìžného èasu mìøení
    etime(tStop,tStrat)
    
end

data_IQBLock(isnan(data_IQBLock)) = [];
IQsamples(isnan(IQsamples)) = [];

