How to convert a accelerometer data to displacements pdr

Question

Yunes Almansoob 2017-2-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/326493-how-to-convert-a-accelerometer-data-to-displacements-pdr

编辑： Yunes Almansoob 2017-3-1

data.txt

Hi friends, i have accelerometer data PDR i want use fft filter to convert it to displacements i use this tow code but the second code some problem and the result not good i am new in mat lab i don't know how use it,the data attached.

%% Frequency domain integration

%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close all hidden

%%%%%%%%%%%%%%%%%%

fni=input('frequency domain integration - input data file name','s');

fid=fopen(fni,'r');

sf=fscanf(fid,'%f',1);%Sampling frequency

fmin=fscanf(fid,'%f',1);%minimum cutoff frequency

fmax=fscanf(fid,'%f',1);%maximum cutoff frequency

c=fscanf(fid,'%f',1);%unit transform coefficients

it=fscanf(fid,'%f',1);%of the number of points

sx=fscanf(fid,'%s',1);%Scaling of the horizontal axis

sy1=fscanf(fid,'%s',1);%of the vertical axis of the input unit

sy2=fscanf(fid,'%s',1);%of the vertical axis of the output unit of the label

fno=fscanf(fid,'%s',1);%output data file name

x=fscanf(fid,'%f',[1,inf]);%The input data is stored as a row vector

status=fclose(fid);

n=length(x);

%Build time vector

t=0:1/sf:(n-1)/sf;

%Is greater than and closest to n, the power of 2 is the FFT length

nfft=2^nextpow2(n);

%FFT transform

y=fft(x,nfft);

%Calculate frequency interval (Hz / s)

df=sf/nfft;

%Calculates the subscript of the specified frequency band corresponding to the frequency array

ni=round(fmin/df+1);

na=round(fmax/df+1);

%Calculate round frequency interval (rad / s) dw=2*pi*df;

%Create a positive discrete circular frequency vector

w1=0:dw:2*pi*(0.5*sf);

%Create a negative discrete circle frequency vector

w2=-2*pi*(0.5*sf-df):dw:-dw;

%Combines positive and negative round frequency vectors into one vector

w=[w1,w2];

%To the number of points for the index, the establishment of circular frequency variable vector

w=w.^it;

%Of the frequency domain transformation

a=zeros(1,nfft); a(2:nfft-1) =y(2:nfft-1)./w(2:nfft-1);

if it == 2

y=-a; %for the second integral phase transformation

else

a1=imag(a); a2=real(a); y=a1-a2*i; %Phase transformation for one-time integration

end

a=zeros(1,nfft);

% Eliminates the frequency components outside the specified positive frequency band

a(ni:na)=y(ni:na);

%Eliminates frequency components outside the specified negative frequency band a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);

y=ifft(a,nfft); %IFFT transform

%Take the inverse of the real part of the n elements and multiplied by the unit transformation coefficient for the integration results

y=real(y(1:n))*c;

subplot(2,1,1); plot(t,x); xlabel(sx); ylabel(sy1); grid on; %draw a few cents of the time curve graph

subplot(2,1,2); plot(t,y); xlabel(sx); ylabel(sy2); grid on; %plot the time- %Open the file output after the integration of the data

fid=fopen(fno,'w');

for k=1:n, fprintf(fid,'%f \n',y(k)); end

status=fclose(fid);

the second code % Frequency domain integration

%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close all hidden

%%%%%%%%%%%%%%%%%%%

fni = input ('frequency domain integration - input data file name:', 's');

fid = fopen (fni, 'r');

sf = 12000;% sampling frequency

fmin = 0.1;% minimum cutoff frequency

fmax = 6000;% Maximum cutoff frequency

c = 1;% unit transform coefficient

it = 2;% integration times

sx = 'Time (s)';% Dimension of the horizontal axis

sy1 = 'Acceleration (m / s ^ 2)';% Vertical axis Enter the label for the unit

sy2 = 'Displacement (m)';% Vertical coordinate axis output unit of the annotation

out.txt;% Output data file name

% Input data is stored as a row vector. X = fscanf (fid, '% f', [1, inf]);

% Acceleration Time Data Separation

for i = 1: 1: (length (x) / 2)

% Time data

t (1i) = x (2 * 1i-1);

% Acceleration data

xx (1i) = x (2 * 1i);

end

status = fclose (fid);

n = length (xx);

The power of% greater than and closest to n is the power of the FFT length

nfft = 2 ^ nextpow2 (n);

% FFT transform

y = fft (xx, nfft);

% Calculated frequency interval (Hz / s)

df = sf / nfft;

% Calculates the index of the frequency array for the specified frequency band

ni = round (fmin / df + 1);

na = round (fmax / df + 1);

% Calculate the circular frequency interval (rad / s)

dw = 2 * pi * df;

% Establish a positive discrete circular frequency vector

w1 = 0: dw: 2 * pi * (0.5 * sf-df);

% Creates a negative discrete circular frequency vector

w2 = 2 * pi * (0.5 * sf-df): -dw: 0;

% Combines positive and negative circular frequency vectors into a single vector

w = [w1, w2];

% The number of integration times is used as the index to construct the circular frequency variable vector

w = w.^it;

% To carry out integral frequency domain transformation

a = zeros (1, nfft); a (2: nfft-1) = y (2: nfft-1) ./w (2: nfft-1);

if it == 2

y = -a;% Phase transformation for quadratic integration

else

a1=imag(a); a2=real(a); y=a1-a2*i;% phase integration for one integration end

a = zeros (1, nfft);

% Eliminates frequency components outside the specified positive band

a (ni: na) = y (ni: na);

% Eliminates frequency components outside the specified negative band

a (nfft-na + 1: nfft-ni + 1) = y (nfft-na + 1: nfft-ni + 1);

y = ifft (a, nfft);% IFFT transform

% The inverse of the real part of the n elements and multiplied by the unit transform coefficient for the integral result

y = real (y (1: n)) * c;

subplot (2,1,1); plot (t, xx); xlabel (sx); ylabel (sy1); grid on;% plotting time- % On drawing the integrated time-history curve graph% open the file after the output of the integrated data (t, y); xlabel (sx); ylabel (sy2)

fid = fopen (fno, 'w');

for k = 1: n, fprintf (fid, '% f \ n', y (k)); end

status = fclose (fid);