Trying to extract data trail by trail, but my data is sampled at different sampling rate.

Velocity, position and torque data has sampling frequency of 1000 and EMG's channels at 2000. So I am not getting the right region of intrest. I think I might need two different time vector or is there any other way to solve this problem ?

%% Intialization
fs=960; % sampling frequency of velocity,position and torque
Fs=1960; % sampling frequency of EMG channels 
fc=10;
% Low pass Filter of 10 Hz
[b,a]= butter(1,fc/(fs/2),"low");
prewindow=-1000;
postwindow=5000;
% Caliberation
%torque=12.80409731*Torque;
%tosition=34.2465753*position;
%velocity=
%% Assigning all variables
%Position
position=V300degpersec_Ch2.values;
%position=-position;
positionoffset=min(position);
Position=position+abs(positionoffset);%% remove offset
Position=filtfilt(b,a,position);
% Torque
torque=V300degpersec_Ch4.values;
Torque=filtfilt(b,a,torque);
%Velocity
velocity=V300degpersec_Ch3.values;
velocity=-velocity;
Velocity=filtfilt(b,a,velocity);
Time=1/fs:1/fs:length(Velocity)/fs;
% Raw EMG
RFp=V300degpersec_Ch5.values;% Proximal Rectus Femoris
%RFp=downsample(RFp,2);
RFd=V300degpersec_Ch6.values;% Distal Rectus Femoris 
%RFd=downsample(RFd,2);
VM=V300degpersec_Ch7.values;% Vastus Medialis
%VM=downsample(VM,2);
BiFem=V300degpersec_Ch8.values;% Biceps Femoris 
%BiFem=downsample(BiFem,2);
SM=V300degpersec_Ch9.values;% Semitendinosus 
%SM=downsample(SM,2);
time=1/Fs:1/Fs:(length(SM))/Fs;% Time vector for EMG channels 

%% Croping ROI
moving=(Velocity>1.5);
%plot(Time,Velocity);
%hold on
%plot(Time, moving*150);
x=diff(moving);
figure(1)
plot(x);
start=find(x>0);
finish=find(x<0);
figure(2)
plot(Velocity);
hold on
plot(Time, moving*150);
%tr=2;
trln=max(finish-start);
trln=trln-10;
ntr=length(start);
%Preloacation of variables
% 
% vel=nan(trln,ntr);
% pos=nan(trln,ntr);
% tor=nan(trln,ntr);
% EMG1=nan(trln,ntr);
% EMG2=nan(trln,ntr);
% EMG3=nan(trln,ntr);
% EMG4=nan(trln,ntr);
% EMG5=nan(trln,ntr);
vel={};
pos={};
tor={};
EMG1={};
EMG2={};
EMG3={};
EMG4={};
EMG5={};
for tr=1:ntr
    vel{tr}=Velocity(start(tr):finish(tr));
    pos{tr}=Position(start(tr):finish(tr));
    tor{tr}=Torque(start(tr):finish(tr));
    EMG1{tr}=RFd(start(tr):finish(tr));
    EMG2{tr}=RFp(start(tr):finish(tr));
    EMG3{tr}=VM(start(tr):finish(tr));
    EMG4{tr}=BiFem(start(tr):finish(tr));
    EMG5{tr}=SM(start(tr):finish(tr));
end
for tr=1:ntr-1
    bounds=[start(tr)-prewindow:finish(tr)+postwindow];
    pos{tr}=Position(bounds);
    vel{tr}=Velocity(bounds);
    tor{tr}=Torque(bounds);
    EMG1{tr}=RFd(bounds);
    EMG2{tr}=RFp(bounds);
    EMG3{tr}=VM(bounds);
    EMG4{tr}=BiFem(bounds);
    EMG5{tr}=SM(bounds);
end
%Time vector position,velocity and torque 
for tr=1:ntr-1
    lengthtr=finish(tr)-start(tr);
    t{tr}=-prewindow:1:lengthtr+postwindow;
end

Answer 1

[Jan]

If one signal is measured with 1000 Hz and another with 2000 Hz, it is trivial to determine timepoints of the first one in the time frame of the second one:

time_2000Hz = (time_1000Hz - 1) * 2 + 1

This should allow you to convert the intervals of interest.

[Kiran Isapure]

Yes, but when I am trying to extract data trail by trail from raw data, the EMG channels does not give right region of intrest

moving=(Velocity>1.5);

x=diff(moving);

start=find(x>0);

finish=find(x<0);

