1. AI技术聚合首页
  2. 技术文章

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

乘风 • • 技术文章 • 阅读 152

雷达回波的多普勒谱提取

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

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

数据集介绍

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

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

这个数据集有七个文件夹
在这里插入图片描述
代表作者在不同时间点采集的数据
数据集的数据不是特别多，本文主要过一下特征提取过程

特征提取

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

  • 距离压缩
      Dechirp体制雷达的距离压缩如上篇轨迹提取的文章，采用加窗FFT的方式，结果如下图所示在这里插入图片描述
      中间有个红线，不是特别清楚原因，可能是因为是三角波的缘故？作者只选取了底下比较强的一半数据进行处理。

  • MTI
      上一篇文章中MTI使用了两脉冲对消，对于该数据集的处理，原作者使用了一个巴特沃斯的高通滤波器，下面是处理后的结果前后对比
    【FMCW雷达人体行为识别——多普勒谱提取】【FMCW雷达人体行为识别——多普勒谱提取】

  • 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 for 2nd 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([-6 6]); colorbar
colormap; %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: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变换、小波变换等。相关知识可以在时频信号分析相关书籍中找到。下面再给出一种同样比较简单有效的方法：流程图如下：
【FMCW雷达人体行为识别——多普勒谱提取】
下面给出示例代码：

 %%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就是最后的结果

文章出处登录后可见！

已经登录？立即刷新

共计人评分，平均分

到目前为止还没有投票！成为第一位评论此文章。

赞 (0)
乘风的头像乘风管理团队
0
上一篇 2023年2月25日 下午9:14
下一篇 2023年2月25日 下午9:15

相关推荐