【FMCW雷达人体行为识别——多普勒谱提取】

雷达回波的多普勒谱提取

  之前写过一个基于FMCW雷达的目标轨迹的提取，感觉看的人还是蛮多的，这周准备写一下关于多普勒谱提取的相关内容。主要内容为英国格拉斯哥大学公开的一个人体行为的数据集。数据集以及示例代码可以访问下方链接，如果访问不了可以下载如下压缩包获取。

• 官网地址：https://researchdata.gla.ac.uk/848/
• 压缩包：https://pan.baidu.com/s/1rW0OfuUrYc7kC9NZ6HHClQ 提取码：kw9d

  当然，官网也给出了提取多普勒谱的示例代码，下面将结合代码进行分析，最终构建图片数据集用于后续的识别分类。

数据集介绍

数据集采集自C波段(5.8GHz)的FMCW雷达,带宽为400MHz,调频周期为1ms。
人体行为的种类包括：行走、坐下、起立、俯身捡东西、喝水、跌倒。
数据集构建的思路包括：不同行为、左/右手、男/女性、身高/年龄、重复次数。

作者给出了不同数据获取的详细信息，但是吧，我觉得雷达信号对这些信息敏感度可能不是很高
同一个动作不同的人、年龄、身高甚至不同的手做差别有很大吗，这应该属于雷达中的微动领域吧

这个数据集有七个文件夹

代表作者在不同时间点采集的数据
数据集的数据不是特别多，本文主要过一下特征提取过程

特征提取

  流程比较简单，主要包括距离压缩、MTI、STFT。下面以一个行走的数据举例。

• 距离压缩   Dechirp体制雷达的距离压缩如上篇轨迹提取的文章，采用加窗FFT的方式，结果如下图所示   中间有个红线，不是特别清楚原因，可能是因为是三角波的缘故？作者只选取了底下比较强的一半数据进行处理。
• MTI   上一篇文章中MTI使用了两脉冲对消，对于该数据集的处理，原作者使用了一个巴特沃斯的高通滤波器，下面是处理后的结果前后对比
• STFT   通过短时傅里叶变换提取目标的多普勒信息，由于目标大致运动范围在10到30个距离单元之间，所以STFT采用这一区间内的数据进行。短时傅里叶变换的示意图如下   可以看出,对一个距离bin进行多次STFT后得到一个横向为时间，纵向为多普勒信息的多普勒谱。对于多个距离bin，采用非相干叠加的方式得到最终的多普勒谱图如下：

下面是上面三个步骤的代码：

%% Data reading part
clear
close all
[filename,pathname]=uigetfile('*.dat');
file=[pathname,filename];%%网站下载的程序没有这句话，所以源程序不可以选择任意文件夹下的图片
fileID =fopen(file,'r');
dataArray =textscan(fileID,'%f');fclose(fileID);
radarData = dataArray{1};
clearvars fileID dataArray ans;
fc =radarData(1);% Center frequency
Tsweep =radarData(2);% Sweep time in ms
Tsweep=Tsweep/1000;%then in sec
NTS =radarData(3);% Number of time samples per sweep
Bw =radarData(4);% FMCW Bandwidth. For FSK, it is frequency step;% For CW, it is 0.
Data =radarData(5:end);% raw data in I+j*Q format
fs=NTS/Tsweep;% sampling frequency ADC
record_length=length(Data)/NTS*Tsweep;% length of recording in s
nc=record_length/Tsweep;% number of chirps
%% Reshape data into chirps and plot Range-Time
Data_time=reshape(Data,[NTS nc]);
win =ones(NTS,size(Data_time,2));%Part taken from Ancortek code for FFT and IIR filtering
tmp =fftshift(fft(Data_time.*win),1);
Data_range=zeros(size(tmp,1)/2,size(tmp,2));%% ！！原始代码中无此句，此句说明见本文最后部分
Data_range(1:NTS/2,:)=tmp(NTS/2+1:NTS,:);
ns =oddnumber(size(Data_range,2))-1;
Data_range_MTI =zeros(size(Data_range,1),ns);[b,a]=butter(4,0.0075,'high');[h, f1]=freqz(b, a, ns);for k=1:size(Data_range,1)Data_range_MTI(k,1:ns)=filter(b,a,Data_range(k,1:ns));
end
freq =(0:ns-1)*fs/(2*ns); 
range_axis=(freq*3e8*Tsweep)/(2*Bw);
Data_range_MTI=Data_range_MTI(2:size(Data_range_MTI,1),:);
Data_range=Data_range(2:size(Data_range,1),:);%% Spectrogram processing for2nd FFT to get Doppler
% This selects the range bins where we want to calculate the spectrogram
bin_indl =10;
bin_indu =30;
MD.PRF=1/Tsweep;
MD.TimeWindowLength =200;
MD.OverlapFactor =0.95;
MD.OverlapLength =round(MD.TimeWindowLength*MD.OverlapFactor);
MD.Pad_Factor =4;
MD.FFTPoints = MD.Pad_Factor*MD.TimeWindowLength;
MD.DopplerBin=MD.PRF/(MD.FFTPoints);
MD.DopplerAxis=-MD.PRF/2:MD.DopplerBin:MD.PRF/2-MD.DopplerBin;
MD.WholeDuration=size(Data_range_MTI,2)/MD.PRF;
MD.NumSegments=floor((size(Data_range_MTI,2)-MD.TimeWindowLength)/floor(MD.TimeWindowLength*(1-MD.OverlapFactor)));
    
