使用Matlab接收Wifi Csi并实时分析呼吸速率_csi 实时数据处理

使用Matlab接收Wifi Csi并实时分析呼吸速率

作者当前研究的问题

在室内环境下，研究通过普通的Wifi信号评估人的呼吸速率。

最后实现的效果

左图是穿戴传感器测量的实时呼吸率，3秒更新一次，右图是实时分析Wifi的Channel State Information的评估结果。呼吸速率的评估存在一些误差，但本文侧重点是如何在Matlab上进行数据的实时处理。因为本文是作者目前的毕设研究内容，关键代码部分暂不公开，同时，由于作者并不了解Matlab相关工具的使用，有些代码显得很冗余，请各位见谅。

数据采集的设置

数据采集的设置环境十分重要，开始的几个月由于没有注意到这点，困扰了好长时间。在 “Human Respiration Detection with Commodity WiFi Devices: Do User Location and Body Orientation Matter”中，作者提到了呼吸信号可检测性的3个因素，可以概括较理想的检测条件是：
1、呼吸引起的胸腔位移发生在菲涅尔区的中央
2、位移方向垂直于LOS
3、利用Intel5300可测量的30条载波的频率分集实现频率选择性检测

最后的数据采集设置为：被检测的位移发生在垂直于LOS路径距离为120cm处。使用Csi tool的Monitor采集数据，利用接收天线对的相位差做数据分析。

原始的相位差时域波性。

Matlab的实时处理

作者是基于 “Linux 802.11n CSI Tool下csi数据的实时可视化-Realtime-processing-for-csitool”工具做的改进。改进点如下，

1、基于生产者-消费者模型，对连续采集的csi数据分析数据。
2、基于实时性要求，使用Matlab的本地多核任务并行进行编程。
3、将消费者产生的呼吸频率和强度通过tcpip发送给另一个Matlab运行实例进行显示。

通过暂停发送数据，接收端会在延迟1S停止绘图，可以验证工具的实时性。

Matlab程序

1、read_bf_socket.m
clear all
delete(gcp)
parpool(‘local’,4)

