Converting tremor data in to frequency and filter data

Can anyone help in apply fft and filtering to my IMU sensor data for detecting tremor

[Jeffrey Clark]

@Paras Shaikh, the csv you provided doesn't provide time readings or units for the measurements, and should we just assume from the last column that there are three separate datasets? If the samples are at the exact same spacing of time please provide that and the acceleration and position units. Are the three sets related and to be evaluated over the same time period or just stuck in one file for convieniance?

[Star Strider]

I need a corresponding time vector (preferably) and if the data are regularly-sampled, the sampling frequency.

Also, what is the ‘label’ variable, and is it important?

[Paras Shaikh]

I took each sample per 1sec.. so time is exact same as samples.. 4500 samples.. So samplings frequency must be 4500

[Star Strider]

‘I took each sample per 1sec.. ’

Then the sampling frequency is 1 sample/sec or 1 Hz.

I have changed my code in my Answer and subsequent Comment to reflect that.

Answer 1

[Star Strider]

The fft plots are straightforward, however only the acceleration signals make sense. The ‘Gyr’ signals exhibit broadband noise.

How do you want to process them?

One approach —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1124040/labeled%20final%20data.csv')

T1 = 4502×7 table
    xAcc     yAcc     zAcc    xGyro    yGyro    zGyro    label
    _____    _____    ____    _____    _____    _____    _____

     0.02    -0.06    0.99     0.12    -0.67    -0.24      1  
     0.02    -0.08    0.98     2.93    -2.56    -7.14      1  
     0.02    -0.09    0.98     7.32    -0.79     1.59      1  
    -0.01    -0.09    0.99    -3.97    -0.31     0.49      1  
    -0.01    -0.08    0.99    -3.54    -2.08    -1.71      1  
    -0.02    -0.09    0.99    -1.22     0.55    -1.95      1  
    -0.04    -0.08    0.99    -0.18    -0.24     0.12      1  
    -0.03    -0.09    0.99     0.49    -0.06     0.24      1  
    -0.04    -0.09    0.99     0.12    -0.43     0.12      1  
    -0.04    -0.09    0.99     0.24    -0.12    -0.12      1  
    -0.04    -0.09    0.99     0.24     0.37     0.12      1  
    -0.04    -0.09    0.99     0.12     0.18    -0.37      1  
    -0.04    -0.09    0.99     0.06     0.31     0.24      1  
    -0.04    -0.09    0.99    -0.31     0.06     0.18      1  
    -0.04    -0.08    0.99    -0.79    -0.24     0.06      1  
    -0.05    -0.09    0.99     0.12    -0.06    -0.18      1  

Acc = T1{:,1:3}

Acc = 4502×3

0.0200 -0.0600 0.9900 0.0200 -0.0800 0.9800 0.0200 -0.0900 0.9800 -0.0100 -0.0900 0.9900 -0.0100 -0.0800 0.9900 -0.0200 -0.0900 0.9900 -0.0400 -0.0800 0.9900 -0.0300 -0.0900 0.9900 -0.0400 -0.0900 0.9900 -0.0400 -0.0900 0.9900

Gyr = T1{:,4:6};

Data = [Acc Gyr];

Fs = 1;                                                     % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
NrSp = size(Data,2);                                        % Number Of Subplots
L = size(Data,1);                                           % Length Of Data Vectors
t = linspace(0, L-1, L)/Fs;                                 % Time Vector

NFFT = 2^nextpow2(L);                                       % For Efficiency
FTData = fft(Data - mean(Data),NFFT)/L;                     % Fopurier Transform
Fv = linspace(0, 1, NFFT/2-1)*Fn;                           % Frequency Vector (For Plots)
Iv = 1:numel(Fv);                                           % Index Vector (For Plots)

figure
sp = [1:2:NrSp 2:2:NrSp];                                   % Order 'subplot' Plots
for k = 1:NrSp
    subplot(3,2,sp(k))
    plot(Fv, abs(FTData(Iv,k))*2)                           % Plot Fourier Transforms
    grid
    xlabel('Frequency')
    ylabel('Magnitude')
    xlim([0 Fn/20])
