update dl814_z7 and dl814_cal_coarse

change frame header of nulling algorithm data in chip-to-chip communication, add several sets of noise data sampled with 814 board
2023-04-05 15:17:14 +08:00
parent d2c22dbeb5
commit 5402f12cf2
122 changed files with 1966442 additions and 24591 deletions
@@ -0,0 +1,109 @@
+%% File Header/文件说明
+%
+% Beamforming with phased array derives signal characteristics from phase
+% differences between sensors. Therefore the amplitude-phase consistence is
+% vital to phased systems.
+%
+% The script analyzes the amp-phase relations of ILA captured data, the format
+% of which is square-shaped pseudorandorm bit sequence. Each channel is
+% processed seperately with different frequency shifts and different optimum
+% sampling points. The estimation of amplitude and phase are based on
+% acculumated correlation peaks.
+%
+clc; clear;
+%% initial parameters/初始参数
+interp = 8; % interpolation factor
+channelNum = 4;
+dataLen = 131072;
+fs = 48e6;
+dataFile = "iladata_ant_3-1.csv";
+
+% standard pn sequence
+pnSequence = comm.PNSequence('Polynomial','x^7+x^6+1', ...
+  'InitialConditions',[zeros(1,6),1],'SamplesPerFrame',2^7-1);
+PRBS7 = pnSequence();
+PRBS7_polar = 1-2*[PRBS7;0];
+
+chipLen = length(PRBS7_polar);    % pn chip length
+symNum = dataLen/chipLen/interp;  % total symbol number
+symNumPerAccum = 1;               % symbol nbumber of each accumulation
+accNum = symNum/symNumPerAccum;   % total accumulation numbers, i.e. peak numbers
+
+%% data analysis
+
+% read data, and swap ch1 and ch2, because loop data locate at ch2 and ch4
+rxdata = readmatrix(dataFile);
+seqRx(:,1) = (rxdata(:, 6) + 1i*rxdata(:, 7)).';
+seqRx(:,2) = (rxdata(:, 4) + 1i*rxdata(:, 5)).';
+seqRx(:,3) = (rxdata(:, 8) + 1i*rxdata(:, 9)).';
+seqRx(:,4) = (rxdata(:,10) + 1i*rxdata(:,11)).';
+
+% compensate frequecny deviation and downsample at optimum sample position based
+% on the varience of constellation
+if dataFile.contains("loop")
+  [seqRxDown,fd,samplePos] = suqareShapeFreqEst(seqRx,interp,fs,-50,50,5);
+else
+  [seqRxDown,fd,samplePos] = suqareShapeFreqEst(seqRx,interp,fs,24900,25100,1);
+end
+
+% split [seqRxDown] into multiple groups, each group contains several symbols to
+% be correlated with PRBS7 and accumulated
+seqRxDownSplit = reshape(seqRxDown(1:dataLen/interp,:),chipLen*symNumPerAccum, ...
+  accNum,channelNum);
+
+seqRxAcc = zeros(chipLen,accNum,channelNum);
+desiredPnt = zeros(accNum,channelNum);
+ampEstSym = zeros(accNum,channelNum);
+angleEstSym = zeros(accNum,channelNum);
+
+for i = 1:accNum
+  for j = 1:channelNum
+    % correlation & multi-symbol accumulation
+    groupCorrAcc = sum(reshape(conv(flip(seqRxDownSplit(:,i,j)), ...
+      PRBS7_polar,'same'),chipLen,symNumPerAccum),2)';
+    seqRxAcc(:,i,j) = groupCorrAcc;
+    % peak position
+    desiredPnt(i,j) = groupCorrAcc(abs(groupCorrAcc)==max(abs(groupCorrAcc)));
+  end
+end
+
+for i = 1:channelNum
+  ampEstSym(:,i) = abs(desiredPnt(:,i))./abs(desiredPnt(:,1));
+  angleEstSym(:,i) = mod((angle(desiredPnt(:,i))-angle(desiredPnt(:,1)))/pi*180,360);
+end
+
+%% figures
+figure(1);
+if dataFile.contains("loop")
+
+  subplot(121);
+  plot(ampEstSym(:,[1,4]));
+  legend({"ch1","ch4"});
+  title("Amplitude Estimation");
+  ylim([min(ampEstSym(:,[1,4]),[],'all')-0.1,max(ampEstSym(:,[1,4]),[],'all')+0.1]);
+  grid on; grid minor;
+
+  subplot(122);
+  plot(angleEstSym(:,[1,4]));
+  legend({"ch1","ch4"});
+  title("Angle Estimation");
+  ylim([min(angleEstSym(:,[1,4]),[],'all')-1,max(angleEstSym(:,[1,4]),[],'all')+1]);
+  grid on; grid minor;
+
+else
+
+  subplot(121);
+  plot(ampEstSym);
+  legend({"ch1","ch2","ch3","ch4"});
+  title("Amplitude Estimation");
+  ylim([min(ampEstSym,[],'all')-0.2,max(ampEstSym,[],'all')+0.2]);
+  grid on; grid minor;
+
+  subplot(122);
+  plot(angleEstSym);
+  legend({"ch1","ch2","ch3","ch4"});
+  title("Angle Estimation");
+  ylim([min(angleEstSym,[],'all')-10,max(angleEstSym,[],'all')+10]);
+  grid on; grid minor;
+
+end
@@ -0,0 +1,50 @@
+function [dataOut, fdTarget, samplePosTarget] = suqareShapeFreqEst(dataIn, interp, fs, fmin, fmax, fd)
+%suqareShapeFreqEst - 利用星座图的集中程度在一定频率范围内使用频率步进估计方波成形信号的频偏
+%
+% data - input complex data, row vector
+% dataLen - data length
+% interp - interpolation rate
+% fs - sampling rate, in Herts
+% fmin - minimum estimation frequency, in Herts
+% fmax - maximum estimation frequency, in Herts
+% fd - estimation frequency step, in Herts
+  
+% dimension
+[dataLen,channelNum] = size(dataIn);
+
+freqList = fmin:fd:fmax; 
+freqNum = length(freqList);
+varRecord = zeros(freqNum,interp);
+
+dataOut = zeros(dataLen/interp,channelNum);
+fdTarget = zeros(channelNum,1);
+samplePosTarget = zeros(channelNum,1);
+
+for channel = 1:channelNum
+
+  for freq = 1:freqNum
+
+    % carrier
+    cxr = exp(-1i*2*pi*freqList(freq)/fs*(0:dataLen-1)');
+
+    % apply frequency shift onto input data
+    seqRx = dataIn(:,channel).*cxr;
+  
+    for samplePos = 0:interp-1
+      seqRxDown = downsample(seqRx,interp,samplePos); % 下采样
+      seqRxDownQuaterOne = abs(real(seqRxDown)) + 1i*abs(imag(seqRxDown)); % 集中到第一象限
+      varRecord(freq,samplePos+1) = var(seqRxDownQuaterOne);
+    end
+  
+  end
+
+  [fdIndex, samplePosIndex] = find(varRecord==min(min(varRecord)));
+  fdTarget(channel) = freqList(fdIndex);
+  samplePosTarget(channel) = samplePosIndex-1;
+
+  cxrTarget = exp(-1i*2*pi*fdTarget(channel)/fs*(0:dataLen-1)');
+  seqRxTarget = dataIn(:,channel).*cxrTarget;
+  dataOut(:,channel) = downsample(seqRxTarget,interp,samplePosTarget(channel));
+end
+
+end