trln=max(finish-start);

trln=trln-10;

ntr=length(start);

%Preloacation of variables

% vel=nan(trln,ntr);

% pos=nan(trln,ntr);

% tor=nan(trln,ntr);

% EMG1=nan(trln,ntr);

% EMG2=nan(trln,ntr);

% EMG3=nan(trln,ntr);

% EMG4=nan(trln,ntr);

% EMG5=nan(trln,ntr);

for tr=1:ntr

vel{tr}=Velocity(start(tr):finish(tr));

pos{tr}=Position(start(tr):finish(tr));

tor{tr}=Torque(start(tr):finish(tr));

EMG1{tr}=RFd(start(tr):finish(tr));

EMG2{tr}=RFp(start(tr):finish(tr));

EMG3{tr}=VM(start(tr):finish(tr));

EMG4{tr}=BiFem(start(tr):finish(tr));

EMG5{tr}=SM(start(tr):finish(tr));

end

for tr=1:ntr-1

bounds=[start(tr)-prewindow:finish(tr)+postwindow];

pos{tr}=Position(bounds);

vel{tr}=Velocity(bounds);

tor{tr}=Torque(bounds);

EMG1{tr}=RFd(bounds);

EMG2{tr}=RFp(bounds);

EMG3{tr}=VM(bounds);

EMG4{tr}=BiFem(bounds);

EMG5{tr}=SM(bounds);

end

With the above code I am getting wrong ROI for EMG channel, I was thinking if the diffrent sampling rate has issue.

[Jan]

@Kiran Isapure: I do not see, that you have considered my advice to scale the time for the EMG data anywhere in the code.

Your code looks strange:

for tr=1:ntr
    vel{tr}=Velocity(start(tr):finish(tr));
    pos{tr}=Position(start(tr):finish(tr));
    ...
end
for tr=1:ntr-1
    bounds=[start(tr)-prewindow:finish(tr)+postwindow];
    pos{tr}=Position(bounds);
    vel{tr}=Velocity(bounds);
    ...
end

Why do you overwrite vel and pos for the elements 1:ntr-1 ?

By the way, [a:b] is exactly the same as a:b, but the latter is slightly faster. [ and ] is Matlab's operator for concatenation. [a:b] concates the vector a:b with nothing.

Note that

    bounds=[start(tr)-prewindow:finish(tr)+postwindow];
    pos{tr}=Position(bounds);    

is slower than:

    a = start(tr)-prewindow;
    b = finish(tr)+postwindow;
    pos{tr}=Position(a:b);    

In this code the index vector a:b is not created explicitly, but Matlab copies the block of memory directly. Aftzer creating the index vector at first, Matlab has to check each element, if it is inside the bounds of the array. This take a lot of processing time.

[Kiran Isapure]

Thanks for those inputs, where and how in my code I should scale the time for EMG data.

Answer 2

[Star Strider]

I would downsample the EMG signal using the resample function to do the comparisons.

Retain the original signal of course, since the 2kHz sampling frequency has data you do not want to discard.

I would not upsample the 1kHz data, because that creates interpolated data where none previously existed.

[Kiran Isapure]

Yes, thats the main reson I dont have downsample the EMG singnal, I might discard the EMG signal. I want to index position where the EMG relfex starts. How would I retain and index position where the relflex starts. Thanks

[Star Strider]

I would downsample the EMG data however only for the purpos of comparing it with the instrumentation signal. I would not discard the EMG signal if that is part of the experimental design.

‘How would I retain and index position where the relflex starts.’

I am not certain because I do not have your data and I have no idea what your experimental design is. If there are timestamps on both signals, you could use the timestamps. The best (and in my opinion most robust) option would be to use the timestamps and data to create timetable arrays for each (instrumentation and EMG) signal, and then use synchronize or retime (or both) to work with the data. You would likely need to use interpolation of some sort to define the appropriate times for each signal..

.

[Kiran Isapure]

Thanks for the solution, Can you help me with the syntax, I tried writing, but it give error.

[Star Strider]

I cannot help without your data (and I am not certain I could help if I had it if I do not understand what you want to do and your experimental design).

[Kiran Isapure]

We are testing knee flexor muscles with EMG electrods on rectus femoris and biceps femoris muscle. The participant knee is extened 10 time.

My raw data consist of 10 trails and I am trying to exacrt trail by trail data. But as mentioned before the data is sampled differently for EMG's and velocity position torque.

For post processing code, I am trying to index position where the EMG onset happens and also I am seeing the duration of EMG.