%     xlim([5 10])
    title(sprintf('Column %d',k))
end
sgtitle('Fourier Transform of Data')

figure
sp = [1:2:NrSp 2:2:NrSp];
for k = 1:NrSp
    subplot(3,2,sp(k))
    plot(t, Data(:,k))
    grid
    ylim([-3 3])
    xlabel('Time')
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
sgtitle('Original Time-Domain Data')

Wp = [0.0001 0.005]/Fn;                                      % Define Passband In Hz, Normalise To (0,pi)
Ws = Wp.*[0.95 1.05];                                       % Stopband (Normalised)
Rs = 50;                                                    % Stopband Ripple (Attenuation) dB
Rp =  1;                                                    % Passband Ripple dB
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp);                                % Design Filter
[sos,g] = zp2sos(z,p,k);                                    % Implement Filter As Second-Order-Section Representation

figure
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 0.01])
set(subplot(2,1,2), 'XLim',[0 0.01])

Data_Filt = filtfilt(sos,g,Data);

figure
for k = 1:NrSp
    subplot(3,2,sp(k))
    plot(t, Data_Filt(:,k))
    grid
    ylim([-1 1]*0.75)
    xlabel('Time')
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
sgtitle('Time-Domain Plots of Bandpass (0.0001 - 0.005 Hz) Filtered Data')

EDIT — (14 Sep 2022 at 11:40)

Changed ‘Fs’ and filter passbands to conform to the 1 Hz sampling frequency.

.

[Paras Shaikh]

Thankyu for your efforts.. If i want one graph for accelerometer data of first 3 columns.. what changes should i apply?

[Star Strider]

My pleasure!

They will overplot each other if plotted on the same axes, so I would use subplot plots in one figure.

Try something like this —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1124040/labeled%20final%20data.csv')

T1 = 4502×7 table
    xAcc     yAcc     zAcc    xGyro    yGyro    zGyro    label
    _____    _____    ____    _____    _____    _____    _____

     0.02    -0.06    0.99     0.12    -0.67    -0.24      1  
     0.02    -0.08    0.98     2.93    -2.56    -7.14      1  
     0.02    -0.09    0.98     7.32    -0.79     1.59      1  
    -0.01    -0.09    0.99    -3.97    -0.31     0.49      1  
    -0.01    -0.08    0.99    -3.54    -2.08    -1.71      1  
    -0.02    -0.09    0.99    -1.22     0.55    -1.95      1  
    -0.04    -0.08    0.99    -0.18    -0.24     0.12      1  
    -0.03    -0.09    0.99     0.49    -0.06     0.24      1  
    -0.04    -0.09    0.99     0.12    -0.43     0.12      1  
    -0.04    -0.09    0.99     0.24    -0.12    -0.12      1  
    -0.04    -0.09    0.99     0.24     0.37     0.12      1  
    -0.04    -0.09    0.99     0.12     0.18    -0.37      1  
    -0.04    -0.09    0.99     0.06     0.31     0.24      1  
    -0.04    -0.09    0.99    -0.31     0.06     0.18      1  
    -0.04    -0.08    0.99    -0.79    -0.24     0.06      1  
    -0.05    -0.09    0.99     0.12    -0.06    -0.18      1  

Acc = T1{:,1:3};
Gyr = T1{:,4:6};

Data = Acc;

Fs = 1;                                                     % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
NrSp = size(Data,2);                                        % Number Of Subplots
L = size(Data,1);                                           % Length Of Data Vectors
t = linspace(0, L-1, L)/Fs;                                 % Time Vector

NFFT = 2^nextpow2(L);                                       % For Efficiency
FTData = fft(Data - mean(Data),NFFT)/L;                     % Fopurier Transform
Fv = linspace(0, 1, NFFT/2-1)*Fn;                           % Frequency Vector (For Plots)
Iv = 1:numel(Fv);                                           % Index Vector (For Plots)

