% % % % % % % % % % % % % % % % % % % % % % % % % %
% ===============================
% === FIND CALIBRATION VALUES ===
% ===============================
%
% Some general background follows:
%
% Eq(1) shows the general way power travels through the system. Power is
% 'created' at some level TX_POWER at the transmitter, passes through the TX
% antenna and is amplified with TX_ANT of gain, then passes through the
% channel and is attenuated by PL. The RX antenna then amplifies the signal
% by RX_ANT and this new power is attenuated by RX_ATTEN. The correlator
% system maps incoming power to outcoming power PDP.
%
% TX_POWER + TX_ANT <> PL <> RX_ANT + RX_ATTEN <==> SYSTEM <==> PDP : eq(1)
%
% Effectively, during measurements, we can say that four values affect the
% power coming out of the RX antenna: the antenna gains (TX and RX), the
% power of the TX, and the channel's path loss. Three of these values are
% known, and can be represented by the variable PCAL as in eq(2). The power
% that comes out of the antenna PR is related to PCAL by path loss PL as shown
% in eq(3).
%
% TX_POWER + TX_ANT + RX_ANT ==> PCAL : eq(2)
%
% PCAL - PL = OUT_OF_ANTENNA = PR : eq(3)
%
% PR then undergoes a linear transformation in the SYSTEM to create PDP
% (eq(4)). This SYSTEM transformation can be inverted, to find PR from
% PDP (eq(5)).
%
% SYSTEM(PR) = PDP : eq(4)
%
% PR = SYSTEM'(PDP) : eq(5)
%
% Thus, knowing TX_POWER, TX_ANT, RX_ANT, and RX_ATTEN, we can
% easily find PL for given PDP as in eq(6).
%
% PL = TX_POWER + TX_ANT + RX_ANT - SYSTEM'(PDP) - RX_ATTEN : eq(6)
%
% SYSTEM is a linear transformation defined as in eq(7). The inverse of
% SYSTEM, SYSTEM' follows from eq(7) and is defined in eq(8).
%
% SYSTEM( X ) = slope * X + intercept : eq(7)
% SYSTEM'( X ) = ( X - intercept ) / slope : eq(8)
%
% Of course, all above calculations need to be done in dB, as the
% relationship is only linear when using logarithmic scales.
%
% Eq(9) will, assuming all units are in dB, accurately anti-calibrate
% the system to return a value for PL given the inputs as stated below.
%
% PL =
% TX_POWER + TX_ANT + RX_ANT - RX_ATTEN - (PDP - intercept)/slope : eq(9)
%
% Finally: the returned value for PL is in REFERENCE to the path loss at
% 5m (or whatever your default distance is). This is due to the mapping
% method used below, where it is assumed that when RX_ATTEN = 0, total PL
% will be equal to zero. This reference is merely a translation shift,
% which is to say that it can be overcome by adding in PL_REF, which is
% simply the Friis PL value for 5m (etc.). Thus, true total path loss can
% be found by eqs(10) and above.
%
% PL_TRUE = PL + PL_REF
% = PL + 20*log10(4*pi*distance*frequency/3e8) : eq(10)
%
%
% Definitions:
% TX_POWER ~ the power going into the TX antenna
% TX_ANT ~ the gain of the TX antenna
% RX_ANT ~ the gain of the RX antenna
% RX_ATTEN ~ the setting of the RX attenuator
% PDP ~ the integrated value of the measured incoming waveform
% intercept ~ the intercept value returned by this program
% slope ~ the slope value returned by this program
% PL ~ the path loss of the channel
%
% % % % % % % % % % % % % % % % % % % % % % % % % %
clear
clc
%%%%%
% THE FOLLOWING MIGHT NEED TO BE MANUALLY SET
%%%%%
absoluteThreshold = 108; % in dBm/ns
attenSettings = 0:10:70; % in dB
distance = 5; % in metres
frequency = 28e9; % in hertz
isLinear = false; % true / false
relativeThreshold = 20; % in dB
RX_ANT = 24.5; % in dBi
TX_ANT = 24.5; % in dBi
TX_POWER = -8.55; % in dBm
%%%%%%%%%%%%
% EVERYTHING FROM HERE ON OUT IS FINE AS IS
%%%%%%%%%%%%
% Define matrices above so that they can be added to later.
cleanedPowerMatrix = [];
intRPowerMat = [];
offsets = [];
recvPowerMatrix = [];
timeMatrix = [];
truePowerMatrix = [];
% Default path of calibration PDP files
path = 'Folder Number 1\Calibration Data FS\';
% Some simple calculations to save time later
pathLoss = 20*log10(4*pi*distance / (3e8/frequency)); % Friis Path Loss
calPower = TX_POWER + TX_ANT + RX_ANT;
afterRX = calPower - pathLoss;
% In the loop below, load all calibration PDPs for given attenuation
% settings
for attenNum=1:length(attenSettings)
% Get path and file names
thisPath = [path,num2str((attenNum-1)*10),' dB\'];
fname = dir([thisPath,'PDP*']);
filePath = [thisPath,fname.name];
fileHandle = fopen(filePath,'r');
% Populate matrices with the time and volts^2 data from the text files
data = textscan(fileHandle,'%f,%f');
timeMatrix = [timeMatrix; transpose(data{1})];
recvPowerMatrix = [recvPowerMatrix; transpose(data{2})];
% Find what the actual power should be, given Friis PL equation
truePowerMatrix = [truePowerMatrix; afterRX - attenSettings(attenNum)];
end
% Data is converted to linear (from dB as loaded) to find thershold
recvPowerMatrix = 10.^(recvPowerMatrix./10);
% Because the data is assumed to be only one peak, we go 1600 sample points
% (100 ns after for our system) after the max, and assume everything from
% there on out to be noise. If noise peaks above this value, it either will
% not be by a significant amount or the measurement is wrong.
for i=1:length(attenSettings)
% Load one set of powers to work with
theseData = recvPowerMatrix(i,1:end);
% Find peak
subThres = max(theseData(find(theseData == max(theseData))+1600:60000));
% Clean based on this peak, do not allow anything less than 5/8 ns to pass
cleanedPowerMatrix = [cleanedPowerMatrix; removeNoiseSpike(theseData,subThres,10,1e-20)];
end
% Take integrals of power measured and populate a matrix
for j=1:length(attenSettings)
intRPowerMat = [intRPowerMat; sum(cleanedPowerMatrix(j,1:end))];
end
% Plot data and label for operator to decide linear range.
f = figure('visible','on');
plot(attenSettings,10*log10(intRPowerMat))
title('Calibration mapping curve','fontsize',20)
xlabel('Attenuator setting (dB)','fontsize',17)
ylabel('Unscaled power received (dBm)','fontsize',17)
% Pseudo do-while loop to run until user agrees with values entered for
% linear range.
isGood = 0;
while(isGood ~= 1)
disp('Please enter extrema of linear range.')
min = input('Min (dB): ');
max = input('Max (dB): ');
% Convert value in dB to array index (works for our system)
min = round(min/10)+1;
max = round(max/10)+1;
% Find a linear transformation values for selected linear range
fitValues = polyfit(attenSettings(min:max),10*log10(transpose(intRPowerMat(min:max))),1);
slope = fitValues(1);
intercept = fitValues(2);
% Create [a,b] and [c,d] for {x} and {y} plot of selected
% line-of-best-fit
a = attenSettings(1);
b = attenSettings(end);
c = slope*a + intercept;
d = slope*b + intercept;
% Replot and user-verification
plot(attenSettings,10*log10(intRPowerMat),'b',[a,b],[c,d],'g')
title('Calibration mapping curve','fontsize',20)
xlabel('Attenuator setting (dB)','fontsize',17)
ylabel('Unscaled power received (dBm)','fontsize',17)
isGood = input('If new line fits, enter 1: ');
end
% Kill the graph to save memory etc.
set(f,'visible','off')
% The following section offers a statistical analysis of the points chosen
% by the user in addition to a sanity check on the calibration process
disp('This section is meant to provide numerical confirmation that')
disp('the numbers entered are correct.')
disp('For the linear range that you entered, the path')
disp('loss(es) calculated follow(s):')
% Show the user how much of an offset there is between data in linear range
% and the linear fit generated, both in dB and linear
for i=min:max
offset = (10*log10(intRPowerMat(i))-intercept) / slope - attenSettings(i);
disp([num2str(attenSettings(i)),' dB ~ ',num2str(offset)]);
offsets = [offsets,offset];
end
disp('The above numbers are in dB and should be close')
disp('to zero, whereas the following number (ratio) is')
disp(['in linear and should be close to 1: ',num2str(mean(10.^(offsets/10)))])
disp(' ')
disp('If the above values are not satisfactory, rerun the program!')
% Subtract out calPower because TX_POWER, TX_ANT, and RX_ANT potentially
% change from one measurement to the next. We can account for these changes
% in processing when we know the specifics.
intercept = intercept - calPower;
% Save relevant variables
save('calVars2','slope','intercept')
% Housekeeping
clear all