毕业课题 | ISAC通感一体化平台数据处理matlab代码 信道估计 时频域特征提取 有监督机器学习检测分类

-----------------------------------目录-------------------------------------

读取接收机数据(GSR_Rx_tdms_read.m) -->

获取信道信息(用导频做信道估计HestiPhase.m) -->

时域频域特征提取，存入excel表格 (基础特征&基于相关论文的重要特征，timeFeatureExtraction.m; HstftFeatureExtraction.m) -->

从excel中读取数据，训练分类器(newInputGenerator.m) -->

训练机器学习分类器

----------------------------------代码----------------------------------------

读取数据代码 GSR_Rx_tdms_read.m

%%
%matlab版本，Matlab 2022b及以上
clc;clear;
%文件路径(从实验设备/接收机中获取的数据)
fileName = "C:\Users\Z\Desktop\data\0201\5liandun1.tdms";

%GSR默认的通道组名称为USRP Rx，若打开其他TDMS文件，需要确认通道组名称是否正确
ChannelGroup="USRP Rx";

%每次读取的数据长度
readSize = 60000000000;

%读取通道组中的ChannelGroup通道
ds = tdmsDatastore(fileName, ReadSize=readSize, SelectedChannelGroup=ChannelGroup);

%从文件头开始读取通道组中所有通道的数据，每次读取点数为readsize
while(hasdata(ds))
    [IQ_cell,info] = read(ds);  
    info.Offset %
end

IQ=table2array(IQ_cell{1});
chanenel_num = length(IQ(1,:));
sample_num = length(IQ(:,1));
% IQ数据为uint32，高16bit为I路，低16位为Q，拆分后数据类型为int16
channel_data = zeros(sample_num,chanenel_num);
for i = 1:chanenel_num
    Q = typecast((bitshift(IQ(:,i),-16)),'int16');
    I = typecast(bitand(IQ(:,i),uint32(65535)),'int16');
    
    %去偶数行的0和数据对齐，可能不用去？不用就注释掉
    I = I(1:2:end);
    Q = Q(1:2:end);

    channel_data(:,i) = complex(I, Q);    
end

%plot(Q)

信道估计代码 HestiPhase.m

% -------用导频H做phase
% --------------OFDM parameters-----------
Carriers = 3276;
Nfft = 4096; % 进行了fft操作的点数，满足2^n
Np = 4096;
%% --------------------------------获取导频的位置参数loc--------------------------
% 16 - QAM - Modulation Scheme 
Es = 1; 
M = 64;
A = sqrt(3/2/(M-1)*Es);% Signal energy and QAM normalization factor
Nps = 4; % Pilot Spacing
msgint = randi([2 5],1,Carriers); % 模拟生成发送信号
dat_ser = A*qammod(msgint,M); % 信号调制
% serial to parllel conversion
dat_par = dat_ser.'; % 将序列转换为并行形式，存储在dat_par中。
% Pilot Insertion - Comb Type Arrangement
j = 1;
k = 1; 
loc = [];
for i = 1:Nfft
    if or((i<=round(Carriers/2)),(i> Nfft-round(Carriers/2)))
        if mod(i,Nps) == 1 % 余数为1
            X(i) = k;    % 插入导频
            loc = [loc i];   % 记录导频位置
            k=k+1;
        end
    end
end


%% --------------------------------获取发送时真实的导频----------------------------
XP_data = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\Tx_Pilot_Xp.xlsx');
real_pilot = XP_data(1:2048);
imag_pilot = XP_data(2049:end);
Xp = complex(real_pilot,imag_pilot);

Nps = 4; % Pilot Spacing
Np = floor(Carriers/Nps); % 每个OFDM信号的导频载波个数