figure
% sp = [1:2:NrSp 2:2:NrSp];                                   % Order 'subplot' Plots
for k = 1:3
    subplot(3,1,k)
    plot(Fv, abs(FTData(Iv,k))*2)                           % Plot Fourier Transforms
    grid
    xlabel('Frequency')
    ylabel('Magnitude')
    xlim([0 Fn/20])
%     xlim([5 10])
    ylim([0 0.4])
    title(sprintf('Column %d',k))
end
sgtitle('Fourier Transform of Data')

figure
% sp = [1:2:NrSp 2:2:NrSp];
for k = 1:NrSp
    subplot(3,1,k)
    plot(t, Data(:,k))
    grid
    ylim([-3 3])
    xlabel('Time')
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
sgtitle('Original Time-Domain Data')

Wp = [0.0001 0.0075]/Fn;                                    % Define Passband In Hz, Normalise To (0,pi)
Ws = Wp.*[0.95 1.05];                                       % Stopband (Normalised)
Rs = 50;                                                    % Stopband Ripple (Attenuation) dB
Rp =  1;                                                    % Passband Ripple dB
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp);                                % Design Filter
[sos,g] = zp2sos(z,p,k);                                    % Implement Filter As Second-Order-Section Representation

figure
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 0.01])
set(subplot(2,1,2), 'XLim',[0 0.01])
sgtitle('Filter Bode Plot')

Data_Filt = filtfilt(sos,g,Data);

figure
for k = 1:3
    subplot(3,1,k)
    plot(t, Data_Filt(:,k))
    grid
    ylim([-1 1]*0.75)
    xlabel('Time')
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
sgtitle('Time-Domain Plots of Bandpass (1E-4 - 75E-4 Hz) Filtered Data')

Adjust the upper passband (the second value in ‘Wp’, 0.0075 Hz in this example) to get the filtered result you want. Use the fft result to guide that choice. Use the same filter for all the signals.

.

[Paras Shaikh]

Hey my required freq is in between 2 to 8hz for this data. Am i getting it? If im changing my data into freq domain??

[Star Strider]

With a sampling frequency of 1 Hz, as you specified in your earlier Comment, the highest frequency you can even see in this signal is 0.5 Hz (the Nyquist frequency).

However if you meant that each column was 4500 samples in one second, then that makes the sampling frequency is 4500 Hz.

There is not much information in the filtered signal —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1124040/labeled%20final%20data.csv');
Acc = T1{:,1:3};
Gyr = T1{:,4:6};

Data = Acc;

Fs = 4500;                                                  % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
NrSp = size(Data,2);                                        % Number Of Subplots
L = size(Data,1);                                           % Length Of Data Vectors
t = linspace(0, Fs, L)/Fs;                                  % Time Vector

NFFT = 2^nextpow2(L);                                       % For Efficiency
FTData = fft(Data - mean(Data),NFFT)/L;                     % Fourier Transform
Fv = linspace(0, 1, NFFT/2-1)*Fn;                           % Frequency Vector (For Plots)
Iv = 1:numel(Fv);                                           % Index Vector (For Plots)

figure
% sp = [1:2:NrSp 2:2:NrSp];                                   % Order 'subplot' Plots
for k = 1:3
    subplot(3,1,k)
    plot(Fv, abs(FTData(Iv,k))*2)                           % Plot Fourier Transforms
    grid
    xlabel('Frequency')
    ylabel('Magnitude')
    xlim([0 Fn/200])
%     xlim([5 10])
    ylim([0 0.4])
    title(sprintf('Column %d',k))
end
sgtitle('Fourier Transform of Data')

figure
% sp = [1:2:NrSp 2:2:NrSp];
for k = 1:NrSp
    subplot(3,1,k)
    plot(t, Data(:,k))
    grid
    ylim([-3 3])
    xlabel('Time')
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
sgtitle('Original Time-Domain Data')

