Description of analyzePassbandMeasurement

0001 function analyzePassbandMeasurement(varargin)
0002 % analyzePassbandMeasurement(dirname,prefix)
0003 % OR
0004 % analyzePassbandMeasurement(dirname,prefix,timeZone)
0005 %
0006 % This function analyzes the data from passband measurements. It reads in a
0007 % standard set of files from the directory given in dirname. It also reads
0008 % in the appropriate data from the archive. The prefix to the standard file
0009 % names can be blank: ''
0010 %
0011 % The files expected in directory dirname are:
0012 % [prefix]CableTransmission.Txt
0013 % [prefix]0dBmSweep.Txt
0014 % [prefix]InjectedPower.txt
0015 % [prefix]SweepNL.txt
0016 % [prefix]SweepNR.txt
0017 %
0018 % The function will produce some plots and write relevant ones as pdf files
0019 % to directory dirname
0020 %
0021 % Oliver, 3 June 2011
0022 %
0023 
0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0025 % The default file names are:
0026 fn_cable = 'CableTransmission.txt';
0027 fn_0dBmSweep = '0dBmSweep.txt';
0028 fn_InjectedPower = 'InjectedPower.txt';
0029 fn_NLSweep = 'SweepNL.txt';
0030 fn_NRSweep = 'SweepNR.txt';
0031 fn_PolSweep = 'SweepPol.txt';
0032 
0033 if nargin == 2
0034     % current OVRO time zone?
0035     OVRO_time_zone = 7; % hours behind UTC
0036     dirname = varargin{1};
0037     prefix = varargin{2};
0038 else
0039     dirname = varargin{1};
0040     prefix = varargin{2};
0041     OVRO_time_zone = abs(varargin{3});
0042 end
0043 
0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0045 
0046 if ~strcmp(dirname(end),'/')
0047     dirname = [dirname '/'];
0048 end
0049 
0050 % now add the prefix to the dirname string:
0051 dirname = [dirname prefix];
0052 
0053 % The cable loss:
0054 dc = importdata([dirname fn_cable],'\t',1);
0055 fcable = dc.data(:,2);
0056 Pcable = dc.data(:,4)-dc.data(:,3);
0057 
0058 dcn = importdata([dirname fn_0dBmSweep],'\t',1);
0059 fnocable = dcn.data(:,2);
0060 Pnocable = dcn.data(:,4)-dcn.data(:,3);
0061 
0062 Gcable = Pcable-Pnocable;
0063 
0064 figure
0065 plot(fcable,Gcable,'k-')
0066 xlabel('Frequency [GHz]')
0067 ylabel('Cable transmission [dB]')
0068 
0069 % The actual power being put out by the sig gen:
0070 dsg = importdata([dirname fn_InjectedPower],'\t',1);
0071 
0072 fsiggen = dsg.data(:,2);
0073 SetPower = dsg.data(:,3);
0074 MeasuredPower = dsg.data(:,4);
0075 
0076 figure
0077 subplot(2,1,1)
0078 plot(fsiggen,SetPower,'k.',fsiggen,MeasuredPower,'r.')
0079 xlabel('Frequency [GHz]')
0080 ylabel('Sig gen power [dB]')
0081 legend('Set Power','Measured Power')
0082 subplot(2,1,2)
0083 plot(fsiggen,SetPower-MeasuredPower,'k.')
0084 xlabel('Frequency [GHz]')
0085 ylabel('Sig gen power error [dB]')
0086 
0087 saveas(gcf,[dirname 'SignalGeneratorBehavior.pdf'])
0088 
0089 % So, the actual power arriving at the input to the cryostat is:
0090 Pin = MeasuredPower+Gcable-30; % assume a flat -30dB of coupling
0091 
0092 figure
0093 plot(fcable,Pin,'k-')
0094 xlabel('Frequency')
0095 ylabel('Power in [dBm]')
0096 
0097 saveas(gcf,[dirname 'PowerInProfile.pdf'])
0098 
0099 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0100 % Now load the actual data from the telescope and the signal generator:
0101 
0102 % The signal generator data:
0103 dleft = importdata([dirname fn_NLSweep],'\t',1);
0104 tleft = dleft.data(:,1)+OVRO_time_zone/24;
0105 fleft = dleft.data(:,2);
0106 Pleft = dleft.data(:,3);
0107 dright = importdata([dirname fn_NRSweep],'\t',1);
0108 tright = dright.data(:,1)+OVRO_time_zone/24;
0109 fright = dright.data(:,2);
0110 Pright = dright.data(:,3);
0111 
0112 dpol = importdata([dirname fn_PolSweep],'\t',1);
0113 tpol = dpol.data(:,1)+OVRO_time_zone/24;
0114 fpol = dpol.data(:,2);
0115 Ppol = dpol.data(:,3);
0116 
0117 % The data to read in from the archive are:
0118 tstart = min([tleft; tright; tpol])-2/60/24; % 2 minute buffer
0119 tend = max([tleft; tright; tpol])+2/60/24;
0120 
0121 d = pipe_read(mjd2string(tstart),mjd2string(tend));
0122 
0123 % Calibrate the linearity of the measurements:
0124 d = calibrate_linearity(d);
0125 tdata = d.antenna0.receiver.utc;
0126 
0127 % The total power received at each diode:
0128 Pdata = d.antenna0.receiver.switchData(:,1:2:23) + d.antenna0.receiver.switchData(:,2:2:24);
0129 
0130 % For measuring the alpha value:
0131 P1 = d.antenna0.receiver.switchData(:,1)+d.antenna0.receiver.switchData(:,4); % "sky"
0132 P2 = d.antenna0.receiver.switchData(:,2)+d.antenna0.receiver.switchData(:,3); % "load"
0133 
0134 P11 = d.antenna0.receiver.switchData(:,21)+d.antenna0.receiver.switchData(:,24); % "sky"
0135 P12 = d.antenna0.receiver.switchData(:,22)+d.antenna0.receiver.switchData(:,23); % "load"
0136 
0137 % there are some throw-away events at the beginning and the end
0138 tlist = 5:(length(tleft)-4);
0139 
0140 % Declare some vectors we want the measured diode response to the injection
0141 % into LCP, and then into RCP.
0142 
0143 PmeasOL = zeros(length(tlist),12);
0144 PmeasOR = zeros(length(tlist),12);
0145 PmeasOP = zeros(length(tlist),12);
0146 
0147 P1L = zeros(length(tlist),1);
0148 P1R = zeros(length(tlist),1);
0149 P2L = zeros(length(tlist),1);
0150 P2R = zeros(length(tlist),1);
0151 
0152 % P3L = zeros(length(tlist),1);
0153 % P3R = zeros(length(tlist),1);
0154 % P4L = zeros(length(tlist),1);
0155 % P4R = zeros(length(tlist),1);
0156 %
0157 % P5L = zeros(length(tlist),1);
0158 % P5R = zeros(length(tlist),1);
0159 % P6L = zeros(length(tlist),1);
0160 % P6R = zeros(length(tlist),1);
0161 
0162 P11L = zeros(length(tlist),1);
0163 P11R = zeros(length(tlist),1);
0164 P12L = zeros(length(tlist),1);
0165 P12R = zeros(length(tlist),1);
0166 
0167 
0168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0169 % The time keeping is not perfect, so match up the injected power times and
0170 % the data times. We are looking for a maximum correlation between a test
0171 % signal set we generate and the data.
0172 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0173 
0174 disp('analyzePassbandMeasurement:: Identifying power injection events')
0175 tDel = 0.35; % the duration of the power events
0176 tOffsetList = -1.5:0.05:1.5;
0177 for k=1:length(tOffsetList)
0178     % tdata: the data of the time samples
0179     % tlist: the list of time events
0180     IselectionL = zeros(size(tdata));
0181     IselectionR = zeros(size(tdata));
0182     IselectionP = zeros(size(tdata));
0183     for m=1:2:(length(tlist)-1)
0184         % Every two events: the on events
0185         ttL = tleft(tlist(m)); % start time of left event
0186         ttR = tright(tlist(m)); % start time of right event
0187         ttP = tpol(tlist(m)); % start time of pol event
0188         IselectionL = IselectionL | (tdata-ttL) > tOffsetList(k)/60/60/24 & (tdata-ttL) < (tOffsetList(k)+tDel)/60/60/24;
0189         IselectionR = IselectionR | (tdata-ttR) > tOffsetList(k)/60/60/24 & (tdata-ttR) < (tOffsetList(k)+tDel)/60/60/24;
0190         IselectionP = IselectionP | (tdata-ttP) > tOffsetList(k)/60/60/24 & (tdata-ttP) < (tOffsetList(k)+tDel)/60/60/24;
0191     end
0192     
0193     % sum up all the data
0194     SumR(k) = sum(Pdata(IselectionR,1),1);
0195     SumL(k) = sum(Pdata(IselectionL,11),1);
0196     SumP(k) = sum(Pdata(IselectionP,5),1);
0197 end
0198 
0199 % what is the best offset for each?
0200 [y,IR] = max(SumR);
0201 [y,IL] = max(SumL);
0202 [y,IP] = max(SumP);
0203 
0204 offsetR = tOffsetList(IR);
0205 offsetL = tOffsetList(IL);
0206 offsetP = tOffsetList(IP);
0207 
0208 disp('analyzePassbandMeasurement::  Reading in the data')
0209 % and now read in the data with these offsets
0210 IsL = zeros(size(tdata));
0211 IsR = zeros(size(tdata));
0212 IsP = zeros(size(tdata));
0213 for k=1:length(tlist)
0214     ttL = tleft(tlist(k)); % start time of event
0215     ttR = tright(tlist(k)); % start time of event
0216     
0217     IselL = (tdata-ttL) > offsetL/60/60/24 & (tdata-ttL) < (offsetL+tDel)/60/60/24;
0218     IselR = (tdata-ttR) > offsetR/60/60/24 & (tdata-ttR) < (offsetR+tDel)/60/60/24;
0219     IselP = (tdata-ttP) > offsetP/60/60/24 & (tdata-ttP) < (offsetP+tDel)/60/60/24;
0220     
0221     % for plotting purposes
0222     IsL = IsL | IselL;
0223     IsR = IsR | IselR;
0224     IsP = IsP | IselP;
0225     
0226     PmeasOL(k,:) = mean(Pdata(IselL,:),1);
0227     PmeasOR(k,:) = mean(Pdata(IselR,:),1);
0228     PmeasOP(k,:) = mean(Pdata(IselP,:),1);
0229     
0230     P1L(k,:) = mean(P1(IselL,:),1);
0231     P1R(k,:) = mean(P1(IselR,:),1);
0232     
0233     P2L(k,:) = mean(P2(IselL,:),1);
0234     P2R(k,:) = mean(P2(IselR,:),1);
0235     
0236     P11L(k,:) = mean(P11(IselL,:),1);
0237     P11R(k,:) = mean(P11(IselR,:),1);
0238     
0239     P12L(k,:) = mean(P12(IselL,:),1);
0240     P12R(k,:) = mean(P12(IselR,:),1);
0241     
0242 end
0243 
0244 figure
0245 subplot(2,1,1)
0246 plot(tdata,Pdata(:,1),'k-',tdata(IsR),Pdata(IsR,1),'r.',tdata(IsP),Pdata(IsP,1),'g.')
0247 xlabel('Time [mjd]')
0248 ylabel('Channel 1')
0249 legend('Data','Identified events','Pol')
0250 subplot(2,1,2)
0251 plot(tdata,Pdata(:,11),'k-',tdata(IsL),Pdata(IsL,11),'r.',tdata(IsP),Pdata(IsP,11),'g.')
0252 xlabel('Time [mjd]')
0253 ylabel('Channel 11')
0254 legend('Data','Identified events','Pol')
0255 saveas(gcf,[dirname 'IdentifiedEvents.pdf'])
0256 
0257 % Now process the data to extract passbands: (on - off)
0258 PmeasuredL = PmeasOL(1:2:end-1,:) - PmeasOL(2:2:end,:);
0259 PmeasuredR = PmeasOR(1:2:end-1,:) - PmeasOR(2:2:end,:);
0260 PmeasuredP = PmeasOP(1:2:end-1,:) - PmeasOP(2:2:end,:);
0261 
0262 dP1L = P1L(1:2:end-1,:) - P1L(2:2:end,:);
0263 dP1R = P1R(1:2:end-1,:) - P1R(2:2:end,:);
0264 
0265 dP2L = P2L(1:2:end-1,:) - P2L(2:2:end,:);
0266 dP2R = P2R(1:2:end-1,:) - P2R(2:2:end,:);
0267 
0268 dP11L = P11L(1:2:end-1,:) - P11L(2:2:end,:);
0269 dP11R = P11R(1:2:end-1,:) - P11R(2:2:end,:);
0270 
0271 dP12L = P12L(1:2:end-1,:) - P12L(2:2:end,:);
0272 dP12R = P12R(1:2:end-1,:) - P12R(2:2:end,:);
0273 
0274 Pnom = Pleft(tlist);
0275 Pnominal = Pnom(1:2:end-1);
0276 
0277 fE = fleft(1:2:end-1); % the frequency of the event.
0278 
0279 % Now exclude measurements 51, 52, 153, and 154
0280 Iselect = [1:50 53:152 155:length(Pnominal)];
0281 
0282 mV2mW = 10^(-23/10)/3.6; % approximate conversion factor
0283 DU2mV = 2.5/10^4*1000/6.5; % convert DU to mV at the diode output
0284 dPL = PmeasuredL(Iselect,:)*DU2mV*mV2mW; % diode input in mW
0285 dPR = PmeasuredR(Iselect,:)*DU2mV*mV2mW; % diode input in mW
0286 dPP = PmeasuredP(Iselect,:)*DU2mV*mV2mW;
0287 
0288 dP1L = dP1L(Iselect,:)*DU2mV*mV2mW;
0289 dP2L = dP2L(Iselect,:)*DU2mV*mV2mW;
0290 dP1R = dP1R(Iselect,:)*DU2mV*mV2mW;
0291 dP2R = dP2R(Iselect,:)*DU2mV*mV2mW;
0292 
0293 dP11L = dP11L(Iselect,:)*DU2mV*mV2mW;
0294 dP12L = dP12L(Iselect,:)*DU2mV*mV2mW;
0295 dP11R = dP11R(Iselect,:)*DU2mV*mV2mW;
0296 dP12R = dP12R(Iselect,:)*DU2mV*mV2mW;
0297 
0298 Psg = 10.^(Pin/10); % the actualy signal generator power, in mW
0299 
0300 % And hence the gains are:
0301 GL = (abs(dPL./repmat(Psg(:),1,12)));
0302 GR = (abs(dPR./repmat(Psg(:),1,12)));
0303 GPol = (abs(dPP./repmat(Psg(:),1,12)));
0304 
0305 % If you want to exclude a particular data point from the plots because of
0306 % RFI, you can use the Ip vector to do so.
0307 [y,Ir] = min(abs(fcable-4.78));
0308 
0309 Ip = true(size(fcable));
0310  Ip(50) = false;
0311  Ip(150) = false;
0312 
0313 % figure
0314 % plot(fcable(Ip),GdL(Ip,2),'k-',fcable(Ip),GdR(Ip,2),'r-')
0315 % xlabel('Frequency [GHz]')
0316 % ylabel('Gain at output 2')
0317 % legend('From L','From R')
0318 
0319 % the alpha values are
0320 % aL1 = (dP1L./dP2L-1)./(dP1L./dP2L+1);
0321 aR1 = (dP1R./dP2R-1)./(dP1R./dP2R+1);
0322 %aR1 = (dP2R./dP1R-1)./(dP2R./dP1R+1);
0323 
0324 aL2 = (dP11L./dP12L-1)./(dP11L./dP12L+1);
0325 %aL2 = (dP12L./dP11L-1)./(dP12L./dP11L+1);
0326 % aR2 = (dP11R./dP12R-1)./(dP11R./dP12R+1);
0327 
0328 D3 = dPR(:,3)./dPL(:,3);
0329 
0330 gainR = (GR(:,1)+GR(:,2));
0331 gainL = (GL(:,11)+GL(:,12));
0332 % What is the average alpha value?
0333 aRm = sum(gainR.*aR1)/sum(gainR);
0334 aLm = sum(gainL.*aL2)/sum(gainL);
0335 
0336 Isel = gainR>max(gainR/8) & Ip(:);
0337 
0338 f = fcable(Ip);
0339 GI1 = (GR((Ip),1)+GR((Ip),2)+GL((Ip),1)+GL((Ip),2))/4;
0340 GI2 = (GL((Ip),11)+GL((Ip),12)+GR((Ip),11)+GR((Ip),12))/4;
0341 GP = (GR((Ip),3)+GR((Ip),4)+GL((Ip),3)+GL((Ip),4))/4;
0342 % BWI1 = mean(GI1/max(GI1))*(max(f)-min(f));
0343 BWI1v2 = (max(f)-min(f))/length(f)*(sum(GI1))^2/(sum(GI1.^2));
0344 % BWI2 = mean(GI2/max(GI2))*(max(f)-min(f));
0345 BWI2v2 = (max(f)-min(f))/length(f)*(sum(GI2))^2/(sum(GI2.^2));
0346 % BWP = mean(GP/max(GP))*(max(f)-min(f));
0347 BWPv2 = (max(f)-min(f))/length(f)*(sum(GP))^2/(sum(GP.^2));
0348 fcI1 = sum(f.*GI1)/(sum(GI1));
0349 fcI2 = sum(f.*GI2)/(sum(GI2));
0350 fcP = sum(f.*GP)/(sum(GP));
0351 
0352 
0353 % fprintf('I1: f_centre = %4.3f, BW = %4.3f\n',fcI1,BWI1)
0354 % fprintf('I2: f_centre = %4.3f, BW = %4.3f\n',fcI2,BWI2)
0355 % fprintf('P: f_centre = %4.3f, BW = %4.3f\n',fcP,BWP)
0356 fprintf('I1: f_centre = %4.3f, BW = %4.3f\n',fcI1,BWI1v2)
0357 fprintf('I2: f_centre = %4.3f, BW = %4.3f\n',fcI2,BWI2v2)
0358 fprintf('P: f_centre = %4.3f, BW = %4.3f\n',fcP,BWPv2)
0359 
0360 % Now fit to aR1 and aL2.
0361 disp('analyzePassbandMeasurement::  Fitting to alpha spectrum')
0362 % Try a number of initial values for the lag error:
0363 Loff = 30*(1:10);
0364 Lstart = [-fliplr(Loff) 0 Loff]/1000;
0365 
0366 % Fit for phase too?
0367 % 0 = fit lag and amplitude ratios separately
0368 % 1 = fit lag, phase, and amplitude ratios separately
0369 % 2 = fit lag, amplitude ratios, and diode 3 simultaneously
0370 fitFlag = 0;
0371 
0372 cR = {fcable(Isel)*1e9 aR1(Isel)};
0373 cL = {fcable(Isel)*1e9 aL2(Isel)};
0374 cS = {fcable(Isel)*1e9 aR1(Isel) aL2(Isel) D3(Isel)};
0375 xR = [];
0376 xL = [];
0377 xS = [];
0378 fLbest = 1e10;
0379 fRbest = 1e10;
0380 fSbest = 1e10;
0381 for k=1:length(Lstart)
0382     switch fitFlag
0383         case 1
0384             [xkR,fvalR] = fminsearch(@(x) ssqerrP(x,cR),[1;Lstart(k);0]);
0385             [xkL,fvalL] = fminsearch(@(x) ssqerrP(x,cL),[1;Lstart(k);0]);
0386         case 0
0387             [xkR,fvalR] = fminsearch(@(x) ssqerr(x,cR),[1;Lstart(k)]);
0388             [xkL,fvalL] = fminsearch(@(x) ssqerr(x,cL),[1;Lstart(k)]);
0389         case 2
0390             % x(1) = xR
0391             % x(2) = LR
0392             % x(3) = xL
0393             % x(4) = LL
0394             % x(5) = xP
0395             [xk,fval] = fminsearch(@(x) ssqerrSim(x,cS),[1;Lstart(k);1;Lstart(k);1]);
0396     end
0397     if fitFlag == 2
0398         if fval < fSbest
0399             xS = xk;
0400             fSbest = fval;
0401         end
0402     else
0403         if fvalR < fRbest
0404             xR = xkR;
0405             fRbest = fvalR;
0406         end
0407         if fvalL < fLbest
0408             xL = xkL;
0409             fLbest = fvalL;
0410         end
0411     end
0412 end
0413 
0414 % Now assign the derived best fits to xR and xL if we are in the
0415 % simultaneous case:
0416 switch fitFlag
0417     case 2
0418         fitPhase = 0;
0419         xR = xS(1:2);
0420         xL = xS(3:4);
0421     case 1
0422         fitPhase = 1;
0423     case 0
0424         fitPhase = 0;
0425 end
0426 
0427 if fitPhase
0428     fprintf('I1: Path error = %3.1f mm, Gain ratio = %3.2f, Phase error = %3.2f deg\n',xR(2)*1000,xR(1),xR(3)*180/pi)
0429     fprintf('I2: Path error = %3.1f mm, Gain ratio = %3.2f, Phase error = %3.2f deg\n',xL(2)*1000,xL(1),xL(3)*180/pi)
0430     aRfit = alpha(xR(1),xR(2),fcable*1e9,xR(3));
0431     aLfit = alpha(xL(1),xL(2),fcable*1e9,xL(3));
0432 else
0433     fprintf('I1: Path error = %3.1f mm, Gain ratio = %3.2f\n',xR(2)*1000,xR(1))
0434     fprintf('I2: Path error = %3.1f mm, Gain ratio = %3.2f\n',xL(2)*1000,xL(1))
0435     aRfit = alpha(xR(1),xR(2),fcable*1e9,0);
0436     aLfit = alpha(xL(1),xL(2),fcable*1e9,0);
0437 end
0438 
0439 
0440 figure
0441 plot(fcable(Isel),aR1(Isel),'k.',fcable(Isel),ones(size(fcable(Isel)))*aRm,'k--',...
0442     fcable(Isel),aL2(Isel),'r.',fcable(Isel),ones(size(fcable(Isel)))*aLm,'r--',...
0443     fcable,aRfit,'k-',fcable,aLfit,'r-',...
0444     fcable,aR1,'k-.',fcable,aL2,'r-.')
0445 xlabel('Frequency [GHz]')
0446 ylabel('\alpha [goal = 1]')
0447 legend('\alpha_{I1}','mean','\alpha_{I2}','mean','Location','East')
0448 ylim([-1.05 1.05])
0449 if fitPhase
0450     sI1 = sprintf('I1: L=%3.1fmm, x=%3.2f, theta=%3.2fdeg',xR(2)*1000,xR(1),xR(3)*180/pi);
0451     sI2 = sprintf('I2: L=%3.1fmm, x=%3.2f, theta=%3.2fdeg',xL(2)*1000,xL(1),xL(3)*180/pi);
0452 else
0453     sI1 = sprintf('I1: L=%3.1fmm, x=%3.2f',xR(2)*1000,xR(1));
0454     sI2 = sprintf('I2: L=%3.1fmm, x=%3.2f',xL(2)*1000,xL(1));
0455 end
0456 
0457 title([sI1 '; ' sI2])
0458 
0459 saveas(gcf,[dirname 'alpha_Spectrum.pdf'])
0460 
0461 % and write the CW alpha values to disk:
0462 fn = [dirname 'alpha_cw.csv'];
0463 f1 = fopen(fn,'w');
0464 fprintf(f1,'# f [GHz], alpha_cw(I1), alpha_cw(I2), gain_I1 [dB], gain_I2 [dB]\n');
0465 fprintf(f1,'%4.3f, %4.3f, %4.3f, %4.3f, %4.3f\n',[fcable aR1 aL2 10*log10(gainR) 10*log10(gainL)]');
0466 fclose(f1);
0467 
0468 % Write the passbands to disk:
0469 fpoints = fcable(Ip);
0470 G1 = (GR(Ip,1)+GR(Ip,2)+GL(Ip,1)+GL(Ip,2))/4;
0471 G2 = (GR(Ip,3)+GR(Ip,4)+GL(Ip,3)+GL(Ip,4))/4;
0472 G3 = (GR(Ip,5)+GR(Ip,6)+GL(Ip,5)+GL(Ip,6))/4;
0473 G4 = (GR(Ip,7)+GR(Ip,8)+GL(Ip,7)+GL(Ip,8))/4;
0474 G5 = (GR(Ip,9)+GR(Ip,10)+GL(Ip,9)+GL(Ip,10))/4;
0475 G6 = (GL(Ip,11)+GL(Ip,12)+GR(Ip,11)+GR(Ip,12))/4;
0476 fn = [dirname 'Passband.csv'];
0477 f1 = fopen(fn,'w');
0478 fprintf(f1,'# f [GHz], Gain I1, Gain Q1, Gain U1, Gain Q2, Gain U2, Gain I2\n');
0479 fprintf(f1,'%4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f\n',[fpoints G1 G2 G3 G4 G5 G6]');
0480 fclose(f1);
0481 
0482 % Plot the data!
0483 ylims = [20 90];
0484 
0485 figure
0486 plot(fcable(Ip),10*log10((GR(Ip,1)+GR(Ip,2)+GL(Ip,1)+GL(Ip,2))/4),'k-',...
0487     fcable(Ip),10*log10((GR(Ip,3)+GR(Ip,4)+GL(Ip,3)+GL(Ip,4))/4),'r-',...
0488     fcable(Ip),10*log10((GR(Ip,5)+GR(Ip,6)+GL(Ip,5)+GL(Ip,6))/4),'b-',...
0489     fcable(Ip),10*log10((GR(Ip,7)+GR(Ip,8)+GL(Ip,7)+GL(Ip,8))/4),'r--',...
0490     fcable(Ip),10*log10((GR(Ip,9)+GR(Ip,10)+GL(Ip,9)+GL(Ip,10))/4),'b--',...
0491     fcable(Ip),10*log10((GL(Ip,11)+GL(Ip,12)+GR(Ip,11)+GR(Ip,12))/4),'k--')
0492 ylim(ylims)
0493 xlabel('Frequency [GHz]')
0494 ylabel('Gain [dB]')
0495 legend('Channel 1','Channel 2','Channel 3','Channel 4','Channel 5',...
0496     'Channel 6','Location','SouthEast')
0497 title([sprintf('f_c [I1] = %4.3f GHz, BW = %4.3f GHz',fcI1,BWI1v2) ';' sprintf('f_c [I2] = %4.3f GHz, BW = %4.3f GHz',fcI2,BWI2v2) ';' sprintf('f_c [P] = %4.3f GHz, BW = %4.3f GHz',fcP,BWPv2)])
0498 grid on
0499 grid minor
0500 
0501 
0502 
0503 figure
0504 plot(fcable(Ip),10*log10((GPol(Ip,1)+GPol(Ip,2))/2),'k-',...
0505     fcable(Ip),10*log10((GPol(Ip,3)+GPol(Ip,4))/2),'r-',...
0506     fcable(Ip),10*log10((GPol(Ip,5)+GPol(Ip,6))/2),'b-',...
0507     fcable(Ip),10*log10((GPol(Ip,7)+GPol(Ip,8))/2),'r--',...
0508     fcable(Ip),10*log10((GPol(Ip,9)+GPol(Ip,10))/2),'b--',...
0509     fcable(Ip),10*log10((GPol(Ip,11)+GPol(Ip,12))/2),'k--')
0510 ylim(ylims)
0511 xlabel('Frequency [GHz]')
0512 ylabel('Gain [dB]')
0513 legend('Channel 1','Channel 2','Channel 3','Channel 4','Channel 5',...
0514     'Channel 6','Location','SouthEast')
0515 title('Pure U input')
0516 grid on
0517 grid minor
0518 
0519 saveas(gcf,[dirname 'C-BASS_Passband.pdf'])
0520 
0521 % Plot the response of all the diodes:
0522 figure
0523 subplot(3,4,1)
0524 plot(fcable(Ip),10*log10(GR(Ip,1)),'k-',fcable(Ip),10*log10(GL(Ip,1)),'r-')
0525 ylabel('Amplitude [dB]')
0526 ylim(ylims)
0527 title('Diode 1')
0528 subplot(3,4,2)
0529 plot(fcable(Ip),10*log10(GR(Ip,2)),'k-',fcable(Ip),10*log10(GL(Ip,2)),'r-')
0530 ylim(ylims)
0531 title('Diode 2')
0532 subplot(3,4,3)
0533 plot(fcable(Ip),10*log10(GR(Ip,11)),'k-',fcable(Ip),10*log10(GL(Ip,11)),'r-')
0534 ylim(ylims)
0535 title('Diode 11')
0536 subplot(3,4,4)
0537 plot(fcable(Ip),10*log10(GR(Ip,12)),'k-',fcable(Ip),10*log10(GL(Ip,12)),'r-')
0538 ylim(ylims)
0539 title('Diode 12')
0540 
0541 subplot(3,4,5)
0542 plot(fcable(Ip),10*log10(GR(Ip,3)),'k-',fcable(Ip),10*log10(GL(Ip,3)),'r-')
0543 ylim(ylims)
0544 ylabel('Amplitude [dB]')
0545 title('Diode 3')
0546 subplot(3,4,6)
0547 plot(fcable(Ip),10*log10(GR(Ip,4)),'k-',fcable(Ip),10*log10(GL(Ip,4)),'r-')
0548 ylim(ylims)
0549 title('Diode 4')
0550 subplot(3,4,7)
0551 plot(fcable(Ip),10*log10(GR(Ip,7)),'k-',fcable(Ip),10*log10(GL(Ip,7)),'r-')
0552 ylim(ylims)
0553 title('Diode 7')
0554 subplot(3,4,8)
0555 plot(fcable(Ip),10*log10(GR(Ip,8)),'k-',fcable(Ip),10*log10(GL(Ip,8)),'r-')
0556 ylim(ylims)
0557 title('Diode 8')
0558 
0559 subplot(3,4,9)
0560 plot(fcable(Ip),10*log10(GR(Ip,5)),'k-',fcable(Ip),10*log10(GL(Ip,5)),'r-')
0561 ylim(ylims)
0562 ylabel('Amplitude [dB]')
0563 xlabel('Frequency [GHz]')
0564 title('Diode 5')
0565 subplot(3,4,10)
0566 plot(fcable(Ip),10*log10(GR(Ip,6)),'k-',fcable(Ip),10*log10(GL(Ip,6)),'r-')
0567 ylim(ylims)
0568 xlabel('Frequency [GHz]')
0569 title('Diode 6')
0570 subplot(3,4,11)
0571 plot(fcable(Ip),10*log10(GR(Ip,9)),'k-',fcable(Ip),10*log10(GL(Ip,9)),'r-')
0572 ylim(ylims)
0573 xlabel('Frequency [GHz]')
0574 title('Diode 9')
0575 subplot(3,4,12)
0576 plot(fcable(Ip),10*log10(GR(Ip,10)),'k-',fcable(Ip),10*log10(GL(Ip,10)),'r-')
0577 ylim(ylims)
0578 xlabel('Frequency [GHz]')
0579 legend('RCP input','LCP input','Location','South')
0580 title('Diode 10')
0581 
0582 saveas(gcf,[dirname 'C-BASS_Passband_Diode_Responses.pdf'])
0583 
0584 end
0585 
0586 function err = ssqerr(x,c)
0587 % Function to fit the alpha value curve to:
0588 v = c{1}; % frequencies of the alpha points in Hz
0589 a = c{2}; % the alpha values
0590 
0591 % x(1) is ratio of voltage gains, x(2) is path error
0592 afit = alpha(x(1),x(2),v,0);
0593 err = sum((afit-a).^2);
0594 end
0595 
0596 function err = ssqerrP(x,c)
0597 % Function to fit the alpha value curve to:
0598 v = c{1}; % frequencies of the alpha points in Hz
0599 a = c{2}; % the alpha values
0600 
0601 % x(1) is ratio of voltage gains, x(2) is path error, x(3) is phase offset
0602 afit = alpha(x(1),x(2),v,x(3));
0603 err = sum((afit-a).^2);
0604 end
0605 
0606 function err = ssqerrSim(x,c)
0607 % Simultaneous fit to left, right, and pol channels
0608 v = c{1};
0609 aR = c{2};
0610 aL = c{3};
0611 D3 = c{4};
0612 
0613 % x(1) = xR
0614 % x(2) = LR
0615 % x(3) = xL
0616 % x(4) = LL
0617 % x(5) = xP
0618 
0619 aRfit = alpha(x(1),x(2),v,0);
0620 aLfit = alpha(x(3),x(4),v,0);
0621 D3fit = D3ratio(x(5),x(3),x(1),x(4),x(2),v);
0622 
0623 err = sum((aR-aRfit).^2+(aL-aLfit).^2+(D3-D3fit).^2);
0624 
0625 end
0626 
0627 function y = alpha(x,L,v,toff)
0628 % x = |g1|/|g2|
0629 % L = path length difference in m
0630 % v = frequency in Hz
0631 c = 3e8;
0632 Lc = L/c;
0633 
0634 y = 2*x*cos(2*pi*v*Lc+toff)/(1+x^2);
0635 end
0636 
0637 function y = D3ratio(xP,xL,xR,LL,LR,v)
0638 % the ratio of the diode 3 output powers
0639 y = xP^2*(1+xR^2+2*xR*cos(2*pi*v*LR/3e8))./(1+xL^2+2*xL*cos(2*pi*v*LL/3e8));
0640 
0641 end
analyzePassbandMeasurement

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
