フィルターのクリア

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 Center および File ExchangeDigital Filter Analysis についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by