Wp = [2  8]/Fn;                                             % Define Passband In Hz, Normalise To (0,pi)
Ws = Wp.*[0.95 1.04];                                       % Stopband (Normalised)
Rs = 50;                                                    % Stopband Ripple (Attenuation) dB
Rp =  1;                                                    % Passband Ripple dB
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp);                                % Design Filter
[sos,g] = zp2sos(z,p,k);                                    % Implement Filter As Second-Order-Section Representation

figure
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 10])
set(subplot(2,1,2), 'XLim',[0 10])
sgtitle('Filter Bode Plot')

Data_Filt = filtfilt(sos,g,Data);

figure
for k = 1:3
    subplot(3,1,k)
    plot(t, Data_Filt(:,k))
    grid
    ylim([-1 1]*0.75)
    xlabel('Time')
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
sgtitle(sprintf('Time-Domain Plots of Bandpass (%.1f - %.1f Hz) Filtered Data',Wp*Fn))

FTData_Filt = fft(Data_Filt - mean(Data_Filt),NFFT)/L;      % Fourier Transform

figure
% sp = [1:2:NrSp 2:2:NrSp];                                   % Order 'subplot' Plots
for k = 1:3
    subplot(3,1,k)
    plot(Fv, abs(FTData_Filt(Iv,k))*2)                      % Plot Fourier Transforms
    grid
    xlabel('Frequency')
    ylabel('Magnitude')
    xlim([0 Fn/200])
%     xlim([5 10])
    ylim([0 0.4])
    title(sprintf('Column %d',k))
end
sgtitle('Fourier Transform of Filtered Data')

The filter is operating correctly and the Fourier transforms of the filtered data are as expected.

EDIT — (1 Oct 2022 at 17:15)

Added this plot, as requested —

figure
hold on
for k = 1:3
    plot(t, Data_Filt(:,k), 'DisplayName',sprintf('Column %2d',k))
end
hold off
grid
ylim([-1 1]*0.75)
xlabel('Time')
ylabel('Amplitude')
title(sprintf('Time-Domain Plots of Bandpass (%.1f - %.1f Hz) Filtered Data',Wp*Fn))
legend('Location','best')

.

[Paras Shaikh]

Can u pleaseee make one graph of filtered data instead of 3 columns . I want in one graph.. but cant do myself as im getting error in my code

[Star Strider]

Of course!

I added it as an EDIT at the end of my previous Comment.

[Paras Shaikh]

Hey cant we do like when first column is finishing.. other is starting from there.. like i want 3 of them together as one ends 2nd starts from there in horizontal line.. same as my data.. first column is of normal tremor . 2nd is of moderate . 3rd is of severe.. so want 3 of them in one line..

[Star Strider]

That is not the way your data are organised.

That plots all three accelometer axes over their shared time vector. The original data all go from 0 to about 1 second, because you have one second of data in the files provided. They are different axes of the acclerometer, all collected over the same time interval at the same instants in time. Putting them together sequentially would be meaningless.

[Paras Shaikh]

Which bandpass filter is used in this code.. and why ?

[Star Strider]

This one, and because it meets your requirements —

Fs = 4500; % Sampling Frequency

Fn = Fs/2; % Nyquist Frequency

Wp = [2 8]/Fn; % Define Passband In Hz, Normalise To (0,pi)

Ws = Wp.*[0.95 1.04]; % Stopband (Normalised)

Rs = 50; % Stopband Ripple (Attenuation) dB

Rp = 1; % Passband Ripple dB

[n,Wn] = ellipord(Wp,Ws,Rp,Rs); % Calculate Filter Order

[z,p,k] = ellip(n,Rp,Rs,Wp); % Design Filter

[sos,g] = zp2sos(z,p,k); % Implement Filter As Second-Order-Section Representation

figure

freqz(sos, 2^16, Fs) % Filter Bode Plot

set(subplot(2,1,1), 'XLim',[0 10])

set(subplot(2,1,2), 'XLim',[0 10])

sgtitle('Filter Bode Plot')

Data_Filt = filtfilt(sos,g,Data);

.

[Paras Shaikh]

Ok but what is the name of filter... Is there any name of filter which is used here

[Star Strider]

It is called an Elliptic filter.