%% ---------发送信号做接收，用于验证导频位置-------------------------------
% xdata = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\xdata.xlsx');
% real_xdata = xdata(1:4096);
% imag_xdata = xdata(4097:end);
% qam_x_data = complex(real_xdata,imag_xdata); % 全部信号：1*4096
% % inverse discret Fourier transform (IFFT)
% X_ifft = ifft(qam_x_data,Nfft);
% % 加循环前缀
% guard = X_ifft(:,end-1024+1:end); % 取出后1/4loop_num
% ofdm = [guard X_ifft];
% zero_tianchong = zeros(1,1074);
% zero_buchong = zeros(1,1998);
% rx_ofdm = [zero_tianchong ofdm zero_buchong];

%% ---------------------接收信号处理，取出导频，得H矩阵------------------
loop_num = 80;
% ---获取第一通道的数据, 1638400*1
Rx_Total_signal = channel_data(:,1); 
% Rx_Total_signal = rx_ofdm; % 收发一致时验证位置
% ---变为矩阵形式，每一列代表一次loop, 8192*200
Rx_Total_matrix = reshape(Rx_Total_signal, [], size(channel_data, 1)/8192);
% ---提取有效的loop，可以通过cpcc.m相关性验证，一般为第30次-200次, 8192*170
Rx_eff_loop_matrix = Rx_Total_matrix(:,31:end); % 8192*170
% Rx_eff_loop_matrix = Rx_Total_matrix; % 8192*170

% ---在每个loop的8192中点中，提取相关的点。cpcc.m中的峰值代表响应的-起始-
% ---第一个峰值一般在1075出现，则1075-2098是循环前缀，2099-6194是发送信号。(有效信号长4096，取后1/4做循环前缀)
% ---4096*170
Rx_eff_point_matrix = Rx_eff_loop_matrix(2099:6194,:);
% ---fft调制到频域, 819*170
%Rx_fre_matrix = fft(Rx_eff_point_matrix)/sqrt(4096);
Rx_fre_matrix = fft(Rx_eff_point_matrix,size(Rx_eff_point_matrix, 1))/sqrt(4096);
% 取出接收端导频 819*170?
Rx_pilot = Rx_fre_matrix(loc,:);
Xp = Xp(1:Np);

%% -------------------------------计算H-------------------------------------
for i = 1:size(Rx_pilot,2)
    H_esti(:,i) = Rx_pilot(:,i)./Xp; % 819*170
end

%% ----------------------------画图，可视化信号-----------------------------
% ---看一下信号强度
abs_H = abs(H_esti);
abs_H_vector = reshape(abs_H, size(abs_H,1) * size(abs_H,2), 1);
x_len = 1:size(H_esti,1)*size(H_esti,2);
figure;
scatter(x_len, abs_H_vector);
title('H abs'); % 每组第一个值为什么最强：表示直流分量，0Hz处的响应。
xlabel('Subcarrier Index');
ylabel('signal power');
hold;

phase_angle = angle(H_esti);
% ------------------------用虚部/实部=phase，画全部有效loops
figure;
phase_angle_vector = reshape(phase_angle, size(phase_angle,1) * size(phase_angle,2), 1);
x_aix = 1:length(phase_angle_vector);
% 打分组标签
group_labels = ceil(x_aix / 819); % 生成分组标签
scatter(x_aix,phase_angle_vector,'.');
hold on;
set(gca, 'XTick', 819:819:length(phase_angle_vector));
set(gca, 'XTickLabel', unique(group_labels));
xlabel('Subcarrier Index');
ylabel('Phase');
grid on;


% -------------------------用虚部/实部=phase，画单个loop
figure;
num=41;
vector_angle = phase_angle(819*(num-1)+1 : 819*(num-1)+819);
x_aix_sigle = 1:length(vector_angle);
scatter(x_aix_sigle,vector_angle,'.');
ylabel('phase');
xlabel('Subcarrier Index');
ylim([-4, 4]);

提取时间域特征，存入excel。timeFeatureExtraction.m

% phase_angle: 819*170
% phase_angle_vector: 139230*1

% -----------------------sample extraction------------
% features_sample: 样本数*特征数
samples1 = phase_angle; % C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature
%samples1 = phase_angle_filtered; % C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\filteredH
%samples1 = H_esti;
pilot_per_ofdm = 819;
features = table;           %特征表
sample_number = 50;        %sample_number为样本个数
loops_eachgroup = 5;        

