雷达回波的多普勒谱提取

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

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

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

数据集介绍

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

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

这个数据集有七个文件夹

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

特征提取

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

  • 距离压缩
     Dechirp体制雷达的距离压缩如上篇
     中间有个红线，不是特别清楚原因，可能是因为是三角波的缘故？作者只选取了底下比较强的一半数据进行处理。

  • MTI
     上一篇文章中MTI使用了两脉冲对消，对于该数据集的处理，原作者使用了一个巴特沃斯的高通滤波器，下面是处理后的结果前后对比

  • STFT
     通过短时傅里叶变换提取目标的多普勒信息，由于目标大致运动范围在10到30个距离单元之间，所以STFT采用这一区间内的数据进行。短时傅里叶变换的示意图如下

    可以看出,对一个距离bin进行多次STFT后得到一个横向为时间，纵向为多普勒信息的多普勒谱。对于多个距离bin，采用非相干叠加的方式得到最终的多普勒谱图如下：

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

%% Data reading partclearclose 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 frequencyTsweep = radarData(2); % Sweep time in msTsweep=Tsweep/1000; %then in secNTS = radarData(3); % Number of time samples per sweepBw = 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 formatfs=NTS/Tsweep; % sampling frequency ADCrecord_length=length(Data)/NTS*Tsweep; % length of recording in snc=record_length/Tsweep; % number of chirps%% Reshape data into chirps and plot Range-TimeData_time=reshape(Data, [NTS nc]);win = ones(NTS,size(Data_time,2));%Part taken from Ancortek code for FFT and IIR filteringtmp = 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));endfreq =(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 for 2nd FFT to get Doppler% This selects the range bins where we want to calculate the spectrogrambin_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);endMD.TimeAxis=linspace(0,MD.WholeDuration,size(Data_spec_MTI2,2));Data_spec_MTI2=flipud(Data_spec_MTI2);figureimagesc(MD.TimeAxis,MD.DopplerAxis.*3e8/2/5.8e9,20*log10(abs(Data_spec_MTI2))); colormap('jet'); axis xyylim([-6 6]); colorbarcolormap; %xlim([1 9])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',[0 0 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:lengthFilesfilename=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就是最后的结果