first commit after restructuring

projects/dl814_k7/dl814_k7_cal_coarse/data/20230326-天线幅相测试/multiBitCorrInterSampleDiff.m

@@ -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

projects/dl814_k7/dl814_k7_cal_coarse/data/20230326-天线幅相测试/suqareShapeFreqEst.m

@@ -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