34 views (last 30 days)

Show older comments

Hi i need matlab code for frequency weighted acceleration as ISO 2631-1 plz.

Attached: sample data sheet (x, y, z axes m/s^2), sample rate 0.005 sec.

William Rose
on 16 Aug 2021

The help you are requesting requires a lot more than the usual amount of effort for answers on this site.

I downloaded the vibrationdata package (1_1_2021 vrsion) from Tom Irvine's site. When I unzipped it it was over 100 MB. I ran vibrationdata in Matlab. The GUI opens. You said "i get this warrning (nnn>50000) when press on calculater botton". However, I do not see a calculator button. I see the screen below.

I assume you loaded some of your own data, and then selected various options, or clicked on some the buttons and menus on the screen. But I don;t know what you did, and I do not see a calcualtor button.

I think that if you have quesitons about this package, you should ask its author, Mr. Tom Irvine.

I tried the link you posted for ISO 2631-1. The link goes to a long list of ISO standards, one of which is ISO 2631-1:1997. I clicked on that link. I would need to pay 138 CHF to see the standard. If you have access to the standard, I would check the following: Does the standard recommend a specific sampling rate, or a range for sampling rate? Does the standard recommend a specific sample duration, or a range for sampling duration? Does the standard recommend adding up all the powers in thespectrum over a certain range of frequencies, with different wieghts for different frequencies? If so, what is the range? What are the weights? Acceleration is a vector. How does the standard reocmmend dealing with the vectorial nature of acceleration? For example, does the standard recommend measuring x, y, and z components of acceleration separately, and analyzing each separately? Or does the standard recommend measuring the instantaneous magnitude, and ignoring the direction? Or some other approach?

Do you have a file of time-domain vibration data that you can post? If so, please epxlain the organization of the file: identity and units of measure for each column, and sampling rate, if there is not a time column.

William Rose
on 16 Aug 2021

Thank you to @Simon Chan for reading your quesiotn carefully and for noticing that you had posted a data file and that you have posted the smapling rate. Sinc you are new to Matlab, I will describe some steps to import and view your data and prepare it for analysis in the 'vibrationdata' package.

accelData=importdata('acceleration.csv');

disp(accelData);

data: [1155×3 double]

textdata: {'X' 'Y' 'Z'}

colheaders: {'X' 'Y' 'Z'}

Since acceleration.csv includes text and data, Matlab imports it into a "structure", which can contain more than one type of data. I used the "disp" command to see what is inside "accelData". I learned from the reuslts of the disp command that there is a 1155x3 array of double precision numbers inside accelData, and the array name is accelData.data.

dt=0.005; %time step size

t=0:dt:1154*dt; %create a vector of times

accelX=accelData.data(:,1); %X component of acceleration=column 1 of array accelData.data

accelY=accelData.data(:,2); %Y component=column 2 of accelData.data

accelZ=accelData.data(:,3); %Z component=column 3

accelR=sqrt(accelX.^2+accelY.^2+accelZ.^2); %compute resultant acceleration

figure; %open a blank figure window

subplot(2,1,1); %divide figure into 2 rows by 1 column; prepare to plot in top panel

plot(t,accelX,'r',t,accelY,'g',t,accelZ,'b'); %plot, X, Y, Z versus time

grid on; legend('X','Y','Z'); %add grid and legend

ylabel('Accel (m/s^2)'); %add axis label

subplot(2,1,2); %divide figure into 2 rows by 1 column; prepare to plot in bottom panel

plot(t,accelR,'k'); %plot resultant acceleration versus time

grid on; legend('Resultant'); %add grid and legend

xlabel('Time (s)'); ylabel('Accel (m/s^2)'); %add axis labels

The commands above generate the figure below.

I learned yesterday that 'vibrationdata' accepts an array that has a time column and a data column. Therfore I create two such arrays, for X acceleration and resultant acceleration:

