How may I perform an accelerometric signal analysis with MATLAB?

1 vue (au cours des 30 derniers jours)
Guglielmo Giambartolomei
Guglielmo Giambartolomei le 5 Oct 2016
Goodmorning everyone, I was asked by the head of the research group to perform an analysis of an accelerometric signal and I want to do it with MATLAB. May anyone recommend me a suitable script? Here there is my small and simple code:
clc;
clear all;
file='SPRG_1.xlsx';
x0=xlsread(file);
%%Time domain analysis of channel 1
x1=x0(:,2);
N1=length(x1);
Fs=2400;
T=1/Fs;
t1=linspace(0,N1*T,N1);
plot(t1,x1)
grid on;
title('Hammer (ch. 1)')
xlabel('t (s)')
ylabel('x1(t), Force (N)')
%%Frequency domain analysis of channel 1
% xdft1=fft(x1-mean(x1))/N;
% xdft1=xdft1(1:N/2+1);
% psdx1=(1/(Fs*N1))*abs(xdft1).^2;
% psdx1(2:end-1)=2*psdx1(2:end-1);
% freq1=0:Fs/length(x1):Fs/2;
%
% figure;
% plot(freq1,abs(xdft1))
% grid on
% title('Single-Sided Amplitude Spectrum of x1(t) (Hammer)')
% xlabel('f (Hz)')
% ylabel('X1(f)')
%%Time domain analysis of channel 2
x2=x0(:,3);
N2=length(x2);
Fs=2400;
T=1/Fs;
t2=linspace(0,N2*T,N2);
figure
plot(t2,x2)
grid on
title('Acc. on the sparger (ch. 2)')
xlabel('t (s)')
ylabel('x2(t), Acceleration (g)')
%%Frequency domain analys of channel 2
xdft2=fft(x2-mean(x2))/N2;
xdft2=xdft2(1:N2/2+1);
psdx2=(1/(Fs*N2))*abs(xdft2).^2;
psdx2(2:end-1)=2*psdx2(2:end-1);
freq2=0:Fs/length(x2):Fs/2;
figure;
plot(freq2,abs(xdft2))
grid on
title('Single-Sided Amplitude Spectrum of x2(t) (Hammer)')
xlabel('f (Hz)')
ylabel('X2(f)')
figure;
plot(freq2,10*log10(psdx2))
grid on
title('Periodogram Using FFT of x2(t)')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
%%Plot the Power Spectral Density of the Signal
figure;
[pxx,fx] = pwelch(x2,[],[],[],Fs);
plot(fx,pxx);
grid on
title('Power Spectral Density of x2(t)')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
%%Design a Lowpass IIR Filter
N=7;
Fp=300;
Ap=1;
h=fdesign.lowpass('N,Fp,Ap', N, Fp, Ap, Fs);
d=design(h, 'cheby1');
%%Apply the filter to to Smooth out the Signal
xfilter = filter(d,x2);
%%Overlay the filtered signal on the original signal.
% Filtered signal is delayed
figure;
plot(t2,x2,'b',t2,xfilter,'r');
grid on;
legend({'Original Signal','Filtered Signal'});
%set(gcf,'NumberTitle','Off', 'Name','Filtered Signal vs. Actual Signal');
%%Compare the original signal and delay compensated filtered signal
figure;
xfiltfilt = filtfilt(d.sosMatrix,d.ScaleValues,x2);
plot(t2,x2,t2,xfiltfilt);
grid on
legend({'Original Signal','Actual (filtered and delayed signal)'});
%%Frequency domain analysis of the filtered and realigned signal
xdft22=fft(xfiltfilt-mean(xfiltfilt))/N2;
xdft22=xdft22(1:N2/2+1);
psdx2=(1/(Fs*N2))*abs(xdft2).^2;
psdx2(2:end-1)=2*psdx2(2:end-1);
freq2=0:Fs/length(x2):Fs/2;
figure;
plot(freq2,abs(xdft22))
grid on
title('Single-Sided Amplitude Spectrum of x2(t) filtered and realigned Signal')
xlabel('f (Hz)')
ylabel('X2(f)')
I wanted to perform a pwelch in order to justify the choice of the low-pass band frequency but I'm not able to do that. Have you got any suggestions? Thank you very much!
  5 commentaires
Guglielmo Giambartolomei
Guglielmo Giambartolomei le 6 Oct 2016
Modifié(e) : Guglielmo Giambartolomei le 6 Oct 2016
@Jan Simon I tried to format the posted code several times. I clicked on the {}Code button and I made copy/paste of my code within "if true" and "end" but only the first part has been shown. Thank you. Edit: I tried selecting all the code and then press the code button and it worked!
Guglielmo Giambartolomei
Guglielmo Giambartolomei le 6 Oct 2016
Modifié(e) : Guglielmo Giambartolomei le 6 Oct 2016
@dpb I don't know if performing a pwelch with MATLAB in order to see what is the frequency threshold (where can I apply the low-pass filter) is correct. I tried with the following command (as you can see in the thread posted code):
[pxx,fx] = pwelch(x2,[],[],[],Fs);
plot(fx,pxx);

Connectez-vous pour commenter.

Réponses (0)

Catégories

En savoir plus sur AI for Signals dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by