Ideal way to smoothe gyroscope data

Hello,

I have a .csv file that contains EMG and gyroscope data collected from the Delsys Trigno Avanti System. These wireless electrodes collect EMG data at 1924.233Hz and gyroscope data at 74.0741Hz (Honestly i'm a bit confused as to why the sampling frequency is not a whole number). I know how to work with EMG data but this is my first time working with gyroscope data. Each sensor collects gyroscope data in x,y and z directions. I am creating a figure that shows EMG data on left and the corresponding gyroscope data in all three directio on the right. The x, y and z direction plots are overlayed on top of each other. I would like to smooth this gyroscope data and a brief google search made me even more unsure of what the best approach might be for this. I know it depends on what we want to do with the data but i am not sure of that yet. As of now i just want to smoothe the data. Previously for smoothing we have used a LPF at 10Hz. I am just curious what you all would suggest. Thank you for your time.

Answer 1

[Star Strider]

The LPF approach is likely appropriate, and an IIR filter (elliptical is most efficient) is likely the best option.

First, determine if the sampling intervals are constant. This can be done by taking the mean and std of the differences (diff) of the time vector. If the standard deviation is small (on the order of or so times the mean value), the time vector can be used without resampling. If it is larger than that, it might first be necessary to use the resample function to resample it to a reliable time vector. Use a sampling frequency at least 2 times ‘1/mean(diff(t))’ for this. Most signal processing (and all filtering) require a regularly-sampled time vector to work correctly with the signal.

The problem then may be in selecting the passband and stopband frequencies. The best way to do that is to calculate the fft of the signal, and then search the plot of it for the ‘best’ passband frequency. Sometimes, this can be automated, however most times it will require experimenting with different frequencies to determine which one produces the best result when it is used to filter the signal.

The fft will also help to determine if there’s broadband noise in the signal. If there is, the two options to reduce it are to use wavelet denoising on the original signal, or use the Savitzky-Golay filter (sgolayfilt) on the original signal. Again, the sgolayfilt arguments will require experimentation. I usually use a 3-degree polynomial and then choose the ‘framelen’ parameter to produce the best result. Do that before doing frequency-selective lowpass filtering.

This outlines the workflow I usually use in processing experimental data. I can help you if you need my help.

.

[Ines Shekhovtsov]

Thank you for your response. The file in question is 290.493 seconds long. The gyroscope data has 21518 samples in the file total. I used the following code to determine if the sampling intervals are constant:

time_vector = linspace(0, 290.493, 21518);

diff_vector = diff(time_vector);

unique_diff = unique(diff_vector);

if length(unique_diff) == 1

disp('The time vector has constant sampling intervals.')

else

disp('The time vector does not have constant sampling intervals.')

end

And they are NOT. The sampling fequency of gyro data is approx 74Hz as mentioned above. So i would resample to at least 2x that to 150Hz. I can do that. My next question would be how to calculate the fft. I myself have done it different ways and with different code. Could you please let me know how you would do it? And this is all a bit new to me but would i have to look at the fft of all of the gyro data individually or just pick one direction (x, y or z). I have 8 emg sensors and each of those will give me gyro data in x, y, z directions for a total of 24 columns of gyro data. I am assuming i need to plot fft of all of these??

Thank you for your help. I need to get better with frequency analysi in general and i appreciate your help.

[Star Strider]

It would help to have the gyroscope data to work with.

I usually use this approach to calculate and plot the fft and this needs to be done after doing the resample step to display a one-sided Fourier Transform —

Fs = 150;                                           % Sampling Frequency
L = 20000;                                          % Length Of Signal
t = linspace(0, L-1, L).'/Fs;                       % Time Vector (After Using 'resample')
Gs = ([0;1;2]+sin([1;3;5]*2*pi*t.')).';             % Gyroscope Signal

figure
plot(t,Gs)
xlabel('Time')
ylabel('Signal')
legend('x', 'y', 'z', 'Location','best')
xlim([0 5])                                         % Show Detail

Fn = Fs/2;                                              % Nyquist Frequency
L = numel(t);                                           % Signal Length
NFFT = 2^nextpow2(L)                                    % For Efficiency, Also Increases Frequency Resolution

NFFT = 32768

FT_Gs = fft((Gs-mean(Gs)).*hann(L), NFFT)/L;            % Fourier Transform (Windowing Is Optional)
Fv = linspace(0, 1, NFFT/2+1)*Fn;                       % Frequency Vector (One Way To Calculate It)
Iv = 1:numel(Fv);                                       % Index Vector

figure
plot(Fv, abs(FT_Gs(Iv,:))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
legend('x', 'y', 'z', 'Location','best')

figure
plot(Fv, abs(FT_Gs(Iv,:))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
legend('x', 'y', 'z', 'Location','best')
xlim([0 6])                                                 % 'Zoom' To Dhow Detail

The ‘FT_Gs’ assignment calculation needs a bit of explanation.

The Fourier transform plots a D-C (constant) offset at 0 Hz if the signal has an offset. Because of appropriate scaling in the plot function, if a D-C offset is present that is significantly higher than the rest of the signal peaks, this can hide them. For that reason, I subtract the D-C offset (here ‘mean(Gs)’) before calculating the Fourier transform.

Using the window (here the Hanning window hann) compensates for the fact that the fft is not infinite in both directions (to ), and produces an acceptable result as though the Fourier transform was infinite.

Using a power-of-two (‘NFFT’) to zero-pad the vector before calculating the fft significantly improves the efficiency of the fft calculation, and also increases the frequency resolution of the result.

Dividing the result of the fft calculation by the length (‘L’) of the original signal normalises it so the amplitudes of the fft are close to the amplitudes of the original time-domain signals.

I hope that you can use this code essentially ‘as is’ with your signals.

.

[Ines Shekhovtsov]

Thank you! I will look over this code and try to incorporate it with my data. The file I am using is 112MB. I can delete a decent chunk of data but it will still be just over 20MB. Is there a way that i can upload that through this forum? I believe i used websave function before but that documentation looks a bit confusing to me now.

[Star Strider]

As always, my pleasure!

The webwrite function might be more appropriate, however you can preferably upload it using the ‘paperclip’.icon (just to the right of the Σ icon in the top toolstrip) in a Comment or when editing your original post.

A 1 MB or less chunk of the file will likely be enough, considering the relatively low sampling rate, even when resampled to 150 Hz. That should be enough to design an appropriate filter and any other necessary preprocessing for it.

[Ines Shekhovtsov]

Thank you! Was able to find a much smaller file and delete further data from it. It has unnecessary information for you but i wanted you to see the .csv file structure when it is exported from the acquistion software. You will see EMG (for 3 muscles), gyro (for those same 3 muscles) and some analog sensor data.

[Star Strider]

I remember your saying that some of the signals need to be resampled.

Do any of these have associated time vectors, or have they been resampled to the correct sampling frequencies? I’m assuming they’ve already been resampled. If they still need to be, I need time vectors for them. (I’m considering only the gyroscope signals at this point.)

% file = websave('Trial_5','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1273265/Trial_5.csv');
% type(file)

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1273265/Trial_5.csv', 'VariableNamingRule','preserve', 'HeaderLines',5)

T1 = 31500×20 table
    EMG 1 (1926Hz)    GYRO X (74Hz)    GYRO Y (74Hz)    GYRO Z (74Hz)    EMG 1 (1926Hz)_1    GYRO X (74Hz)_1    GYRO Y (74Hz)_1    GYRO Z (74Hz)_1    EMG 1 (1926Hz)_2    GYRO X (74Hz)_2    GYRO Y (74Hz)_2    GYRO Z (74Hz)_2    EMG 1 (2222Hz)    EMG 2 (2222Hz)    EMG 3 (2222Hz)    EMG 4 (2222Hz)    Ana5 (2222Hz)    Ana6 (2222Hz)    Ana7 (2222Hz)    Ana8 (2222Hz)
    ______________    _____________    _____________    _____________    ________________    _______________    _______________    _______________    ________________    _______________    _______________    _______________    ______________    ______________    ______________    ______________    _____________    _____________    _____________    _____________

      -0.0038605         -3.9237          -2.4885          -1.1603           -0.011246           3.6031             4.4427              1.4962           -0.012589            -7.7863            4.5954             -2.0916            5.1748            5.1966            5.1679            5.1369           5.2034           5.1952           1.4115           0.26578   
      -0.0041962         -3.9847          -2.2137          -1.1756           -0.013931           3.6031             4.8092              1.4962           -0.012253             -7.771            4.8855             -1.8779            5.1743            5.1969             5.169             5.138           5.2032           5.1951           1.4111           0.26578   
      -0.0028534         -4.1374          -2.1069          -1.1756           -0.011749           3.7557             4.9008              1.6336            -0.01091            -7.8473            4.9924             -1.8626            5.1753            5.1968            5.1689            5.1375           5.2038           5.1958           1.4112           0.26546   
       -0.002182         -4.0305          -1.9847          -1.2214           -0.011246            3.542             4.6718              1.5573           -0.012924            -7.8779            4.8092             -1.8015            5.1751            5.1969            5.1684            5.1375            5.204           5.1955           1.4112           0.26546   
      -0.0031891              -4          -2.1527          -1.2519           -0.010407           3.4504             4.4275              1.5267           -0.012421            -7.8473             4.626             -1.9389            5.1753            5.1973            5.1687            5.1375            5.204           5.1955           1.4108           0.26546   
      -0.0046998         -3.9847          -2.4122          -1.1603           -0.013596           3.2977             4.2443              1.4198           -0.012253            -7.8473            4.5038                  -2            5.1748            5.1965            5.1678             5.137           5.2037           5.1951           1.4111           0.26562   
      -0.0025177         -4.0611          -2.5191          -1.0992           -0.010239           3.1756             4.2137              1.2824           -0.014603            -7.6947            4.3053             -2.0763            5.1747            5.1966            5.1684            5.1378           5.2032           5.1954           1.4112           0.26594   
      -0.0036927         -4.3206           -2.626         -0.96183           -0.013596           3.2824             4.3817              1.2519           -0.013764            -7.6031            4.3817             -2.1221             5.175            5.1971            5.1689            5.1378           5.2032           5.1955           1.4117           0.26594   
       -0.010575         -4.2595          -2.8244         -0.71756          -0.0095674           3.0229             4.0611               1.084           -0.012924            -7.4198            4.3206             -2.3511            5.1753            5.1969            5.1685             5.137            5.204           5.1958           1.4117           0.26625   
       -0.013931         -4.2137          -3.2824         -0.80916          -0.0099031           2.9466             3.6336              1.0534           -0.010239            -7.4046            3.6794             -2.5649            5.1747            5.1968            5.1682             5.138           5.2032           5.1955           1.4115           0.26625   
       -0.019303         -4.2748          -3.8015         -0.71756          -0.0087282           2.7481             3.4656             0.96183           -0.012924            -7.3588            3.2672             -2.7023            5.1759            5.1977            5.1689            5.1372           5.2051           5.1968           1.4117           0.26531   
       -0.022995         -4.3359          -3.8931         -0.62595            -0.01326           2.6107             3.4809             0.77863           -0.011582            -7.3588            3.1908             -2.8397            5.1747            5.1963            5.1682            5.1373           5.2032           5.1951           1.4114           0.26562   
       -0.021317         -4.4122          -3.8473         -0.59542           -0.011582           2.3969             3.5725             0.82443           -0.010575            -7.2366            3.1298             -2.9313            5.1745            5.1962            5.1682            5.1375           5.2029           5.1946           1.4111           0.26562   
       -0.019303         -4.4733          -3.9389         -0.59542           -0.015442           2.4275             3.4351             0.51908           -0.010239            -7.1756            3.0229             -2.9313            5.1747            5.1973            5.1685            5.1378           5.2037           5.1952           1.4114           0.26546   
       -0.019974         -4.5802          -4.1679         -0.62595           -0.017624           2.3511             3.3435             0.33588           -0.011246             -7.145            2.7328             -2.9008            5.1751            5.1969            5.1685            5.1373           5.2037           5.1954            1.412           0.26609   
        -0.01796         -4.4427          -4.0305         -0.68702           -0.010407           2.3053             3.1756             0.41221           -0.013764            -7.1603            2.8092             -2.7939             5.175            5.1966            5.1684            5.1373           5.2038           5.1951           1.4114           0.26594   

VN = T1.Properties.VariableNames;

% Fs = 150;                                               % Sampling Frequency
Fs = 74;
% L = 20000;                                              % Length Of Signal
L = size(T1,1);
t = linspace(0, L-1, L).'/Fs;                           % Time Vector (After Using 'resample')
% Gs = ([0;1;2]+sin([1;3;5]*2*pi*t.')).';                 % Gyroscope Signal
Gs = T1{:,[2 3 4]};
NrNaN = nnz(isnan(Gs))

NrNaN = 91350

NrNaNCol = [nnz(isnan(Gs(:,1))) nnz(isnan(Gs(:,2))) nnz(isnan(Gs(:,3)))]

NrNaNCol = 1×3

30450 30450 30450

Gs = fillmissing(Gs, 'nearest');                        % Interpolate Missing Values

figure
plot(t,Gs)
grid
xlabel('Time')
ylabel('Signal')
legend(VN{2:4}, 'Location','best')
xlim([0 14])                                             % Show Detail

Fn = Fs/2;                                              % Nyquist Frequency
L = numel(t);                                           % Signal Length
NFFT = 2^nextpow2(L)                                    % For Efficiency, Also Increases Frequency Resolution

NFFT = 32768

FT_Gs = fft((Gs-mean(Gs)).*hann(L), NFFT)/L;            % Fourier Transform (Windowing Is Optional)
% FT_Gs = fft((Gs-mean(Gs)).*kaiser(L), NFFT)/L;          % Fourier Transform (Windowing Is Optional)
Fv = linspace(0, 1, NFFT/2+1)*Fn;                       % Frequency Vector (One Way To Calculate It)
Iv = 1:numel(Fv);                                       % Index Vector

figure
semilogy(Fv, abs(FT_Gs(Iv,:))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
legend(VN{2:4}, 'Location','best')
xlim([0 Fn])

figure
semilogy(Fv, abs(FT_Gs(Iv,:))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
legend(VN{2:4}, 'Location','best')
xlim([0 10])                                                % 'Zoom' To Dhow Detail

NrSp = size(FT_Gs,2);

figure
for k = 1:NrSp
    subplot(NrSp,1,k)
    semilogy(Fv, abs(FT_Gs(Iv,k))*2)
    grid
    ylabel('Amplitude')
    title(VN{k+1})
    xlim([0 10])
end
xlabel('Frequency (Hz)')

% ylabel('Amplitude')
% legend('x', 'y', 'z', 'Location','best')

Gs_filt = lowpass(Gs, 5, Fs, 'ImpulseResponse','iir');      % Create 'ellip' Filter, Filter All Signals Simultaneously

figure
for k = 1:NrSp
    subplot(NrSp,1,k)
    plot(t, Gs_filt(:,k))
    grid
    ylabel('Amplitude')
    title(sprintf('Filtered %s',VN{k+1}))
    xlim([0 14])                                             % Show Detail
end
xlabel('Time (s)')

There were apparently some (30450 in each column) NaN values in the original signal that I discovered when there was no output from the fft call that could be plotted (blank plots). That’s in the original data. The fillmissing call did a 'nearest' intepolation to eliminate them, and then the fft plot appeared.

There doesn’t appear to be any actual signal after about 14 seconds. Up until then however, the signals are relatively well-behaved and the filter produces reasonable results. Setting the lowpass filter passband to 5 Hz produces a reasonably detailed yet relatively noise-free signal. Lower it further to remove more high-frequency noise.

.

[Ines Shekhovtsov]

Thank you for all of your help! The raw data in the .csv will need to be resampled. I believe for a lot of my future analysis i will resample EMG and sync data to 2000Hz. The gyro data will only be reseampled to 150Hz so as to not lose too much of original. If you can see from the .csv file the output doesn't manually create a time vector column since there are three different types of signals (EMG, gyroscope and analog signal) all sampled at different frequencies. In the beginning of my code, i removed the NaN values as they are just there to match data point to data point up to the highest sampled signal (The analog is sampled at 2222.2222Hz). And yes, i gave you a small file that is pretty short which is why there is no more signal after 14s. As far as time vector for this analysis, i am assuming i can manually create that based on the number of seconds in the file (14.175s), the number of data points for each signal type and their acquisition frequencies.

[Star Strider]

As always, my pleasure!

It will be necessary to have the original time vectors for the gyroscope signals in order to resample them correctly to a new, constant-interval, time vector.

I use the syntax:

[sr,tr] = resample(s,t,Fs);

for such problems, where ‘s’ is the original signal,‘t’ the original time vector, ‘Fs’ the desired sampling frequency of the resampled signal (it can be anything that reasonably preserves the original signal data without doing excessive interpolation, here likely 70 Hz to 80 Hz, rethinking my earlier suggestion of 150 Hz), ‘sr’ is the resampled signal, and ‘tr’ is the resampled time vector.

For regularly-sampled signals with a known sampling frequency, you can calculate the time vector the same way I calculated ‘t’ in my code.

The rest of my code should work essentially as written for the resampled signals, with necessary changes to some of the parameters. I will of course help to get it to run with the resampled data if there are any problems with it.

I did not include any of the EMG data in my code, since I understood that the gyroscope data were the primary concern. I can include the EMG data if I understand what you want to do with them.

.