How to convert a accelerometer data to displacements pdr
5 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!