signal processing of real time signal
5 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
Mathieu NOE
2022 年 7 月 19 日
hello
your code does envelope and fft analysis but there is nothing related to conversion from velocity to displacement conversion. Did you mean to do it in the frequency domain and in the time domain ? (fyi, use cumtrapz for that)
採用された回答
Mathieu NOE
2022 年 7 月 19 日
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 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Spectral Measurements についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!