signal processing of real time signal

3 views (last 30 days)
i m doing signal processing of real time signal i have written code for it but i dont known weather it is correct or not ? i had velocity data and i coverted that into the displacement.please check.
clc;
clear;
close all;
%%Time Domain Plot
vib1=readtimetable("sing2.txt","SampleRate",1);
tiledlayout(4,1)
nexttile
plot(vib1.Time,vib1.Sing1)
L=length(vib1.Sing1);
N= 2^nextpow2(L);
xlabel("Time(s)")
ylabel("Ampl(mm/s^2)")
%%Envolop analysis
nexttile
yh=hilbert(vib1.Sing1,N);
yz=abs(yh);
plot(vib1.Time,yz)
% % FFt
fs=1;%Sampling Frequency
delt=seconds(1/fs); %time step
Totaltime=(length(yz)-1);
my_fft=2/N*(fft(yz));%FFT of Signal normaization
abs_fft=abs(my_fft); %Absolute value of fft
delf=1/Totaltime; %Frequency Resolution
n2=0:1:N/2-1; %fft results are plotted for N/2 data points
fk=delf*n2;% frequency values
abs_fft(1:N/2);
nexttile
plot(fk,abs_fft(1:N/2))
[v,p]=findpeaks(abs_fft(1:N/2),fk,'MinPeakHeight',0.01);
findpeaks(abs_fft(1:N/2),fk,'MinPeakHeight',0.01)
xlabel("Frequency (Hz)")
ylabel("Ampl(mm/s^2)")
% Spectral density
nexttile
pspectrum(yz,vib1.Time,"spectrogram","FrequencyLimits",[.1005 .2976])
  2 Comments
GOURANG UPADHYAY
GOURANG UPADHYAY on 19 Jul 2022
i have done that also but in diffrent code.
so yes i have converted my data into displacement.but i just want to know approch is correct or not

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 19 Jul 2022
hello again !!
FYI this is simple code for time domain acceleration to velocity / displacement integration (data file attached)
in your case you would only need one (and not two) stage of integration.
be smart and use the FFT spectrum of the input data to fine tune the high pass filter (if needed ! it can be useful for dynamic signals where we need to keep all signals zero mean - even after integration; if this does not apply to your case then remove the detrend and the high pass filter and keep only the integration with cumtrapz)
clc
clearvars
data = readmatrix('Earthquake.xlsx');
t = data(:,1);
dt = mean(diff(t)); % Average dt
fs = 1/dt; % sampling rate
% Acceleration data
acc = data(:,2)*9.81; % Add gravity factor ( g to m/s²)
acc = detrend(acc); % baseline correction
% some additionnal high pass filtering % baseline correction
N = 2;
fc = 0.25; % Hz
[B,A] = butter(N,2*fc/fs,'high');
acc = filter(B,A,acc);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
velocity = cumtrapz(dt,acc);
velocity = detrend(velocity); % baseline correction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp = cumtrapz(dt,velocity);
disp = detrend(disp); % baseline correction
%
figure(1)
subplot(311),plot(t,acc);
title('accel'); ylabel('m/s²');xlabel('time (s)');
subplot(312),plot(t,velocity);
title('velocity'); ylabel('m/s');xlabel('time (s)');
subplot(313),plot(t,disp);
title('disp'); ylabel('m');xlabel('time (s)');
%%fft
nfft = length(acc);
fft_spectrum = abs(fft(disp,nfft));
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*fs/nfft;
figure(2)
plot(freq_vector,fft_spectrum);
title('disp'); ylabel('m');xlabel('Frequency (Hz)');
xlim([0 10]);
  8 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by