Converting accelerometer data to frequency

Hi, I'm trying to take accelerometer data I recorded with an Adafruit BNO05 IMU and generate two graphs: one that plots acc vs. time and one that plots acc vs. freq. I'm attaching a data file for reference. Our current code doesn't seem to be working quite right for the acc vs. freq part. I'm also attaching that code. Any suggestions for how we can fix the code?

%Reads raw data from the xlsfile.

function [acceleration,time] = TremorReader(filename)



%Open the Excel file and extract information

[num text raw ] = xlsread(filename) ;

[r,c] = size(num) ;



%Delete unecessary information

num(2:2:end,:) = [] ;

num(:,2) = [] ;

num(:,5) = [] ;



%Extract time

time = num(:,1) ;

time = time - time(1) ;



%convert from millisecond to seconds

timesec = time./1000 ;



%replace in data again

num(:,1) = timesec ;



%Extract all axes



ax = num(:,2) ;

ay = num(:,3) ;

az = num(:,4) ;



%Compute the cumulative/magnitude vector

acc = sqrt(ax.^2 + ay.^2 + az.^2); %possibly in G's
acc_cmpersec = acc.*980; %converts to cm/s^2


hold on

subplot(1,2,2);

plot(timesec,acc);



dt = .103;%time between samples????

N1 = length(acc);

  N = 2^nextpow2(N1);

  if N > N1

      acc(N1+1:N) = 0; % pad array with 0's

  end

 df = 1 / (N*dt); % frequency increment

 Nyq = 1 / (2*dt); % Nyquist frequency

 freq = fftshift(fft(acc));


freq = sort(freq);
mag_freq = sqrt((real(freq).^2)+(imag(freq).^2)); %magnitude of the frequency, combines the real and imaginary parts 
%M = abs(freq), this is the same as mag_freq, just a different method
phase_freq = angle(freq);
subplot(1,2,1);
plot(mag_freq, acc);
%hold on

% used plot(freq,acc,'r') before (other night)

%hold off


RMSacc = norm(acc_cmpersec)./sqrt(length(acc_cmpersec)); %root mean square to combine the acceleration vectors into one value for the amplitude and score
Dom_freq = max(mag_freq); %taking the dominant frequency 
T = RMSacc./((2.*pi.*Dom_freq ).^2); %translational displacement
amplitude = abs((T)./(sin(2.*pi*Dom_freq .*4)));


%if amplitude==0;
%        score=0;
%    elseif amplitude<1;
%        score=1;
%    elseif amplitude>=1 & amplitude<3;
%        score=2;
%    elseif amplitude>=3 & amplitude<10;
%        score=3;
%    else
%        score=4;
%    end
%disp(score); %These will be used to find TP,TN,FP, and NP

end

Answer 1

[Star Strider]

Try this:

num = xlsread('test550.csv');
num = num(1:2:end,:);                                           % Eliminate ‘NaN’ Values
num(:,1) = num(:,1)*1E-3;                                       % Convert Milliseconds To Seconds
Fs = 1/mean(diff(num(:,1)));                                    % Sampling Frequency (Convert Milliseconds To Seconds)
Fn = Fs/2;                                                      % Nyquist Frequency
Tr = linspace(num(1,1), num(1,end), size(num,1));               % Time Vector (Regular Samples)
Dr = resample(num(:,2:end), Tr);                                % Resample To Constant Sampling Interval
Dr_mc  = Dr - mean(Dr,1);                                       % Subtract Mean

FDr_mc = fft(Dr_mc, [], 1);                                     % Fourier Transform
Fv = linspace(0, 1, fix(size(FDr_mc,1)/2)+1)*Fn;                % Frequency Vector
Iv = 1:numel(Fv);                                               % Index Vector

figure
plot(Fv, abs(FDr_mc(Iv,:))*2)
grid
hl = legend('B','C','D','E','F', 'Location','NE')
title(hl, 'Excel Columns')
xlabel('Frequency (Hz)')
ylabel('Amplitude')

Experiment to get the result you want.

[Samantha Kench]

The graph looks great! What do the letters correspond to (the 'B', 'C', 'D', etc.)? Thank you for helping!

[Star Strider]

As always, my pleasure!

The letters are the columns in your Excel file. They weren’t labeled as anything else, so I didn¹t know what else to call them.

[Samantha Kench]

Ok! Thank you for responding so fast! Your feedback was very helpful.

[Star Strider]

As always, my pleasure!

Also, if you want to increase the frequency resolution, this works:

N1 = size(num,1);
N = 2^(nextpow2(N1)+3);

FDr_mc = fft(Dr_mc, N);                                         % Fourier Transform

The data are unchanged. The spectrum plot looks smoother. The rest of the code is unchanged.

[Samantha Kench]

Quick question: the file I included was tested on a vortexer at 550 rpm, so we figured the dominant frequency should be around 9 Hz (550/60). Based on what we are seeing it looks like the dominant frequency of the y axis is around 0.5 Hz. Could our data be inaccurate, or would we have to merge the axes for dominant frequency? This is for a biomedical engineering project that measures tremor.

[Star Strider]

Long reply :-)

Looking at the time domain of the mean-corrected signals, the waveform with the highest amplitude has a mean period of about 5.15 s, corresponding to a frequency of about 0.2 Hz. The spectrum shows that your data have two principal peaks, with one signal having a peak at about 0.15 Hz and three others at about 0.6 Hz.

The sampling frequency you used (9.6 Hz) cannot measure a 9 Hz signal, since the highest frequency you can uniquely resolve is half that, about 4.9 Hz (the Nyquist frequency). I’m not certain what tremors you want to measure, however a reasonable sampling frequency would be at least 50 Hz, and at best 250 Hz, because that gives you enough frequency resolution to be able to eliminate line frequency noise and other interference. Most tremor frequencies are 15 Hz or less (the frequency depends on the etiology), so anything more than twice that would work as a sampling frequency. Within reason, higher frequencies are usually best, however anything higher than about 250 Hz for EKG or tremor measurement would be excessive.

So the problem with the signal you posted is that you cannot resolve a 9 Hz frequency no matter what you do. You need to sample at a much higher rate for that. Also, use a hardware anti-aliasing filter (a Bessel filter is usually best because it has a linear phase characteristic) with a cutoff at your Nyquist frequency between your amplifiers your ADC. Biomedical engineering instrumentation texts discuss this in detail.

[Samantha Kench]

Thank you for all of your help! We will look into all of this. This has been extremely helpful :)

[Star Strider]

As always, my pleasure!

Have fun with your project.