% ----------------------计算11个特征值
for i=1:1:sample_number-loops_eachgroup
    % for i=1:1
    % ----当五个ofdm一组
    % v = samples1(1:pilot_per_ofdm,(i-1)*5+1:(i-1)*5+1+4);
    % v = reshape(v, [], 1);

    % -----当一个一组
    % v = samples1(1:pilot_per_ofdm,i);

    % -----当连续的五个ofdm一组
    v = samples1(1:pilot_per_ofdm,i:i+4);
    v = reshape(v, [], 1);

    features.Amplitude(i) = mean(abs(H_esti(:,i)));           % 1信号强度
    features.Mean(i) = mean(v);                         %2平均值
    features.Std(i) = std(v);                           %3标准差
    features.Skewness(i) = skewness(v);                 %4偏度
    features.Kurtosis(i) = kurtosis(v);                 %5峭度
    features.max(i) = max(v);                           %6最大值
    features.min(i) = min(v);                           %7最小值
    features.Peak2Peak(i) = peak2peak(v);               %8峰峰值
    features.RMS(i) = rms(v);                           %9均方根
    features.CrestFactor(i) = max(v)/rms(v);            %10振幅因数
    features.ShapeFactor(i) = rms(v)/mean(abs(v));      %11波形因数
    features.ImpulseFactor(i) = max(v)/mean(abs(v));    %12冲击因数
    features.MarginFactor(i) = max(v)/mean(abs(v))^2;   %13裕度因数
    features.Energy(i) = sum(v.^2);                     %14能量
    [corr_value,lags] = xcorr(v(1:410),v(411:end));                      
    features.auto_corr(i) = max(corr_value);            %15自相关性
    if i ~= sample_number/loops_eachgroup
    this_5 = sum(v(815:819));
    next_5 = sum(i*pilot_per_ofdm+1:i*pilot_per_ofdm+5);
    end
    features.next_group_related(i) = this_5-next_5;        % 16和下一组之间的联系
end

% --------------保存features table到特定文件夹下，并重新命名
%folderName = 'feature_excel_table';
folderName = 'svdestimate';
fileName = '5duzouyuan3.xlsx';
filePath = fullfile(folderName, fileName);
writetable(features, filePath);

提取频率域特征，存入excel。HstftFeatureExtraction.m

% phase_angle: 819*170
% phase_angle_vector: 139230*1

% -----------------------sample extraction------------
% features_sample: 样本数*特征数
% phase没有频率特性，这里取H的fft
samples2 = H_esti;
pilot_per_ofdm = 819;
fre_features = table;           %特征表
sample_number = 50;          %sample_number为样本个数
loops_eachgroup = 5;
frequency_samples= zeros(floor(pilot_per_ofdm/2),sample_number);

 %---------------------把时域每个样本都从时域通过傅里叶变换到频域--------------
for i=1:1:sample_number     
    x = samples2(1:pilot_per_ofdm,i);
    L = length(x);
    y = fft(x);
    y = y/L;
    P2 = abs(fft(x)/L);         % 每一个点的强度值
    P1 = P2(1:L/2);
    P1(2:end-1) = 2*P1(2:end-1);
    frequency_samples(:,i)= P1; %P1为向量，其长度为P1_length
                                %frequency_samples为所有样本的频域数据幅值，同样为矩阵
                                %行数为P1_length，列大小即为样本个数
end

