Power spectral analyses on accelerometer data

Hi there,

I was hoping to get some clarification. I have 3D acceleration data in the x, y, and z-room oriented directions. Data attached (column 1: time; column 2: x-axis; column 3: z-axis; column 4: y-axis)

Below is a description of what I’m trying to achieve:

1. Perform FFT
2. Perform power spectral density using the periodogram function
3. Perform power spectral density using the periodogram function to obtain the peak frequency between 3-15Hz using a sliding window of 1-s over each 3-s window with 50% overlap

I have bolded the areas in which the codes are not working. Is anyone able to identify my errors in the code below or clarify how to correctly perform the psd (psd performed on the FFT, and using a sliding window)?

t = data(:,1);                                           % In seconds
L = size(data,1);                                        % Data time vector length
Accel = data(:,2:end);                                   % Accel data
Fs = 5000;                                               % Sampling frequency (Hz)
Fn = Fs/2;                                               % Nyquist 

%FFT
FT_Accel = fft(Accel)/L;                               
FT_Accel_ampl = 2*abs(fft(Accel)/L);                     
FT_Accel_power = FT_Accel_ampl.^2;

%PSD
Ps = spectrum.periodogram;
psdest = psd(Ps, Accel,'Fs',Fs);
semilogx(psdest.Frequencies,10*log10(psdest.Data));
grid on;

%PSD - calculating peak frequency between 3-15Hz using a sliding window of 1 s over each 3 s window with 50% overlap
N = size(data,1);
window = 3*Fs;                                                   %3s window
[pxx,f] = periodogram(Accel,window,[],Fs,hanning(300000));       %is this not working because the '5000' Fs needs to be replaced with the nyquist? Default overlap of 50%
plot                                                                

Answer 1

[Mathieu NOE]

suggested code below

you can reduce your sampling rate. there is not much signal above 100 Hz , so sampling at 500 Hz is enough

(factor 10 reduction on file size and post processing time)

enjoy !

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096*4;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting : dB(L))
% option_w = 1 : A weighted spectrum  : (dB (A) ) for acoustic measurements
option_w = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[data,head] = readclm('Accel.txt',4,1);
channel = 2;  % 1 = X, 2 = Y , 3 = Z
signal = data(:,1+channel);

samples = length(signal);

Fs = 5000 ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(1),semilogx(freq,sensor_spectrum_dB);grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,w,NOVERLAP);  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end

% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');


function pondA_dB = pondA_function(f)
% dB (A) weighting curve
n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
pondA = n./r;
pondA_dB = 20*log10(pondA(:));
end

function  [outdata,head] = readclm(filename,nclm,skip,formt)

% READCLM Reads numerical data from a text file into a matrix.
%	Text file can begin with a header or comment block.
%	[DATA,HEAD] = READCLM(FILENAME,NCLM,SKIP,FORMAT)
%	Opens file FILENAME, skips first several lines specified
%	by SKIP number or beginning with comment '%'.
%	Then reads next several lines into a string matrix HEAD
%	until the first line with numerical data is encountered
%	(that is until first non-empty output of SSCANF).
%	Then reads the rest of the file into a numerical matrix
%	DATA in a format FORMAT with number of columns equal
%	to number of columns of the text file or specified by
%	number NCLM. If data does not match the size of the
%	matrix DATA, it is padded with NaN at the end.
%
%	READCLM(FILENAME) reads data from a text file FILENAME,
%	skipping only commented lines. It determines number of
%	columns by the length of the first data line and uses
%	the floating point format '%g';
%
%	READCLM uses FGETS to read the first lines and 	FSCANF
%	for reading data.

%  Kirill K. Pankratov,  [e-mail]
%  03/12/94, 01/10/95.

 % Defaults and  parameters ..............................
formt_dflt = '%g';  % Default format for fscanf
addn = nan;         % Number to fill the end if necessary

 % Handle input ..........................................
if nargin<1, error('  File name is undefined'); end
if nargin<4, formt = formt_dflt; end
if nargin<3, skip = 0; end
if nargin<2, nclm = 0; end
if isempty(nclm), nclm = 0; end
if isempty(skip), skip = 0; end

 % Open file ............................
[fid,msg] = fopen(filename);
if fid<0, disp(msg), return, end

 % Find header and first  data line ......................
is_head = 1;
jl = 0;
head = ' ';
while is_head  % Add lines to header.....
  s = fgets(fid);           % Get next line
  jl = jl+1;
  is_skip = jl<=skip;
  is_skip = jl<=skip | s(1)=='%';
  out1 = sscanf(s,formt);   % Try to read this line
   % If unreadable by SSCANF or skip, add to header
  is_head = isempty(out1) | is_skip;
  if is_head & ~is_skip
    head = str2mat(head,s(1:length(s)-1)); end
end
head = head(2:size(head,1),:);

 % Determine number of columns if not specified
out1 = out1(:)';
l1 = length(out1);
if ~nclm, nclm = l1; end

 % Read the rest of the file ..............................
if l1~=nclm  % First line format is different from ncolumns

  outdata = fscanf(fid,formt);
  lout = length(outdata)+l1;
  ncu = ceil(lout/nclm);
  lz = nclm*ncu-lout;
  outdata = [out1'; outdata(:); ones(lz,1)*addn];
  outdata = reshape(outdata,nclm,ncu)';

else              % Regular case

  outdata = fscanf(fid,formt,[nclm inf]);
  outdata = [out1; outdata'];  % Add the first line

end

fclose (fid);     % Close file ..........
end

[KR]

Thank you for sharing such an informative code, Mathieu!

I am running into a couple of errors when I evaluate the:

function pondA_dB = pondA_function(f)
% dB (A) weighting curve
n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
pondA = n./r;
pondA_dB = 20*log10(pondA(:));
end

and

function  [outdata,head] = readclm(filename,nclm,skip,formt)

Error states:

 function pondA_dB = pondA_function(f)
 ↑
Error: Function definition not supported in this context. Create functions in code file.

I have saved pondA_function and readclm as separate files, and use to call on them:

addpath('E:\Matlab');

Can you identify what I am doing wrong?

Many thanks in advance!

[Mathieu NOE]

hello

glad it helps you

I tested again the code in my release 2020b and it's working fine. I cannot figure what's the root cause of your errors

but if saving the functions as separate files is solving the problem, let's keep it this way

FYI, that's the way I do with those long functions scripts; otherwise it makes the main script very long and unreadable

[KR]

The issue might be due to me using 2019b. Thanks for your help!

[Mathieu NOE]

you're welcome !