%% ========================================================================
%  DIAGNOSTYKA ŁOŻYSKA 6201, UC201 i 6205-2RS JEM SKF
%  ------------------------------------------------------------------------
%  Skrypt:
%    1.  Wczytuje sygnał przyspieszeń (100 s, fs = 48 kHz) z pliku
%        „dane_Lud3_1.mat" lub „dane_Lud3_2.mat" 
%    2.  Filtruje sygnał FIR‑em dolnoprzepustowym 0‑8 kHz
%        (impulsy łożyskowe najczęściej leżą w tym paśmie).
%    3.  Liczy widmo (FFT) oraz zaznacza:
%          •  linię prędkości obrotowej (fo) i jej harmoniczne,
%          •  częstotliwości charakterystyczne BPFO, BPFI, BS, FT 
%             wraz z harmonicznymi.
%    4.  Wyznacza obwiednię (metoda RMS w oknie WL = 10 próbek) i rysuje
%        widmo obwiedni tuż nad widmem sygnału, aby łatwiej ocenić linie.
%  ------------------------------------------------------------------------
%  UWAGA:  Skrypt ilustruje zasadę – nie optymalizuje parametrów filtracji
%  ani okna obwiedni; w praktycznej diagnostyce warto je dostroić.
%  ========================================================================
close all;  clear;  clc;

%% 0.  PODSTAWOWE PARAMETRY
fp  = 48e3;      % [Hz]  – częstotliwość próbkowania
tk  = 100;       % [s]   – całkowity czas sygnału

% % dla eksperymentu 1
yup = 0.0005;    % przesunięcie wykresu
xlimitob = 100;
ylimitob = 1e-3;

% %-----------Eksperyment 1, łożysko 6201 - 6 kulek ----------------
fo = 8.6975; % [Hz]  – prędkość obrotowa wału
BPFOq=2.143;
BPFIq=3.857;
BSq=1.608;
FTq=0.357;
load dane_Lud3_1                  
% 
% %--------------------------------------------------------

% % Do eksperymentów 2 i 3-----------------
% fo  = ;
% n  = ;        % liczba elementów tocznych            [szt.]
% d  = ;        % średnica kulki                       [m]
% D  = ;        % średnica podziałowa                  [m]
% alpha = ;     % kąt kontaktu (poprzeczne – 0 rad)    [rad]
%               
% 
% BPFOq = (n/2) * (1 - (d/D)*cosd(alpha));           % bieżnia zewn.
% BPFIq = (n/2)  * (1 + (d/D)*cosd(alpha));           % bieżnia wewn.
% BSq  = (D/(2*d)) * (1 - ((d/D)*cosd(alpha))^2);   % obrót kulki
% FTq  = 0.5 * (1 - (d/D)*cosd(alpha));             % koszyk
%
% yup = 0.003;
% xlimitob = 600;
% ylimitob = 6e-3;
%%-------------------------------------------------------------------------------

% Częstotliwości charakterystyczne łożyska 
BPFO = BPFOq*fo;  % bieżnia zewn.
BPFI = BPFIq*fo;  % bieżnia wewn.
BS   = BSq*fo;  % element toczny
FT   = FTq*fo;  % koszyk

%% 1.  WCZYTANIE SYGNAŁU --------------------------------------------------
z_raw = yt(:,2);                 % dokładniejszy kanał przyspieszeń
t     = yt(:,3);                 % wektor czasu z pliku (≈ linspace)

%% 2.  FILTRACJA DOLNOPRZEPUSTOWA 0–8 kHz (liniowa faza) ------------------
N   = 100;   Fp = 8e3;   Rp = 0.01;  % rząd, f_odcięcia, zafalowanie
lpFilt = dsp.LowpassFilter('DesignForMinimumOrder',false,'FilterOrder',N, 'PassbandFrequency',Fp, 'SampleRate',fp,'PassbandRipple',Rp,'StopbandAttenuation',80);
z = lpFilt(z_raw);               % zero‑phase filtering (dsp obj. = linear)

figure;                           % szybki podgląd filtracji
plot(t(1:5000),z_raw(1:5000)); hold on;
plot(t(1:5000),z(1:5000));
legend('surowy','po LP 0‑8 kHz'); grid;
title('Fragment sygnału przed i po filtracji');

%% 3.  WIDMO SYGNAŁU  -----------------------------------------------------
L = 2^23;                         % punkty FFT (≈ 8.4 mln → Δf ≈ 5.7 mHz)
Zspec = abs( fft(z, L) ) / (L/2); % moduł FFT z normalizacją amplitudy
f  = (0:L/2-1).' * fp/L;          % wektor częstotliwości (0 → fp/2)

figure('Name','Widmo sygnału przyspieszeń');
plot(f, Zspec(1:L/2), 'LineWidth',2); hold on; grid on;
xlabel('częstotliwość [Hz]'); ylabel('|X(f)|  [a.u.]');
xlim([0 xlimitob]);                   % oglądamy tylko do 100 Hz

% -- linie prędkości obrotowej (1× … 6×) -----------------
for k = 1:6
    xline(k*fo, 'k--', '1X','LabelOrientation','horizontal');
end

% -- linie łożyskowe + harmoniczne -----------------------
drawLines = @(base,col,label) arrayfun(@(h) xline(h*base, ['--' col], sprintf('%s·%d',label,h),'LabelOrientation','horizontal'), 1:6); % funkcja pomocnicza

drawLines(BPFO,'g','BPFO');
drawLines(BPFI,'m','BPFI');
drawLines(BS  ,'b','BS');
drawLines(FT  ,'c','FT');

%% 4.  OBWIEDNIA (RMS, okno WL)  +  WIDMO OBWIEDNI ------------------------
WL = 10;                                     % okno 10 próbek WL (~0.21 ms)
[yupper,~] = envelope(abs(z), WL, 'rms');    % górna obwiednia RMS
Zenv = abs( fft(yupper, L) ) / (L/2);        % widmo obwiedni

% dorysuj widmo obwiedni o yup wyżej, by nie zasłonić widma bazowego
plot(f, Zenv(1:L/2) + yup, 'r', 'LineWidth',1);
ylim([0 ylimitob]);
f=get(gca,'Children');
legend([f(32),f(1)],'widmo przyspieszeń','widmo obwiedni')
set(gcf,'position',[0,150,1900,800])


%% 5.  KRÓTKI WYKRES CZASOWY (SYGNAŁ + OBWIEDNIA) -------------------------
figure('Name','10 ms sygnału + obwiednia');
idx = 1:round(0.01*fp);                 % pierwsze 10 ms
plot(t(idx), z(idx)); hold on;
plot(t(idx), yupper(idx), 'k', 'LineWidth',1);
xlabel('czas [s]'); ylabel('przyspieszenie  /  obwiednia');
legend('sygnał filtrowany','obwiednia RMS'); grid;