My post processing code is below

%% postprocessEMG

%% Find Region of Intrest(ROI) limits and shift position to start at 0 when ROI starts

%find index where time = 0 (ROI starts)

time_0 = cellfun(@(time_2000Hz) find(time_2000Hz == 0), time_2000Hz, 'UniformOutput',false);

%position when ROI starts

for tr=1:ntr-1

position_0(:,tr) = pos{tr}(1:time_0{tr});

end

%% cell to martix for pos

maxRows = 1;

for xx = 1:length(pos)

if length(pos{xx}) > maxRows

maxRows = length(pos{xx});

end

end

Pos= nan(maxRows,length(pos)); %make NaN array

%loop through your cell array and paste data into NaN array

for yy = 1:length(pos)

Pos(1:length(pos{yy}),yy) = pos{yy};

end

%% shift position to start at 0 when ROI starts

for tr=1:ntr-1

b1(:,tr)= Pos(tr) - position_0(:,tr);

end

%% EMG activation threshold

%create vector of logicals flagging baseline window based on time column

baselineROI=cellfun(@(time_2000Hz) time_2000Hz>=-1000 & time_2000Hz<0, time_2000Hz, 'UniformOutput', false);

%create array with baseline window from EMG stretched column

for tr=1:ntr-1

ebaseline(:,tr) = Rec_filtered_data{tr}(baselineROI{tr});

end

%find mean baseline EMG

for tr=1:ntr-1

ebaseline_mean(:,tr) = mean(ebaseline(tr));

end

%calculate baseline EMG threshold value (mean baseline EMG + 2.5 SD)

STD=std(ebaseline);% STD

for tr=1:ntr-1

ethreshold(:,tr) = (ebaseline_mean(:,tr) + 2 * STD(:,tr));

end

%% Cell to matrix

%loop through to find the largest data set in your cell array

maxRows = 1;

for xx = 1:length(Movav_data)

if length(Movav_data{xx}) > maxRows

maxRows = length(Movav_data{xx});

end

end

Movav= nan(maxRows,length(Movav_data)); %make NaN array

%loop through your cell array and paste data into NaN array

for yy = 1:length(Movav_data)

Movav(1:length(Movav_data{yy}),yy) = Movav_data{yy};

end

%% EMG onset and Offset

for tr=1:10

activation(:,tr)= find(Movav(:,tr) >= ethreshold(:,tr),1,"First");

deactivation(:,tr)= find(Movav(:,tr) >= ethreshold(:,tr),1,"Last");

end

%time index where EMG moving average first increases and decreases over threshold

for tr=1:10

activation_time(:,tr) = time_2000Hz{tr}(activation(:,tr));

end

for tr=1:10

deactivation_time(:,tr)=time_2000Hz{tr}(deactivation(:,tr));

end

% EMG Burst duration

burst_time = deactivation_time - activation_time;

figure(11)

bar(burst_time);

%find EMG activation index

for i=1:10

activation_time_index(i) = find(Rec_filtered_data{i} >= activation_time(:,tr),1,'first');

end

%%position where EMG first increases over threshold

for i=1:10

activation_position(i) = b1(activation_time_index(:,i),1);

end

[Star Strider]

I have no way of determining how well that works (or if it works) without the data it is used with.

[Kiran Isapure]

I have attached position, velocity and EMG data. Hopefully now we can work through this issue. Thanks for looking into this.

[Kiran Isapure]

oh, I forgot to attach the EMG data in pervious chat. Thanks

[Star Strider]

Before I go any farther just now, I note that the the row size of the Position vector is not the same as the row size of the Velocity vector. I assume they should be the same.

How were they calculated?

What should I do with the size information?

I assume that the instrumentation signals were sampled at 1kHz and the EMG at 2kHz, so I’m using those to construct the time vectors for each.

What else do I need to know about these signals?

LD1 = load(websave('Position','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1268180/Position.mat'))

LD1 = struct with fields:
    Pos: [6762×11 double]

Pos = LD1.Pos;

LD2 = load(websave('Velocity','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1268185/Velocity.mat'))

LD2 = struct with fields:
    vel: {[6758×1 double]  [6759×1 double]  [6756×1 double]  [6758×1 double]  [6756×1 double]  [6756×1 double]  [6762×1 double]  [6759×1 double]  [6756×1 double]  [6758×1 double]  [1225×1 double]}

vel = LD2.vel;

LD3 = load(websave('EMG','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1268190/EMG.mat'))