spmd
switch labindex
case 1
subcarrier=30;
fs=100;
csi_entry_BufferSize=24fs;
adsize=24fs;
transfersize=8*fs;

                timestamp=zeros(1,csi_entry_BufferSize);
                Phase_diff=zeros(csi_entry_BufferSize,subcarrier);
                csi_writepos=1;
                csi_readpos=1;
                csi_items=0;                     
                shift_time=3;
                shift_point=shift_time*fs;
                woker=2;
                
                index = -1;                     % The index of the plots which need shadowing
                broken_perm = 0;                % Flag marking whether we've encountered a broken CSI yet
                triangle = [1 3 6];             % What perm should sum to for 1,2,3 antennas
                csi_entry=[];
                
                
           %% Build a TCP Server and wait for connection
                port = 8090;
                t= tcpip('0.0.0.0', port, 'NetworkRole', 'server');
                t.InputBufferSize = 4096;
                t.Timeout = 15;
                fprintf('Waiting for connection on port %d\n',port);
                fopen(t);
                fprintf('Accept connection from %s\n',t.RemoteHost);  
                
                while 1
                % Read size and code from the received packets
                        s = warning('error', 'instrument:fread:unsuccessfulRead');
                        try
                            field_len = fread(t, 1, 'uint16');
                        catch
                            warning(s);
                            disp('Timeout, please restart the client and connect again.');
                            break;
                        end
                        
                        code= fread(t,1);    
                        % If unhandled code, skip (seek over) the record and continue
                        if (code == 187) % get beamforming or phy data
                            bytes = fread(t, field_len-1, 'uint8');
                            bytes = uint8(bytes);
                            if (length(bytes) ~= field_len-1)
                                fclose(t);
                            end
                        else if field_len <= t.InputBufferSize  % skip all other info
                            fread(t, field_len-1, 'uint8');
                            continue;
                            else
                                continue;
                            end
                        end

                if (code== 187) % (tips: 187 = hex2dec('bb')) Beamforming matrix -- output a record
                    csi_entry= read_bfee(bytes);      
                    perm = csi_entry.perm;
                    Nrx = csi_entry.Nrx;
                    if Nrx > 1 % No permuting needed for only 1 antenna
                        if sum(perm) ~= triangle(Nrx) % matrix does not contain default values
                            if broken_perm == 0
                                broken_perm = 1;
                                fprintf('WARN ONCE: Found CSI (%s) with Nrx=%d and invalid perm=[%s]\n', filename, Nrx, int2str(perm));
                            end
                        else
                            csi_entry.csi(:,perm(1:Nrx),:) = csi_entry.csi(:,1:Nrx,:);
                        end
                    end
                end
                if code==187
                %% Initialize variable
                    timestamp(csi_writepos)=csi_entry.timestamp_low; 
                    csi = get_scaled_csi(csi_entry);    
                    csi1=squeeze(csi).';
                    csi2=angle(csi1)';   
                    Phase_diff (csi_writepos,1)=csi2(1,:)-csi2(3,:);       
                    csi_writepos=csi_writepos+1;
                    if(csi_writepos>adsize)
                         csi_writepos=1; 
                    end  
                    csi_items=csi_items+1;

                    if csi_items>=transfersize+transfersize+shift_point
                        end_pos=csi_readpos+transfersize+shift_point;
                        if end_pos>adsize
                               time_tmp=[timestamp(csi_readpos:adsize),timestamp(1:end_pos-adsize)];
                               pd_tmp=[Phase_diff((csi_readpos:adsize),:)',Phase_diff((1:end_pos-adsize),:)']';                                          
                        else
                               time_tmp=timestamp(csi_readpos:end_pos);
                               pd_tmp=Phase_diff((csi_readpos:end_pos),:);   
                        end                            s2=struct('timestamp',time_tmp(1:transfersize),'pd_data',pd_tmp((1:transfersize),:),'rate',fs,'sub_cnt',subcarrier);                            s3=struct('timestamp',time_tmp(shift_point+1:end),'pd_data',pd_tmp((shift_point+1:end),:),'rate',fs,'sub_cnt',subcarrier);   
                        labSend(s2,2);
                        labSend(s3,3);                             
                        csi_items=csi_items-2*shift_point;
                        csi_readpos=csi_readpos+2*shift_point;
                        if csi_readpos>adsize
                             csi_readpos=csi_readpos-adsize;
                        end
                    elseif csi_items>=transfersize
                           end_pos=csi_readpos+transfersize;
                                if end_pos>adsize
                                       time_tmp=[timestamp(csi_readpos:adsize),timestamp(1:end_pos-adsize)];
                                       pd_tmp=[Phase_diff((csi_readpos:adsize),:)',Phase_diff((1:end_pos-adsize),:)']';                                          
                                else
                                       time_tmp=timestamp(csi_readpos:end_pos);
                                       pd_tmp=Phase_diff((csi_readpos:end_pos),:);         
                                end                                   s2=struct('timestamp',time_tmp(1:transfersize),'pd_data',pd_tmp((1:transfersize),:),'rate',fs,'sub_cnt',subcarrier);
                                worker=5-worker;
                                labSend(s2,worker);                           
                                csi_items=csi_items-shift_point;
                                csi_readpos=csi_readpos+shift_point;
                                if csi_readpos>adsize
                                     csi_readpos=csi_readpos-adsize;
                                end       
                    end
                end
                end
        case 2
            while 1
                 cardio_process();    %数据的分析处理，此部分程序暂不公开
            end
        case 3
            while 1
                 cardio_process();
            end
        case 4   	
                t = tcpip('localhost', 30000, 'NetworkRole', 'client');
                fopen(t);
           while 1               
                  data=labReceive('any',10);
                  sender=[data.Breathrate  data.Breathalp data.Heartrate data.heartalp];
                  fwrite(t, sender,'double');
            end
           end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

