How to convert a accelerometer data to displacements pdr
5 views (last 30 days)
Show older comments
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);
0 Comments
Answers (0)
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!