LD3 = struct with fields:
    ALLEMG: {[517468×1 double]}

ALLEMG = LD3.ALLEMG;

return

FsInst = 1E+3;
LInst = 1;
FsEMG = 2E+3;
LEMG = size(EMG,1);
t1 = linspace(0, LInst-1, Linst)/FsInst;                    % Instrumentation Time Vector
t2 = linspace(0, LEMG-1, LEMG)/FsEMG;                    % Instrumentation Time Vector

return

%% postprocessEMG
%% Find Region of Intrest(ROI) limits and shift position to start at 0 when ROI starts
%find index where time = 0 (ROI starts)
time_0 = cellfun(@(time_2000Hz) find(time_2000Hz == 0), time_2000Hz, 'UniformOutput',false);
%position when ROI starts
for tr=1:ntr-1
    position_0(:,tr) = pos{tr}(1:time_0{tr});
end
%% cell to martix for pos
maxRows = 1;
for xx = 1:length(pos)
    if length(pos{xx}) > maxRows
        maxRows = length(pos{xx});
    end
end
Pos= nan(maxRows,length(pos)); %make NaN array
%loop through your cell array and paste data into NaN array
for yy = 1:length(pos)
    Pos(1:length(pos{yy}),yy) = pos{yy};
end
%% shift position to start at 0 when ROI starts
for tr=1:ntr-1
    b1(:,tr)= Pos(tr) - position_0(:,tr);
end
%% EMG activation threshold
%create vector of logicals flagging baseline window based on time column
baselineROI=cellfun(@(time_2000Hz) time_2000Hz>=-1000 & time_2000Hz<0, time_2000Hz, 'UniformOutput', false);
%create array with baseline window from EMG stretched column
for tr=1:ntr-1
    ebaseline(:,tr) = Rec_filtered_data{tr}(baselineROI{tr});
end
%find mean baseline EMG
for tr=1:ntr-1
    ebaseline_mean(:,tr) = mean(ebaseline(tr));
end
%calculate baseline EMG threshold value (mean baseline EMG + 2.5 SD)
STD=std(ebaseline);% STD
for tr=1:ntr-1
    ethreshold(:,tr) = (ebaseline_mean(:,tr) + 2 * STD(:,tr));
end
%% Cell to matrix
%loop through to find the largest data set in your cell array
maxRows = 1;
for xx = 1:length(Movav_data)
    if length(Movav_data{xx}) > maxRows
        maxRows = length(Movav_data{xx});
    end
end
Movav= nan(maxRows,length(Movav_data)); %make NaN array
%loop through your cell array and paste data into NaN array
for yy = 1:length(Movav_data)
    Movav(1:length(Movav_data{yy}),yy) = Movav_data{yy};
end
%% EMG onset and Offset
for tr=1:10
    activation(:,tr)= find(Movav(:,tr) >= ethreshold(:,tr),1,"First");
    deactivation(:,tr)= find(Movav(:,tr) >= ethreshold(:,tr),1,"Last");
end
%time index where EMG moving average first increases and decreases over threshold
for tr=1:10
    activation_time(:,tr) = time_2000Hz{tr}(activation(:,tr));
end
for tr=1:10
    deactivation_time(:,tr)=time_2000Hz{tr}(deactivation(:,tr));
end
% EMG Burst duration
burst_time = deactivation_time - activation_time;
figure(11)
bar(burst_time);
%find EMG activation index
for i=1:10
    activation_time_index(i) = find(Rec_filtered_data{i} >= activation_time(:,tr),1,'first');
end
%%position where EMG first increases over threshold
for i=1:10
    activation_position(i) = b1(activation_time_index(:,i),1);
end

[Kiran Isapure]

Thanks helping out with time vector. The position is matrix and velcotiy is cell.

%find EMG activation index
for i=1:10
    activation_time_index(i) = find(Rec_filtered_data{i} >= activation_time(:,tr),1,'first');
end
%%position where EMG first increases over threshold
for i=1:10
    activation_position(i) = b1(activation_time_index(:,i),1);
end

this part of code is not working

[Jan]

@Kiran Isapure: Please use the tools above the section for typing comments to format code as code. This improves the readability.

Explain "not working" with any details. It is easier to fix a problem than to guess, what the problem is.

[Star Strider]

I atill cannot run the posted code.

See the error message —

LD1 = load(websave('Position','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1268180/Position.mat'))

LD1 = struct with fields:
    Pos: [6762×11 double]

