Home > data-import > read_pos_sensor.m

read_pos_sensor

PURPOSE ^

read_pos_sensor Reads position sensor data, then puts it in timeseries

SYNOPSIS ^

function [s_pos,s_orient] = read_pos_sensor(sid, sensor_id_list, new_freq)

DESCRIPTION ^

 read_pos_sensor   Reads position sensor data, then puts it in timeseries 
                   form for easy manip.  also resamples to 60 Hz or 
                   new_freq Hz if specified. The function will try to read
                   data from 'position_sensor_r/chopped.txt'. It this file
                   does not exist, it will try the original data file
                   'all-*.txt'.
 
 [s_pos s_orient] = read_pos_sensor(sid, sensor_id_list,new_freq)
 INPUT
     sid:              subject id.
     sensor_id_list:   list of sensor ID. Default is all sensor, and empty [] means all sensors also.
     new_freq:         resampling frequency of data. Default is 60 Hz, and empty [] means 60 Hz also.
     adjust_para:      this parameter is used for unwrap orientation data
                       for moving event detection only.

 OUTPUT
     s_pos:            cell array; s_pos{i}: postions of sensor i. if i is
                       not a member of sensor_id_list, s_pos{i} will be empty.
     s_orient:         cell array; s_orient{i}: orientation of sensor i. if i is
                       not a member of sensor_id_list, s_orient{i} will be empty.


EXAMPLE 1
   [s_pos s_ori] = read_pos_sensor(1701);
       read all sensors' data and resample them in 60HZ

EXAMPLE 2
   [s_pos s_ori] = read_pos_sensor(1701, [1 3 5]);
       read data of sensor 1, 3, and 5, and resample them in 60HZ

