finish pa phase calibration, add sample data

测试了三台对飞节点的相对幅相,测试了一台814节点的相对幅相
2023-05-16 09:14:53 +08:00
parent c32d684223
commit 986a2ea1e1
256 changed files with 1273 additions and 43 deletions
@@ -0,0 +1,123 @@
+%% 幅相测试结果离线导出
+clc; clear;
+%% 初始参数
+captureNum = 2;     % 同一频点捕获数
+channelNum = 4;
+% freqnum    = 90;    % 频点
+% dataLen = captureNum*freqnum;  
+dataFile = "0508/iladata_freq_comp_3.csv";
+sys_clk = 48e6;     % 系统时钟
+phase_width =   32; % 频率控制字总长度
+
+
+%% data analysis
+rx_string = readmatrix(dataFile,'OutputType','string'); % 含有字符串的矩阵，便于根据STATE确定数据位置
+rxdata = readmatrix(dataFile);
+ch1_dmf = [];
+ch2_dmf = [];
+ch3_dmf = [];
+ch4_dmf = [];
+ch1_begin = [];
+ch2_begin = [];
+ch3_begin = [];
+ch4_begin = [];
+pinc = [];
+
+% 找出ch1~ch4进行captureNum次DMF输出的相关峰序列、每次DMF开始时对应的采样序号
+% 和每次DMF补偿频偏的频率控制字
+for k = 1 : length(rx_string(:,1))
+    if (rx_string(k,13) == 'STATE_CH1')
+        ch1_dmf = [ch1_dmf rxdata(k,12)];
+        if(rx_string(k-1,13) == 'STATE_CH_SW')
+            ch1_begin = [ch1_begin;k];
+            pinc = [pinc;rxdata(k+1,17)];
+        end
+    elseif (rx_string(k,13) == 'STATE_CH2')
+        ch2_dmf = [ch2_dmf rxdata(k,12)];
+        if(rx_string(k-1,13) == 'STATE_CH_SW')
+            ch2_begin = [ch2_begin;k];
+        end
+    elseif (rx_string(k,13) == 'STATE_CH3')
+        ch3_dmf = [ch3_dmf rxdata(k,12)];
+        if(rx_string(k-1,13) == 'STATE_CH_SW')
+            ch3_begin = [ch3_begin;k];
+        end
+    elseif (rx_string(k,13) == 'STATE_CH4')
+        ch4_dmf = [ch4_dmf rxdata(k,12)];
+        if(rx_string(k-1,13) == 'STATE_CH_SW')
+            ch4_begin = [ch4_begin;k];
+        end
+    end
+end
+% 对DMF序列以捕获次数分行
+ch1_dmf = (reshape(ch1_dmf,[length(ch1_dmf)/captureNum,captureNum])).';
+ch2_dmf = (reshape(ch2_dmf,[length(ch2_dmf)/captureNum,captureNum])).';
+ch3_dmf = (reshape(ch3_dmf,[length(ch3_dmf)/captureNum,captureNum])).';
+ch4_dmf = (reshape(ch4_dmf,[length(ch4_dmf)/captureNum,captureNum])).';
+
+
+ch1_max = zeros(captureNum,1);
+ch2_max = zeros(captureNum,1);
+ch3_max = zeros(captureNum,1);
+ch4_max = zeros(captureNum,1);
+ch1_pos = zeros(captureNum,1);
+ch2_pos = zeros(captureNum,1);
+ch3_pos = zeros(captureNum,1);
+ch4_pos = zeros(captureNum,1);
+% 找出DMF输出相关峰后半段序列中最大峰对应序号
+for kk = 1:captureNum
+    [ch1_max(kk),ch1_pos(kk)] = max(ch1_dmf(kk,ceil(length(ch1_dmf(kk,:))/2):end));
+    [ch2_max(kk),ch2_pos(kk)] = max(ch2_dmf(kk,ceil(length(ch2_dmf(kk,:))/2):end));
+    [ch3_max(kk),ch3_pos(kk)] = max(ch3_dmf(kk,ceil(length(ch3_dmf(kk,:))/2):end));
+    [ch4_max(kk),ch4_pos(kk)] = max(ch4_dmf(kk,ceil(length(ch4_dmf(kk,:))/2):end));
+    ch1_pos(kk) = ch1_pos(kk) + ceil(length(ch1_dmf(kk,:))/2) - 1;
+    ch2_pos(kk) = ch2_pos(kk) + ceil(length(ch2_dmf(kk,:))/2) - 1;  
+    ch3_pos(kk) = ch3_pos(kk) + ceil(length(ch3_dmf(kk,:))/2) - 1;
+    ch4_pos(kk) = ch4_pos(kk) + ceil(length(ch4_dmf(kk,:))/2) - 1;
+end
+% 计算ch2~ch4相关峰相对于ch1偏移的位置。ch_begin为DMF序列在ila中的开始位置，
+% pos为DMF序列中最大值的位置，相加为最大值在整个ila中的决定位置。插值为载波偏移的采样点数
+ch2_pos_shift = (ch2_pos + ch2_begin) - (ch1_pos + ch1_begin);
+ch3_pos_shift = (ch3_pos + ch3_begin) - (ch1_pos + ch1_begin);
+ch4_pos_shift = (ch4_pos + ch4_begin) - (ch1_pos + ch1_begin);
+
+% 计算频偏补偿载波的(频率/fs)
+f_out   = 2 * pi * pinc / (2^phase_width);
+
+% 计算串行捕获带来的相位偏移
+ch2_phase_comp = ch2_pos_shift .* f_out;
+ch3_phase_comp = ch3_pos_shift .* f_out;
+ch4_phase_comp = ch4_pos_shift .* f_out;
+
+PeakValue_ch1 = [];
+PeakValue_ch2 = [];
+PeakValue_ch3 = [];
+PeakValue_ch4 = [];
+
+% 根据ch_peak_valid信号找到硬件捕获点的信号
+for kkk = 1:length(rx_string(:,1))
+    if(rxdata(kkk,16) == 1)
+        PeakValue_ch1 = [PeakValue_ch1 (rxdata(kkk+1, 4) + 1i*rxdata(kkk+1, 5))]; % 用rxdata(kkk, 4)+1i*rxdata(kkk, 5)也可以
+        PeakValue_ch2 = [PeakValue_ch2 (rxdata(kkk+1, 6) + 1i*rxdata(kkk+1, 7))];
+        PeakValue_ch3 = [PeakValue_ch3 (rxdata(kkk+1, 8) + 1i*rxdata(kkk+1, 9))];
+        PeakValue_ch4 = [PeakValue_ch4 (rxdata(kkk+1, 10) + 1i*rxdata(kkk+1, 11))];
+    end
+end
+
+% 计算相对ch1偏移的相位（考虑串行捕获导致的频偏校准载波相位偏移）
+for kkkk = 1:captureNum
+  angleEstSym(kkkk,1) = mod((angle(PeakValue_ch2(kkkk))-angle(PeakValue_ch1(kkkk))-ch2_phase_comp(kkkk))/pi*180,360);
+  angleEstSym(kkkk,2) = mod((angle(PeakValue_ch3(kkkk))-angle(PeakValue_ch1(kkkk))-ch3_phase_comp(kkkk))/pi*180,360);
+  angleEstSym(kkkk,3) = mod((angle(PeakValue_ch4(kkkk))-angle(PeakValue_ch1(kkkk))-ch4_phase_comp(kkkk))/pi*180,360);
+end
+disp("考虑串行捕获导致的频偏校准载波相位偏移时")
+angleEstSym
+
+% 计算相对ch1偏移的相位（不考虑串行捕获导致的频偏校准载波相位偏移）
+for kkkk = 1:captureNum
+  angleEstSym(kkkk,1) = mod((angle(PeakValue_ch2(kkkk))-angle(PeakValue_ch1(kkkk)))/pi*180,360);
+  angleEstSym(kkkk,2) = mod((angle(PeakValue_ch3(kkkk))-angle(PeakValue_ch1(kkkk)))/pi*180,360);
+  angleEstSym(kkkk,3) = mod((angle(PeakValue_ch4(kkkk))-angle(PeakValue_ch1(kkkk)))/pi*180,360);
+end
+disp("不考虑串行捕获导致的频偏校准载波相位偏移时")
+angleEstSym
@@ -0,0 +1,108 @@
+%% 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 = 8192;
+fs = 48e6;
+dataFile = "0508/iladata_rx_data.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(:, 4) + 1i*rxdata(:, 5)).';
+seqRx(:,2) = (rxdata(:, 6) + 1i*rxdata(:, 7)).';
+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,-30000,30000,100);
+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])
+
+  subplot(122);
+  plot(angleEstSym(:,[1,4]));
+  legend({"ch1","ch4"});
+  title("Angle Estimation");
+  ylim([min(angleEstSym(:,[1,4]),[],'all')-10,max(angleEstSym(:,[1,4]),[],'all')+10])
+
+else
+
+  subplot(121);
+  plot(ampEstSym);
+  legend({"ch1","ch2","ch3","ch4"});
+  title("Amplitude Estimation");
+  ylim([min(ampEstSym,[],'all')-0.2,max(ampEstSym,[],'all')+0.2])
+
+  subplot(122);
+  plot(angleEstSym);
+  legend({"ch1","ch2","ch3","ch4"});
+  title("Angle Estimation");
+  ylim([min(angleEstSym,[],'all')-10,max(angleEstSym,[],'all')+10])
+
+end
+
+phase_sym = (360 - mean(angleEstSym))/180;
+phase = fix(phase_sym * (2^9));
+angleEstSym
@@ -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
--- a/Show More
+++ b/Show More