%  MATLAB example script using the PECS server
%  Macro writtern 10/12/2004 by Iain Charlton
%  Copyright (C) 2004 OMG Plc. All rights reserved

% Global matrices containing the peak knee flexion and the % gait cycle at
% which it occurs


% Get the PECS server and refresh it (NB refresh isn't required, it just
% updates the info on display)

global  hServer;
global  hTrial;
global  hProcessor;
global  hParamStore;
global  hEvStore;

hServer = actxserver( 'PECS.Document' );
invoke( hServer, 'Refresh' );

hTrial = get(hServer, 'Trial' );
hProcessor = get(hServer, 'Processor' );
hParamStore = get(hTrial, 'ParameterStore' );
hEvStore = get(hTrial, 'EventStore' );

save hTrial hTrial

%trial_description = invoke( hTrial, 'Description');
trial_name = invoke( hTrial, 'ReferenceName');

% Declare the interfaces for the events we will need


dat_path = invoke( hTrial, 'DataPath');
analog_rate = invoke( hTrial, 'AnalogBaseRate');

% Get the EMG data
emglabels_list = {'01' '02' '03' '04' '05' '06' '07' '08' '09'};
nemgChannels = size(emglabels_list,2);

% Determine size of emg data as well as first and last data points
hemg=invoke( hTrial, 'FindAnalogChannel',emglabels_list{3});
emgFirst=invoke( hemg, 'FirstSampleNum')+1;
emgLast=invoke( hemg, 'LastSampleNum')+1;
nemgPoints = emgLast - emgFirst+1;

% Create empty arrays for the emg data
emg_data = zeros(nemgChannels,nemgPoints);
temp_data = zeros(1,nemgPoints);

% Find each channel and save data into emg_data array
for i = 1:nemgChannels
    hemg=invoke( hTrial, 'FindAnalogChannel',emglabels_list{i});

    for j = emgFirst+1:emgLast
        temp_data(j+1)=invoke( hemg, 'GetSample',j);
    end

    % get rid of any offset in the EMG signal
    temp_data=temp_data*10000;
    
    if i == 1
        emg_data = temp_data;
    else
        emg_data = [emg_data;temp_data];
    end

end

%Create emg time array
emg_sample_rate = 1/analog_rate;
emg_Time_end = (length(emg_data)-1) * emg_sample_rate;
emg_time = (0:emg_sample_rate:emg_Time_end);

%now pass a 30Hz highpass butterworth filter to the data
Wn=30/(analog_rate/2);
[B,A]=butter(4,Wn,'high');
for i = 1:nemgChannels
   emg_data(i,:) = filtfilt(B,A,emg_data(i,:));
end

raw_emg = emg_data;

% full wave rectify EMG data
emg_data = abs(emg_data);

%now pass a 6Hz lowpass butterworth filter to the emg
Wn=6/(analog_rate/2);
[B,A]=butter(2,Wn);
for i = 1:nemgChannels
   emg_data(i,:) = filtfilt(B,A,emg_data(i,:));
end

complete_data = [emg_time;emg_data]';
complete_data2 = [emg_time;raw_emg]';

[m,n]=size(complete_data);
% Now output the data to a csv file
new_file = [dat_path trial_name 'EMG_filtered.csv'];
new_file2 = [dat_path trial_name 'EMG_raw.csv'];
fid = fopen(new_file,'w');
fid2 = fopen(new_file2,'w');
% Print header information
fprintf(fid,'Time(sec)')
fprintf(fid2,'Time(sec)')
for i = 1:length(emglabels_list)
    fprintf(fid,',Channel %s',emglabels_list{i});
    fprintf(fid2,',Channel %s',emglabels_list{i});
end
fprintf(fid,'\n');
fprintf(fid2,'\n');

% Print data
for i= 1:m
    fprintf(fid,'%f,',complete_data(i,:));
    fprintf(fid2,'%f,',complete_data2(i,:));
    fprintf(fid,'\n');
    fprintf(fid2,'\n');
end
fclose(fid);
fclose(fid2);

release( hEvStore );
release( hParamStore );
release( hTrial );
release( hServer );


