필터 지우기
필터 지우기

How may I perform an accelerometric signal analysis with MATLAB?

조회 수: 2 (최근 30일)
Guglielmo Giambartolomei
Guglielmo Giambartolomei 2016년 10월 5일
편집: Guglielmo Giambartolomei 2016년 10월 6일
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
Guglielmo Giambartolomei
Guglielmo Giambartolomei 2016년 10월 6일
편집: Guglielmo Giambartolomei 2016년 10월 6일
@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 2016년 10월 6일
편집: Guglielmo Giambartolomei 2016년 10월 6일
@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);

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Digital Filter Analysis에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by