end
matlabpool close

2、Plot.m
clf;
fs=100;
shift_time=3;
shift_point=shift_timefs;
% 绘图相关
x1=0;
y1=0;
frame_cnt=0;
figure1=figure(1);
ax1=subplot(2,1,1);
h=animatedline(x1,y1);
axis([x1,x1+10,-10,10]);
title(‘Breath’);
ylabel('amp10e3/fre*3’);
grid on;

ax2= subplot(2,1,2);
h1=animatedline(x1,y1);
axis([x1,x1+10,-10,10]);
title(‘Heartbeat’);
ylabel(‘amp10e3/fre3’);
xlabel(‘time/seconds’);
grid on;
drawnow
% Create textbox
dim1 = [0.00278571428571429 0.940476190476191 0.257928571428571 0.0547619047619057];
anno=annotation( figure1,‘textbox’,dim1,‘LineStyle’,‘none’,‘FitBoxToText’,‘on’);

dim3=[0.00854716981132076 0.946946946946947 0.0841752021563343 0.0180180180180181];
anno3=annotation( figure1,‘textbox’,dim3,‘LineStyle’,‘none’,‘FitBoxToText’,‘on’);

dim2= [0.00278571428571429 0.469047619047619 0.263285714285714 0.061904761904763];
anno2=annotation( figure1,‘textbox’,dim2,‘LineStyle’,‘none’,‘FitBoxToText’,‘on’);

dim4=[0.00854716981132076 0.485485485485485 0.0820188679245283 0.0190190190190191];
anno4=annotation( figure1,‘textbox’,dim4,‘LineStyle’,‘none’,‘FitBoxToText’,‘on’);

t = tcpip(‘0.0.0.0’, 30000, ‘NetworkRole’, ‘server’);

%Open a connection. This will not return until a connection is received.

fopen(t);

%绘制反映时域中呼吸及心率变化图
while 1
% if exist(auxfilename,‘file’)~=0 %文件存在
% ! rm 1.txt
% [f1,s1,f2,s2]=textread(filesname,‘%f %f %f %f’) ;
if get(t,‘BytesAvailable’)>1
data = fread(t, 4,‘double’);
f1=data(1);
s1=data(2);
f2=data(3);
s2=data(4);
tint=0:1/fs:shift_time;
tint1=tint+frame_cntshift_time;
%将呼吸和心率在画图的时候频率扩大了三倍
curv1=s110e2sin(2pi3f1tint); %呼吸信号的曲线,
curv2=s210e2sin(2pi3f2*tint); %心率信号的曲线

    str1=['Breathrate  ',num2str(f1*60),'bpm'];
1

% str2=['Ground-Breathrate ',num2str(100),‘bpm’];
str3=['Heartrate ',num2str(f2*60),‘bpm’];
% str4=['Ground- Heartrate ',num2str(100),‘bpm’];

    anno.String=str1;
1

% anno3.String=str2;
anno2.String=str3;
% anno4.String=str4;

    for k=1:length(tint1)
         addpoints(h,tint1(k),curv1(k));
         addpoints(h1,tint1(k),curv2(k));
         if(frame_cnt>=3)
             x1=x1+1/fs;
             axis([ax1,ax2],[x1,x1+10,-10,10]);
         end
         drawnow update
    end
    %saveas(figure1,['Plot' num2str(frame_cnt) '.jpg']);
    frame_cnt=frame_cnt+1;
1
2
3
4
5
6
7
8
9
10
11

end
end

3、cardio_process.m

4、完整的matlab代码如链接(资源需付费下载)： https://pan.baidu.com/s/1q4fXQgK4UesLXAfA2VkZ7g?pwd=9531
提取码：9531

视频

http://player.youku.com/embed/XNDM5MzczOTc1Ng==