dataX=[t' accelX]; %transpose t, to make it a column vector, to match accelX

dataR=[t' accelR];

Then I run 'vibrationdata'

vibrationdata

I make the selections shown below and click "Begin Signal Analysis".

In the next pop-up window, I select "ISO 2631 Human Vibration", and I click "Analyze".

In the next pop-up window, I select Units=m/s^2; Mean Removal=Yes; Weight=Wk, Z-axis, Seat; Subdivide into Segments=No. I specify input array name=data1, start time=0, end time=6. I click "Calculate". This produces a plot like the one shown by @Simon Chan:

You can repeat for the resultant aceleration, and so on. As Mr. Chan correctly notes, the 'vibrationdata' package has many options.

William Rose
on 16 Aug 2021

Edited: William Rose
on 16 Aug 2021

I decided to test the vibrationdata package of Mr. Tom Irvine by creating some simulated acceleration data. My data is random noise with a white (i.e. flat) spectrum frequency range from 1 to 20 Hz. The simulated data is recorded at 160 Hz for 1 minute. Here is a script to make the simulated data with the specified properties:

clear;

fs=160; %sampling rate (Hz)

Td=60; %duration (s)

N=Td*fs; %number of samples

dt=1/fs; %sampling interval

df=1/Td; %frequency spacing of the signal, in the frequency domain

t=0:dt:Td-dt; %vector of times

f=0:df:fs-df; %vector of frequencies for the discrete Fourier transform

Xf=zeros(1,N); %initialize array Xf(), for the signal's discrete Fourier transform

%next line makes the spectrum have magnitude 1, and random phase, from 1 to 20 Hz

Xf(f>=1 & f<=20)=ones(1,length(f(f>=1 & f<=20))).*exp(1i*2*pi*rand(1,length(f(f>=1 & f<=20))));

%Next line makes the spectrum from 140 to 159 Hz be the complex conjugate of the spectrum from 20 to 1 Hz.

%This is required by the mathematics of the discrete Fourier transform, in order to obtain a time domain signal that is real and not complex.

Xf(f>=140 & f<=159)=conj(flip(Xf(f>=1 & f<=20)));

Plot the magnitude and phase of the signal's DFT to confirm that it has the desired properties: constant magnitude, and random phase angles, from 1 to 20 Hz, and Xf(f=140 to 159 Hz) is complex conjugate of X(f=20 to 1 Hz).

figure;

subplot(2,1,1); plot(f, abs(Xf),'.r'); ylabel('Magnitude');

subplot(2,1,2); plot(f, angle(Xf),'.r'); ylabel('Phase Angle (radians)'); xlabel('Frequency (Hz)');

Compute the time domain signal by inverse FFT:

xt=ifft(Xf);

Plot the entire time domain signal, and the first and last 1 seconds:

figure;

subplot(2,1,1); plot(t,xt,'.b'); xlabel('Time (s)'); ylabel('Amplitude');

subplot(2,2,3); plot(t(t<1),xt(t<1),'.-b'); xlabel('Time (s)');

subplot(2,2,4); plot(t(t>=59),xt(t>=59),'.-b'); xlabel('Time (s)');

The time domain signal xt looks reasonable. It is real and randm-looking and there are not fluctuaitons with frequencies that are obviousy higher than 20 Hz. Save it in a text file so we can read it in to Tom Irvine's 'vibrationdata' package.

save('simVibData1.txt','xt','-ascii');

William Rose
on 16 Aug 2021

I tried using the file generated by the code above in'vibriton data'. After clicking sevral buttons I reached a screen that says the data file should have two columns: time and acceleration. And after experimenting, I found that it 'vibraitondata' looks for an array in memory, rather than a file on disk. Therefore I am changing the last line of hte code I posted above, from

save('simVibData1.txt','xt','-ascii');

to

data1=[t',xt'];

Then I make selections in 'vibrationdata': Under "Signal Analysis Functions", I select "Time History" and "Acceleration" and "ISO Generic VC, 2631, 10816'. Then I click the button "Begin Signal Analysis". A new window pops up. In that window, I select "ISO 2631 Human Vibration". Then I click "Analyze". A new window pops up. In that window, I enter Input Array Name=data1 (because I have alrady created that array with my script), Unit=G, Mean Removal=Yes, start time=0 sec, end time=60 sec, Weight='Wk, Z-axis, seat', Segment Duration=10 sec, Percent Overlap=50. Then I click the button "Calculate". Three figures pop up see below, and a window pops up that says "Results written to command window.". I click the "OK" button in that window and it disappears.

I had the same result as you ("Warning: nnn>500000") when I left start time, end time, and duation fields blank, and then I clicked "Calculate".

William Rose
on 16 Aug 2021

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!