%EXAMPLE 3
   [s_pos s_ori] = read_pos_sensor(1701, [ ], 30);
       read data of all sensors, and resample them in 30HZ

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [s_pos,s_orient] = read_pos_sensor(sid, sensor_id_list, new_freq)
0002 % read_pos_sensor   Reads position sensor data, then puts it in timeseries
0003 %                   form for easy manip.  also resamples to 60 Hz or
0004 %                   new_freq Hz if specified. The function will try to read
0005 %                   data from 'position_sensor_r/chopped.txt'. It this file
0006 %                   does not exist, it will try the original data file
0007 %                   'all-*.txt'.
0008 %
0009 % [s_pos s_orient] = read_pos_sensor(sid, sensor_id_list,new_freq)
0010 % INPUT
0011 %     sid:              subject id.
0012 %     sensor_id_list:   list of sensor ID. Default is all sensor, and empty [] means all sensors also.
0013 %     new_freq:         resampling frequency of data. Default is 60 Hz, and empty [] means 60 Hz also.
0014 %     adjust_para:      this parameter is used for unwrap orientation data
0015 %                       for moving event detection only.
0016 %
0017 % OUTPUT
0018 %     s_pos:            cell array; s_pos{i}: postions of sensor i. if i is
0019 %                       not a member of sensor_id_list, s_pos{i} will be empty.
0020 %     s_orient:         cell array; s_orient{i}: orientation of sensor i. if i is
0021 %                       not a member of sensor_id_list, s_orient{i} will be empty.
0022 %
0023 %
0024 %EXAMPLE 1
0025 %   [s_pos s_ori] = read_pos_sensor(1701);
0026 %       read all sensors' data and resample them in 60HZ
0027 %
0028 %EXAMPLE 2
0029 %   [s_pos s_ori] = read_pos_sensor(1701, [1 3 5]);
0030 %       read data of sensor 1, 3, and 5, and resample them in 60HZ
0031 %
0032 %%EXAMPLE 3
0033 %   [s_pos s_ori] = read_pos_sensor(1701, [ ], 30);
0034 %       read data of all sensors, and resample them in 30HZ
0035 %
0036 %
0037 
0038     if nargin < 1 || isempty(sid)
0039         help read_pos_sensor;
0040         error('You must provide the subject ID. Read help please!')
0041     end
0042 
0043     if nargin < 2 
0044         sensor_id_list = [];
0045     end
0046 
0047     if nargin < 3 || isempty( new_freq )     %change 2 to 3 by yiwen, because new_freq is the third parameter
0048         new_freq = 60;
0049     end
0050     s_pos = []; s_orient = [];
0051 
0052     pos_file = fullfile(get_subject_dir(sid), 'position_sensor_r/chopped.txt');
0053     if exist(pos_file, 'file')
0054         very_raw_pos = load(pos_file);
0055     else
0056         wildcard = fullfile(get_subject_dir(sid),'position_sensor_r','all*.txt');
0057         dir_entry = dir(wildcard);
0058         if ~isempty(dir_entry)
0059         pos_file = fullfile(get_subject_dir(sid),'position_sensor_r',dir_entry.name);
0060         very_raw_pos = dlmread(pos_file, '', 2, 0);
0061         disp('Skipping first two lines of position sensor file');
0062         else
0063             fprintf('Sensor .txt file does not exist for subject %d\n', sid);
0064             return
0065         end
0066     end
0067 
0068     timing_info = get_timing(sid);
0069     % offset = desired start time - actual start time
0070     % so timestamp + offset = timestamp - start time + desired start time.
0071     offset = timing_info.motionTime - very_raw_pos(1, 9);
0072 
0073     %to check the sensor_id_list
0074     all_sensor_id_list = unique( very_raw_pos(:,1) )';
0075     if isempty(sensor_id_list)
0076         sensor_id_list = all_sensor_id_list;
0077     end
0078     not_existed_sensor = setdiff(sensor_id_list, all_sensor_id_list);
0079     if ~isempty( not_existed_sensor )
0080         disp('The following sensor can not be found in data file!');
0081         disp( not_existed_sensor );
0082     end
0083 
0084     sensor_num = max(sensor_id_list);
0085 
0086     pos_out    = cell(sensor_num,1); 
0087     orient_out = cell(sensor_num,1); 
0088     for sensor = sensor_id_list
0089 
0090         raw_pos = very_raw_pos(very_raw_pos(:, 1) == sensor, :);
0091         raw_orient = raw_pos(:, 5:7);
0092 %         check_orientation_data(raw_orient);
0093 %         if nargin < 4 || isempty( adjust_para )
0094             unwrapped_orient = unwrap(raw_orient * (pi / 180)) * (180 / pi); 
0095 %         else
0096 %             %This is for run_moving_event_detection. I don't know why it use a
0097 %             %different parameter. Now, I think I have fixed this problem.
0098 %             unwrapped_orient = unwrap(raw_orient * (pi / 180)) * adjust_para;
0099 %         end
0100 
0101         pos = timeseries(raw_pos(:, 2:4), raw_pos(:, 9) + offset, 'Name', 'position');
0102         orient = timeseries(unwrapped_orient, raw_pos(:, 9) + offset, 'Name', 'orient');
0103 
0104         timeProps = get(pos, 'TimeInfo');
0105         new_basis = timeProps.Start:1/new_freq:timeProps.End;
0106 
0107     %    disp 'running linear interp'
0108         pos_lin = resample(pos, new_basis);
0109         orient_lin = resample(orient, new_basis);
0110 
0111 
0112     %     pchip_interp = @(new_Time,Time,Data)...
0113     %                interp1(Time,Data,new_Time,...
0114     %                        'pchip');
0115     %     pos = setinterpmethod(pos, pchip_interp);
0116     %     orient = setinterpmethod(orient, pchip_interp);
0117     %
0118     %     disp 'running pchip interp'
0119     %     pos_pchip = resample(pos, new_basis);
0120     %     orient_pchip = resample(orient, new_basis);
0121     %     disp 'done'
0122     %     hold off
0123     %     plot(pos, ':')
0124     %     hold on
0125     %     plot(pos_lin, '--')
0126     %     plot(pos_pchip, '-.')
0127 
0128 
0129         pos_out{sensor} = pos_lin;
0130         orient_out{sensor} = orient_lin;    
0131     end
0132 
0133     s_pos = cell(sensor_num,1);
0134     s_orient = cell(sensor_num, 1);
0135     for id = sensor_id_list
0136         s_pos{id}    = ts2cont(pos_out{id});
0137         s_pos{id} = filter_convert(s_pos{id}, new_freq);
0138         s_orient{id} = ts2cont(orient_out{id});
0139     end
0140 end
0141   
0142 
0143 function check_orientation_data(orient) 
0144     for i = 1 : size(orient,1)-1
0145         isDisp = 0;
0146         for j = 1 : size(orient,2)
0147             if orient(i+1,j)-orient(i,j) > 180 || orient(i+1,j)-orient(i,j)<-180
0148                 isDisp = 1;
0149             end
0150         end
0151         if isDisp
0152             disp(orient(i:i+1,:));
0153         end
0154     end
0155 end
0156 
0157 function result = filter_convert(data, freq)
0158 SF = freq; % sampling frequency
0159 NF = SF/2; % Nyquist frequency
0160 CF = 12; % Cut-off frequency
0161 % initalize normalized cut-off frequency Wn with a value between 0 and 1
0162 Wn = CF/NF; % == 9Hz/30Hz = 0.3
0163 % run butter
0164 [b,a] = butter(2, Wn, 'low'); % 2nd order, low-pass filter
0165 
0166 for ind = 2:4
0167     data(:,ind) = filtfilt(b,a,data(:,ind)); % b == numerator coefficients of filter, a == denominator coefficients, x == input data
0168 end
0169 data(:, 2:end) = data(:,2:end)*25.4;
0170 result = data;
0171 end
0172 
0173 % % The following two functions have been added to development lib
0174 % seperately
0175 %
0176 % function ts = cont2ts(cont)
0177 % ts = timeseries(cont(:, 2:end), cont(:, 1));
0178 % end
0179 %
0180 % function cont = ts2cont(ts)
0181 % cont = horzcat(get(ts, 'Time'), get(ts, 'Data'));
0182 % end
0183

Generated on Tue 23-May-2017 03:00:58 by m2html © 2005