60 GHz FMCW 雷达手势数据集及MATLAB仿真代码

clc;close all;clear ;%% 作者：调皮连续波，微信公众号%% 时间：20241101%% 功能：实现FMCW MIMO雷达雷达手势数据读取和处理，输出距离、速度和角度信息，并用于手势识别%% 数据集：https://ieee-dataport.org/documents/60-ghz-fmcw-radar-gesture-dataset
%% 添加MATLAB读取npy数据的函数路径addpath('..\fulldata_zipped\npy-matlab') ;addpath('..\fulldata_zipped\fulldata');
%% 读取npy数据filename ='SwipeUp_e4_u1_s394.npy';your_data = readNPY(filename);
%% 第一维表示每个手势的帧长，第二维表示虚拟天线的数量，第三维表示一帧中的chirps数量，第四维表示样本数量frame =size(your_data,1);VirtualArray =  size(your_data,2);chirps =size(your_data,3);AdcSamples = size(your_data,4);
%% 波形参数B = ( 62.5- 58.5) *1e9;c =3*10^8;f = 60*10^9;lambda = c/f;res =c/(2*B);FrameRate =33;Tf = 1/FrameRate;Tc = 300*10^-6;vres = lambda/(2*chirps*Tc);fftpoint =AdcSamples/2;
%% 坐标轴Rax = (0:fftpoint-1)*res;Vex = (-chirps/2:chirps/2-1)*vres;
Rangefft = complex(zeros(VirtualArray,chirps,AdcSamples/2));Doppler = complex(zeros(VirtualArray,chirps,AdcSamples/2));
RangefftMTI = complex(zeros(VirtualArray,chirps,AdcSamples/2));

MtiSize =6; %MTI滑动窗长度MtiBuffer = complex(zeros(VirtualArray,AdcSamples/2,MtiSize)); 
%%gifName = '..\fulldata_zipped\fulldata\test.gif';   % 打开或者创建一个gif文件% 打开GIF文件以写入  if exist(gifName, 'file')      delete(gifName); % 如果文件已存在，则删除  end  fid = fopen(gifName, 'w');  
for i =1:frame    RadarData = squeeze(your_data(i,:,:,:));    for j =1:VirtualArray        %时域信号        adcDataRowRx = squeeze(RadarData(j,:,:));        %%RangeFFT        for k = 1:chirps            temp= fft(adcDataRowRx(k,:));            Rangefft(j,k,:)  =temp(:,1:AdcSamples/2);        end        %%静态杂波滤除算法        MtiBuffer(j,:,1:end -1) =MtiBuffer(j,:,2:end);        MtiBuffer(j,:,end) = mean(squeeze(Rangefft(j,:,:)),1);        RangefftMTI(j,:,:) = squeeze(Rangefft(j,:,:))-mean(squeeze(MtiBuffer(j,:,:)),2).';        %%DopplerFFT        for m =1:AdcSamples/2            Doppler(j,:,m) = fftshift(fft(squeeze(RangefftMTI(j,:,m))));        end    end
    %%三个通道非相参积累    DopplerNCC = squeeze(sum(abs(Doppler),1));
    % figure(1);    % subplot(121);    % mesh(abs(adcDataRowRx));    % xlabel('samples');ylabel('chirps');zlabel('Amp');title([filename,' # ADC data Frmes = ',num2str(i)]);grid on    %    % subplot(122);    % imagesc(10*log10(abs(squeeze(RangefftMTI(1,:,:)))));axis xy    % xlabel('Range(m)');ylabel('chirps');zlabel('Power(dB)');title('RangeFFT data');grid on

    %%caponBFC测水平角    scale = -50:1:50;   % 扫面范围    M=2;           % 阵元单元    d = 0:lambda/2:(M-1)*lambda/2;% 阵列天线    N = size(Doppler,2);    Rx(:,:,1)  =squeeze(Doppler(1,:,:));    Rx(:,:,2)  =squeeze(Doppler(3,:,:));    P_Capon =zeros(length(Rax),length(scale));    idR =0;    for r =1:fftpoint         idR =idR+1;        Y =squeeze( Rx(:,r,:)).';        R = (Y*Y')/N;        R_inv = pinv(R);        idx = 1;        for n=scale            steeringVector  = (exp(-1j*2*pi/lambda*d*sind(n))).'; % 转向矢量转置            P_Capon(idR,idx) =( 1/(steeringVector'*R_inv*steeringVector)); % Capon算法            idx = idx + 1;        end    end
    %%caponBF测俯仰角    scale = -50:1:50;   % 扫面范围    M=2;           % 阵元单元    d = 0:lambda/2:(M-1)*lambda/2;% 阵列天线    N = size(Doppler,2);    Rx(:,:,1)  =squeeze(Doppler(2,:,:));    Rx(:,:,2)  =squeeze(Doppler(3,:,:));    P_CaponEle =zeros(length(Rax),length(scale));    idR =0;    for r =1:fftpoint         idR =idR+1;        Y =squeeze( Rx(:,r,:)).';        R = (Y*Y')/N;        R_inv = pinv(R);        idx = 1;        for n=scale            steeringVector  = (exp(-1j*2*pi/lambda*d*sind(n))).'; % 转向矢量转置            P_CaponEle(idR,idx) =( 1/(steeringVector'*R_inv*steeringVector)); % Capon算法            idx = idx + 1;        end    end
    figure(2);    subplot(131);    imagesc(Rax,Vex,10*log10(DopplerNCC));axis xy    clim([30,40])    xlabel('Range(m)');ylabel('Doppler(m/s)');zlabel('Power(dB)');title(['DopplerFFT data # Frames = ',num2str(i)]);grid on
    subplot(132);    imagesc(scale,Rax,10*log10(abs(P_Capon))); axis xy    clim([50,65])    xlabel('Horizontal Angle(°)');ylabel('Range(m)');zlabel('Power(dB)');title('Horizontal CaponBF DOA');grid on
    subplot(133);    imagesc(scale,Rax,10*log10(abs(P_CaponEle))); axis xy    clim([50,65])    xlabel('Pitch Angle(°)');ylabel('Range(m)');zlabel('Power(dB)');title('Pitch CaponBF DOA');grid on
    % 写入gif过程    frame = getframe(gcf);    im{i} = frame2im(frame);    [A,map] = rgb2ind(im{i},256); %将 RGB 图像转换为索引图像 A。map 最多包含 n 个颜色。    if i == 1        imwrite(A,map,gifName,'gif','LoopCount',Inf,'DelayTime',0.1);%'DelayTime'为每帧图像播放时间    else        imwrite(A,map,gifName,'gif','WriteMode','append','DelayTime',0.1);    end    % pause(0.1)
endfclose(gifName);