Ideal way to smoothe gyroscope data
Show older comments
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.
Accepted Answer
Star Strider
on 24 Jan 2023
0 votes
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.
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.
.
8 Comments
Ines Shekhovtsov
on 24 Jan 2023
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
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.
), 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
on 24 Jan 2023
Star Strider
on 24 Jan 2023
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
on 24 Jan 2023
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)
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))
NrNaNCol = [nnz(isnan(Gs(:,1))) nnz(isnan(Gs(:,2))) nnz(isnan(Gs(:,3)))]
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
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
on 25 Jan 2023
Star Strider
on 25 Jan 2023
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.
.
More Answers (0)
Categories
Find more on Get Started with Signal Processing Toolbox in Help Center and File Exchange
on 24 Jan 2023
on 25 Jan 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