Pos = LD1.Pos;

LD2 = load(websave('Velocity','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1268185/Velocity.mat'))

LD2 = struct with fields:
    vel: {[6758×1 double]  [6759×1 double]  [6756×1 double]  [6758×1 double]  [6756×1 double]  [6756×1 double]  [6762×1 double]  [6759×1 double]  [6756×1 double]  [6758×1 double]  [1225×1 double]}

vel = LD2.vel;

LD3 = load(websave('EMG','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1268190/EMG.mat'))

LD3 = struct with fields:
    ALLEMG: {[517468×1 double]}

ALLEMG = LD3.ALLEMG;

% return

FsInst = 1E+3;
LInst = 1;
FsEMG = 2E+3;
LEMG = size(ALLEMG,1);
t1 = linspace(0, LInst-1, LInst)/FsInst;                    % Instrumentation Time Vector
t2 = linspace(0, LEMG-1, LEMG)/FsEMG;                    % Instrumentation Time Vector

% return

%% postprocessEMG
%% Find Region of Intrest(ROI) limits and shift position to start at 0 when ROI starts
%find index where time = 0 (ROI starts)
time_0 = cellfun(@(time_2000Hz) find(time_2000Hz == 0), time_2000Hz, 'UniformOutput',false);

Unrecognized function or variable 'time_2000Hz'.

%position when ROI starts
for tr=1:ntr-1
    position_0(:,tr) = pos{tr}(1:time_0{tr});
end
%% cell to martix for pos
maxRows = 1;
for xx = 1:length(pos)
    if length(pos{xx}) > maxRows
        maxRows = length(pos{xx});
    end
end
Pos= nan(maxRows,length(pos)); %make NaN array
%loop through your cell array and paste data into NaN array
for yy = 1:length(pos)
    Pos(1:length(pos{yy}),yy) = pos{yy};
end
%% shift position to start at 0 when ROI starts
for tr=1:ntr-1
    b1(:,tr)= Pos(tr) - position_0(:,tr);
end
%% EMG activation threshold
%create vector of logicals flagging baseline window based on time column
baselineROI=cellfun(@(time_2000Hz) time_2000Hz>=-1000 & time_2000Hz<0, time_2000Hz, 'UniformOutput', false);
%create array with baseline window from EMG stretched column
for tr=1:ntr-1
    ebaseline(:,tr) = Rec_filtered_data{tr}(baselineROI{tr});
end
%find mean baseline EMG
for tr=1:ntr-1
    ebaseline_mean(:,tr) = mean(ebaseline(tr));
end
%calculate baseline EMG threshold value (mean baseline EMG + 2.5 SD)
STD=std(ebaseline);% STD
for tr=1:ntr-1
    ethreshold(:,tr) = (ebaseline_mean(:,tr) + 2 * STD(:,tr));
end
%% Cell to matrix
%loop through to find the largest data set in your cell array
maxRows = 1;
for xx = 1:length(Movav_data)
    if length(Movav_data{xx}) > maxRows
        maxRows = length(Movav_data{xx});
    end
end
Movav= nan(maxRows,length(Movav_data)); %make NaN array
%loop through your cell array and paste data into NaN array
for yy = 1:length(Movav_data)
    Movav(1:length(Movav_data{yy}),yy) = Movav_data{yy};
end
%% EMG onset and Offset
for tr=1:10
    activation(:,tr)= find(Movav(:,tr) >= ethreshold(:,tr),1,"First");
    deactivation(:,tr)= find(Movav(:,tr) >= ethreshold(:,tr),1,"Last");
end
%time index where EMG moving average first increases and decreases over threshold
for tr=1:10
    activation_time(:,tr) = time_2000Hz{tr}(activation(:,tr));
end
for tr=1:10
    deactivation_time(:,tr)=time_2000Hz{tr}(deactivation(:,tr));
end
% EMG Burst duration
burst_time = deactivation_time - activation_time;
figure(11)
bar(burst_time);
%find EMG activation index
for i=1:10
    activation_time_index(i) = find(Rec_filtered_data{i} >= activation_time(:,tr),1,'first')
end
%%position where EMG first increases over threshold
for i=1:10
    activation_position(i) = b1(activation_time_index(:,i),1)
end

.

[Kiran Isapure]

Sorry for that, I am trying to index activiation_time_index from b1( position ROI). Hope this will help

[Star Strider]

The error:

Unrecognized function or variable 'time_2000Hz'.

is still unresolved, so I can’t run the code!