|
目录
一、引言
二、毫米波雷达检测呼吸、心跳基本原理
1.TI官方开发资料：
2.博主“调皮连续波”开源资料以及原理讲解：
三、 毫米波雷达提取呼吸、心跳信号Matlab算法处理
1.硬件平台： IWR6843ISKEVM+DCA1000EVM
2.mmavestudio参数设置：
配置说明：
算法流程简介：
（1） 预处理：
（2）粗略的人体定位：距离维FFT
（3）消除静态干扰算法【因为后面用了滑动平均去噪，故这里不做静态干扰算法处理】
(4)经典算法提取相位：相位反正切
（5）相位解缠绕
（6）相位差分
（7）脉冲噪声去除：滑动平均滤波
（8）带通滤波器输出呼吸信号：
带通滤波器的设计可以参考上一篇内容：MATLAB设计滤波器之新版filterDesigner使用_CoCo哥的博客-CSDN博客
（9）带通滤波器输出心跳信号：
（10）提取心跳波的包络线，归一化心跳波
四、参考资料
已上代码资料到csdn的资源区：
五 、总结
一、引言
非雷达科班出身，“入坑”毫米波雷达一年多，现在把入门毫米波雷达检测呼吸、心跳的传统与改进算法进行开源并归纳整理了相关的资料。欢迎交流以及专业人士的指正。
二、毫米波雷达检测呼吸、心跳基本原理
1.TI官方开发资料：
Vital Signs 68xx Users Guide (ti.com)
2.博主“调皮连续波”开源资料以及原理讲解：
干货 | IWR1642EVM呼吸心跳原始数据采集与仿真分析（含MATLAB代码和数据） (qq.com)
（1）线性调频连续波雷达基本原理（第1讲）： 干货-线性调频连续波雷达基本原理（第1讲）
（2）线性调频连续波雷达基本原理（第2讲）： 干货-线性调频连续波雷达基本原理（第2讲）
（3）线性调频连续波雷达基本原理（第3讲）： 干货-线性调频连续波雷达基本原理（第3讲）
3. 驾驶员生命体征检测原理视频说明（1642EVM，77GHZ）中文讲解：
毫米波雷达的应用无处不在- 1.4 对驾驶员心跳呼吸检测的应用-模拟与混合信号在线培训- 德州仪器（TI）官方视频课程培训
三、 毫米波雷达提取呼吸、心跳信号Matlab算法处理
以下matlab代码基于“调皮连续波”的77GHz IWR1642EVM的算法处理代码进行改进，同时参考了期刊论文文献： 《基于毫米波雷达的生命体征检测》（张兰春，顾海潮）重新设计了带通滤波器分离呼吸心跳信号，并且还参考了另一篇论文：《mmHRV_Contactless_Heart_Rate_Variability_Monitoring_Using_Millimeter-Wave_Radio》对提取的心跳信号进行估计包络的归一化处理。
1.硬件平台： IWR6843ISKEVM+DCA1000EVM
2.mmavestudio参数设置：
配置说明：
雷达配置：一帧2个chirp，一帧50ms，ADC采样点为200个，采样率为4Msps，数据处理时仅采用每帧的第一个chirp。 毫米波雷达发射起始频率：f0=60.25GHz 调频斜率：S=64.985MHz/us（~65e17）
有效带宽B=ts*S=3.24925GHz≈3.25GHZ 距离分辨率：ΔR=0.04615m
3.信号处理算法（matlab）
算法流程简介：
（1） 预处理：
%复采样
numChirps = fileSize/2/numADCSamples/numRX; %含有实部虚部，除以2
%共2048个chirps（1024帧*2个chirp）
LVDS = zeros(1, fileSize/2);
%combine real and imaginary part into complex data将实部虚部结合成复数
%read in file: 2I is followed by 2Q adcData数据组成:两个实部，接着是两个虚部
counter = 1;
for i=1:4:fileSize-1 %1T4R
LVDS(1,counter) = adcData(i) + sqrt(-1)*adcData(i+2); %复数形式
LVDS(1,counter+1) = adcData(i+1)+sqrt(-1)*adcData(i+3);
counter = counter + 2;
end
% create column for each chirp：每一列为chirp
LVDS = reshape(LVDS, numADCSamples*numRX, numChirps);
%each row is data from one chirp：每一行为chirp
LVDS = LVDS.';
end
%% 重组数据（4条接收天线的复数数据）
adcData = zeros(numRX,numChirps*numADCSamples);
for row = 1:numRX
for i = 1: numChirps
adcData(row, (i-1)*numADCSamples+1:i*numADCSamples) = LVDS(i, (row-1)*numADCSamples+1:row*numADCSamples);
end
end
%重组数据retVal：200*2048矩阵，每一列为一个chirp
retVal= reshape(adcData(1, :), numADCSamples, numChirps); %取第一个接收天线数据，数据存储方式为一个chirp一列
process_adc=zeros(numADCSamples,numChirps/2);%每帧中的两个chrip取第一个，200*1024
for nchirp = 1:2:numChirps %1T4R （1T1R）只处理单发单收的数据，并且只处理两个chrip取出的第一个
process_adc(:, (nchirp-1)/2+1) = retVal(:,nchirp);
end
（2）粗略的人体定位：距离维FFT
fft1d= zeros(numChirps/2,numADCSamples);
for chirp_fft=1:numChirps/2
fft1d(chirp_fft,:) = abs(fft((process_adc(:,chirp_fft))));%
end
figure;
mesh(X,Y,fft1d);
xlabel('距离(m)');
ylabel('脉冲chrip数');
zlabel('幅度');
title('距离维-1DFFT结果');
（3）消除静态干扰算法【因为后面用了滑动平均去噪，故这里不做静态干扰算法处理】
%% 消除静态干扰
% 使用相量均值相消算法（平均相消算法）:效果较差
%{
fft1d_avg = zeros(1024,256);
avg = sum(fft_data(:,:))/256; %参考接收脉冲
for chirp=1:1024
fft1d_avg(chirp,:) = fft_data(chirp,:)-avg;
end
fft_data=fft1d_avg;
%}
%
%MTI：动目标显示算法
%{
for ii=1:numChirps-1 % 滑动对消，少了一个脉冲
fft_data(ii,:) = fft_data(ii+1,:)-fft_data(ii,:);
end
%在这里增加了最后一个chirp的消除：（有待改进）
fft_data(numChirps,:) =fft_data(numChirps-1 ,:);
%}
%MTD：动目标检测
%{
for ii=1:256
avg = sum(fft_data(:,ii))/1024;
fft_data(:,ii) = fft_data(:,ii) - avg ;
end
%}
(4)经典算法提取相位：相位反正切
提取人体定位的位置的相位信号，即胸腔振动信号x(t）（△R）：
 
%% 提取相位 extract phase(相位反正切)
%实虚部分离（为了提取rangebin的相位）
real_data = real(fft_data);%实部
imag_data = imag(fft_data);%虚部
for i = 1:numChirps
for j = 1:RangFFT %对每一个range bin取相位 extract phase（弧度rad）
angle_fft(i,j) = atan2(imag_data(i, j),real_data(i, j));
end
end
% Range-bin tracking 找出能量最大的点，即人体的位置
for j = 1:RangFFT
if((j*detaR)0.5) % 限定了检测距离为0.5-2.5m
for i = 1:numChirps % 进行非相干积累
fft_data_last(j) = fft_data_last(j) + fft_data_abs(i,j);%通过FFT后的多普勒信号的幅值进行定位
end
if ( fft_data_last(j) > range_max)
range_max = fft_data_last(j);
max_num = j; %最大能量序列号（range bin）maxnum
end
end
end
%% 取出能量最大点的相位 extract phase from selected range bin
angle_fft_last = angle_fft(:,max_num);%1024个chrip的某列range bin的相位
% 提取相位信号（原始）
figure;
plot(angle_fft_last);
xlabel('时间/点数（N）：对应每个chrip');
ylabel('相位');
title('未展开相位信号');
phi=angle_fft_last;
（5）相位解缠绕
%% 进行相位解缠 phase unwrapping(手动解)，自动解可以采用MATLAB自带的函数unwrap()
%或称为相位解卷绕，由于相位值在 [ − π ,π ] 之间，而我们需要相位展开以获取实际的位移曲线，
% 因此每当连续值之间的相位差大于或者小于±π时，通过从相位中减去2π来获得相位展开。
% n = 1;
% for i = 2:numChirps
% diff = angle_fft_last(i) - angle_fft_last(i-1); %连续值之间的相位差
% if diff > pi
% angle_fft_last(i:end) = angle_fft_last(i:end) - 2*pi;
% n = n + 1;
% elseif diff < -pi
% angle_fft_last(i:end) = angle_fft_last(i:end) + 2*pi;
% end
% end
%相位解包方法2：
for i = 2:numChirps
diff = angle_fft_last(i) - angle_fft_last(i-1); %连续值之间的相位差
while((diff>pi) || (diff pi
angle_fft_last(i) = angle_fft_last(i) - 2*pi;
diff = angle_fft_last(i) - angle_fft_last(i-1); %连续值之间的相位差
elseif diff < -pi
angle_fft_last(i:end) = angle_fft_last(i:end) + 2*pi;
diff = angle_fft_last(i) - angle_fft_last(i-1); %连续值之间的相位差
end
end
end
%
figure;
plot(angle_fft_last);
xlabel('时间/点数（N）：对应每个chirp');
ylabel('相位');
title('相位解缠 phase unwrapping后的结果');
%unwrap函数：
% phi=unwrap(phi);
% figure;
% plot(phi);
（6）相位差分
%% phase difference 相位差分
%通过减去连续的相位值，对展开的相位执行相位差运算，
% 这将有利于 ：增强心跳信号并消除硬件接收机存在的相位漂移，抑制呼吸信号及其谐波
angle_fft_last2=zeros(1,numChirps);
% 方法：相位差分是通过不断将当前采样点展开相位与前一采样点做差实现
for i = 2:numChirps
angle_fft_last2(i) = angle_fft_last(i) - angle_fft_last(i-1);
end
figure;
plot(angle_fft_last2);
xlabel('点数（N）：对应每个chrip');
ylabel('相位（rad）');
title('相位差分后信号');
（7）脉冲噪声去除：滑动平均滤波
%% 脉冲噪声去除：滑动平均滤波（选取 0.25 s 的滑动窗口，窗口长度为5）
% 去除由于测试环境引起的脉冲噪声
phi=smoothdata(angle_fft_last2,'movmean',5);
figure;
plot(phi);
title('滑动平均滤波相位信号');
%对相位信号作FFT
N1=length(phi);
FS=20;
FFT = abs(fft(phi)); %--FFT 取模，幅度
f=(0:N1-1)*(FS/N1); %其中每点的频率
%傅里叶变换结果对称
figure;
plot(f(1:N1/8),FFT(1:N1/8)) %取前一部分放大观察
xlabel('频率（f/Hz）');
ylabel('幅度');
title('相位信号FFT ');
（8）带通滤波器输出呼吸信号：
带通滤波器的设计可以参考上一篇内容：MATLAB设计滤波器之新版filterDesigner使用_CoCo哥的博客-CSDN博客
%% IIR带通滤波 Bandpass Filter 0.1-0.5hz，输出呼吸信号
fs =20; %呼吸心跳信号的采样率
%设计IIR，4阶巴特沃斯带通滤波器
load('coe3.mat', 'breath_pass');
breath_data = filter(breath_pass,phi);
figure;
plot(breath_data);
xlabel('时间(s)');
xticklabels({'0','10','20','30','40','50','60'});
ylabel('幅度');
title('呼吸时域波形');
%% 谱估计 -FFT -Peak interval
N1=length(breath_data);
fshift = (-N1/2:N1/2-1)*(fs/N1); %
breath_fre = abs(fftshift(fft(breath_data))); %--FFT 取模，幅度（双边频谱）
breath = abs(fft(breath_data)); %
%傅里叶变换结果对称
figure;
% plot(fshift,breath_fre);
plot(f(1:130),breath(1:130)); %取前一部分放大观察
xlabel('频率（f/Hz）');
ylabel('幅度');
title('呼吸信号FFT ');
breath_fre_max = 0; % 呼吸频率
for i = 1:length(breath_fre) %谱峰最大值搜索
if (breath_fre(i) > breath_fre_max)
breath_fre_max = breath_fre(i);
breath_index=i;
end
end
breath_count =(fs*(numChirps/2-(breath_index-1))/numChirps)*60; %呼吸频率解算
（9）带通滤波器输出心跳信号：
%% IIR带通滤波 Bandpass Filter 0.8-2hz，输出心跳的数据
% 设计IIR，8阶巴特沃斯带通滤波器
load('coe4.mat', 'heart_pass');
heart_data = filter(heart_pass,phi);
figure;
plot(heart_data);
xlabel('时间(s)');
xticklabels({'0','10','20','30','40','50','60'});
ylabel('幅度');
title('心跳时域波形');
N1=length(heart_data);
fshift = (-N1/2:N1/2-1)*(fs/N1); % zero-centered frequency
f=(0:N1-1)*(fs/N1); %其中每点的频率
heart_fre = abs(fftshift(fft(heart_data)));
heart=abs(fft(heart_data));
figure;
% plot(fshift,heart_fre);
plot(f(1:200),heart(1:200));
xlabel('频率（f/Hz）');
ylabel('幅度');
title('心跳信号FFT');
heart_fre_max = 0;
for i = 1:length(heart_fre)/2
if (heart_fre(i) > heart_fre_max)
heart_fre_max = heart_fre(i);
if(heart_fre_max |