% ----------------------------提取频域相关特征-------------------------------
for i=1:1:sample_number-loops_eachgroup  
    fre_v = frequency_samples(:,i:i+4);
    fre_v = reshape(fre_v, [], 1);
    fre_features.Mean(i) = mean(fre_v);                         %平均值
    fre_features.Std(i) = std(fre_v);                           %标准差
    fre_features.Skewness(i) = skewness(fre_v);                 %偏度
    fre_features.Kurtosis(i) = kurtosis(fre_v);                 %峭度
    fre_features.max(i) = max(fre_v);                           %最大值
    fre_features.min(i) = min(fre_v);                           %最小值
    fre_features.Peak2Peak(i) = peak2peak(fre_v);               %峰峰值
    fre_features.RMS(i) = rms(fre_v);                           %均方根
    fre_features.CrestFactor(i) = max(fre_v)/rms(fre_v);            %振幅因数
    fre_features.ShapeFactor(i) = rms(fre_v)/mean(abs(fre_v));      %波形因数
    fre_features.ImpulseFactor(i) = max(fre_v)/mean(abs(fre_v));    %冲击因数
    fre_features.MarginFactor(i) = max(fre_v)/mean(abs(fre_v))^2;   %裕度因数
    fre_features.Energy(i) = sum(fre_v.^2);                     %能量
end
% --------------保存features table到特定文件夹下，并重新命名
folderName = 'HandPhaseFeature\Hfre';
fileName = '5duzouyuan3_conti_5.xlsx';
filePath = fullfile(folderName, fileName);
writetable(fre_features, filePath);

从excel中读取全部特征，准备投入分类器。 newInputGenerator.m

% 5liandun1
du5liandun1_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\\HandPhaseFeature\5liandun1.xlsx');
du5liandun1_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\\HandPhaseFeature\Hfre\only5frefeatures\5liandun1.xlsx');

% 5liandun2
du5liandun2_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5liandun2.xlsx');
du5liandun2_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5liandun2.xlsx');

% 5liandun3
du5liandun3_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5liandun3.xlsx');
du5liandun3_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5liandun3.xlsx');

% 5ss1
du5ss1_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5ss1.xlsx');
du5ss1_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5ss1.xlsx');

% 5ss2
du5ss2_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5ss2.xlsx');
du5ss2_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5ss2.xlsx');

% 5ss3
du5ss3_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5ss3.xlsx');
du5ss3_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5ss3.xlsx');

% 5度 走远1
du5zouyuan1_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5duzouyuan1.xlsx');
du5zouyuan1_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5duzouyuan1.xlsx');

% 5度 走远2
du5zouyuan2_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5duzouyuan2.xlsx');
du5zouyuan2_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5duzouyuan2.xlsx');

% 5度 走远3
du5zouyuan3_time = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\5duzouyuan3.xlsx');
du5zouyuan3_fre = xlsread('C:\Users\Z\Desktop\ISAC\code\phwyh\yh\HandPhaseFeature\Hfre\only5frefeatures\5duzouyuan3.xlsx');

%------------------把每个动作的，Phase时域，H频域特征放在一起
all5du_ss1 = [du5ss1_time du5ss1_fre];
all5du_ss2 = [du5ss2_time du5ss2_fre];
all5du_ss3 = [du5ss3_time du5ss3_fre];
all5du_liandun1 = [du5liandun1_time du5liandun1_fre];
all5du_liandun2 = [du5liandun2_time du5liandun2_fre];
all5du_liandun3 = [du5liandun3_time du5liandun3_fre];
all5du_zouyuan1 = [du5zouyuan1_time du5zouyuan1_fre];
all5du_zouyuan2 = [du5zouyuan2_time du5zouyuan2_fre];
all5du_zouyuan3 = [du5zouyuan3_time du5zouyuan3_fre];
%------------------分一下 5度 走远,静止和liandun
all5du_jingzhi = [all5du_ss1;all5du_ss2;all5du_ss3];
all5du_zouyuan = [all5du_zouyuan1;all5du_zouyuan2;all5du_zouyuan3];
all5du_liandun = [all5du_liandun1;all5du_liandun2;all5du_liandun3];
ft_all = [all5du_jingzhi;all5du_zouyuan;all5du_liandun];
y_label = [zeros(135, 1);2 * ones(135, 1); ones(135, 1)];
x_input = [ft_all,y_label];

matlab分类器的使用：

toolbox：Statistics and Machine Learning Toolbox