Data_spec_MTI2=0;
Data_spec2=0;for RBin=bin_indl:1:bin_indu
    Data_MTI_temp =fftshift(spectrogram(Data_range_MTI(RBin,:),MD.TimeWindowLength,MD.OverlapLength,MD.FFTPoints),1);
    Data_spec_MTI2=Data_spec_MTI2+abs(Data_MTI_temp);                                
    Data_temp =fftshift(spectrogram(Data_range(RBin,:),MD.TimeWindowLength,MD.OverlapLength,MD.FFTPoints),1);
    Data_spec2=Data_spec2+abs(Data_temp);
end
MD.TimeAxis=linspace(0,MD.WholeDuration,size(Data_spec_MTI2,2));

Data_spec_MTI2=flipud(Data_spec_MTI2);

figure
imagesc(MD.TimeAxis,MD.DopplerAxis.*3e8/2/5.8e9,20*log10(abs(Data_spec_MTI2)));colormap('jet'); axis xy
ylim([-66]); colorbar
colormap;%xlim([19])
clim =get(gca,'CLim');set(gca,'CLim',clim(2)+[-40,0]);xlabel('Time[s]','FontSize',16);ylabel('Velocity [m/s]','FontSize',16)set(gca,'FontSize',16)title(filename)

你以为本文到这里就已经完了吗？？？
  哈哈哈哈，怎么可能，我们的目的是获得这个多普勒图的数据集用来后续的识别分类，所以还有两个实际的问题:

图片的保存：保存上图经过色彩风格调整后的伪彩色图片而不是黑白图片，因为黑白图片看起来效果不是很好。
#黑白图片可以使用imshow显示，不过需要归一化；保存使用imwrite
批量化处理：原始的数据集是dat格式的数据而并非图片，所以需要批量生成图片

Solvation

1. 图片保存
  如果想保存图窗显示的伪彩色图，需要做的操作如下：

    去除坐标轴，这没啥好解释的
    去除白边，一般保存图窗都会带白边，这是我们不希望看到的
    调整尺寸，最后保存的图片尺寸会不一样，只能找比例关系再调整一下了

代码如下：

imagesc(out);colormap('jet');
clim =get(gca,'CLim');set(gca,'CLim',clim(2)+[-40,0]);[r,c]=size(out);set(gca,'xtick',[],'ytick',[]);%去除坐标轴
set(gca,'LooseInset',get(gca,'TightInset'))%去除白边
set(gcf,'innerposition',[00 c*16/25 r*16/25])%调整尺寸
savepath=['./img/1/',filename(1:end-4),'.jpg'];%保存地址根据文件名变化
saveas(gca,savepath,'jpg')%图窗保存

2.批量化处理
  直接读取文件夹到一个结构体，再使用for循环遍历文件夹体中的文件，然后保存的图片根据文件名变化就可以实现批量化处理了。部分关键代码如下：

Files=dir(fullfile('.\path\*.dat'));% path是包含dat文件的目录名
lengthFiles=length(Files);for i=1:lengthFiles
    filename=Files(i).name;
    pathname=Files(i).folder;
    file=[pathname,'/',filename];
    fileID =fopen(file,'r');
    dataArray =textscan(fileID,'%f');fclose(fileID);
    radarData = dataArray{1};---
    savepath=['./img/1/',filename(1:end-4),'.jpg'];%保存地址根据文件名变化
    saveas(gca,savepath,'jpg')%图窗保存
end

#这里补充一个小细节： 由于原始代码中没有对一个矩阵初始化，所以在批量化处理时由于数据集中有的数据尺寸不一样会导致出错。所以添加代码段中的 “ ！！”部分。

  最后展示得到的部分图片数据集结果：

  OK，本文的内容就到这里了。主要是基于原作者提供代码上进行的一些小tricks，供大家参考哇！！！
  最后的图片数据集我就不给出了，有需要的跑一下代码就可以了。😄😄哦,对了，最后补充一下，关于多普勒谱的构建方法还有很多：比如WVD变换、小波变换等。相关知识可以在时频信号分析相关书籍中找到。下面再给出一种同样比较简单有效的方法：流程图如下：

下面给出示例代码：

%%range_profile代表一系列距离像，n是距离采样点数，a是帧数，假设64帧，每帧有128个脉冲
 %下面得到64帧距离多普勒图
 for a=1:Z
      for n=1:N
      temp(n,:)=range_profile(n,128*(a-1)+1:128*a).*(doppler_win)';%把一系列脉冲按帧处理加窗
      temp_fft(n,:)=fftshift(fft(temp(n,:),Nc));%然后横向进行FFT，可以得到距离多普勒图，Nc是128
      end
      temp3d(:,:,a)=temp_fft(:,:);%64帧RDM
 end
%下面对每帧多普勒图中的所有距离单元加一块然后转置，拼接到V变量中
 for a=1:Z
     V(:,a)=sum(temp3d(:,:,a))';
 end
 %V就是最后的结果