Apply a windowed short-time discrete Fourier transform to filter out everything but the fundamental frequency?

Hi All,

I am trying to calculate the relative Fourier phase of the trunk/pelvis rotation using the methods described here: https://www.ncbi.nlm.nih.gov/pubmed/29543108

I have identified the fundamental frequency of my signal (4000 samples) and now I would like to determine the 'windowed' Fourier phase angles for it.

They identify the steps as:

Select the window length as 4x the period of the fundamental frequency (1.45Hz, fs: 200Hz) is:

window_length = 4*(1/1.45)*200

Apply a short-time discrete Fourier transfomr to each window to determine the fourier phase angles of the fundamental frequency at each sample (from sample 1 to total number of samples - window length +1).

Run the SDFT for all windows; a new window (j) is formed by shifting data by one sample. Obtain the phase angles from each window s(j).

Transform s(j) into a continuous time series.

I am struggling with the last two steps and could really use a point in the right direction. My script below:

%FFT to calculate the estimated PSD for the segmental angles 

F=ut;                   % Segment (4000 samples at 200Hz)
F(isnan(F)) = []        % Eliminate NaN
Fs = 200;               % Sampling Freq
T = 1/Fs;               % Sampling Period
L = length(F)           % Length of signal 
t = (0:L-1)*T;          % Time vector 
FF = fft(F)             % Fourier Series of Data, Freq Vector
P2 = abs(FF/L)          % two sided spectrum P2
P1 = P2(1:L/2+1);       % One sided spectrum P1
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;     % Frequency Domain

%Find fundamental frequency and first 5 harmonics
[pks, locs]=findpeaks(P1, 'minpeakdistance', 20,  'minpeakheight',0.001);
hold on;
harm = sort(pks, 'descend');
harm = harm(1:6);
for k = 1:6
idx(k) = find(pks == harm(k));
end
harm = horzcat(locs(idx(:,1)), pks(idx(:,1)));
for k = 1:6
    harm(k,3) = f(harm(k,1));
end
plot(harm(:,3), harm(:,2), 'or')

%Bin PSD by +/- 20 samples
for k = 1:6
    if harm(k,1)-20 < 0
        harm(k,4) = sum(P1(1:harm(k,1)+20));
    else
          harm(k,4) = sum(P1(harm(k,1)-20:harm(k,1)+20));
end
end

fund_freq = harm(1,3);

%window length = 4(1/f1max)samplingrate
L = (4*(1/fund_freq)*200);
win = rectwin(L);

%short time discrete fourier transform 
s = stft(F,d, 'window', win, 'OverlapLength', 550, 'FFTlength', 3448);
s = angle(s);
s = rad2deg(s);

I am not sure if I am doing the STFT right.

I am also not sure ohow to transform s into a continuous time series after performing the STFT.

I have attached the 'ut' variable as an example signal.

Your help is very much appreciated.