%ПРОГРАММА ДЛЯ ОТРАБОТКИ АЛГОРИТМОВ ВИБРОДИАГНОСТИКИ ПОДШИПНИКОВ

fdata=input('Введите имЯ файла замеров -','s')
load (fdata); %Загрузка данных из файла dataL
%Константы
Fs=30000;     %Частота дискретизации равна 30кГц
Ff=16384;       %Частота Фурье-преобразования
Mvp=0.0064;     %Масшт коэффициент 
fp=16.66667;    %Частота вращениЯ по умолчанию
%Параметры подшипника 307
%fp=input('Введите частоту вращениЯ ротора (Гц) fp= ')
frol=0.376*fp; fcep=0.376*fp;  fvn=4.37*fp;  ftk=1.888*fp;  fnar=2.63*fp; 

f0p=fp-fp*0.08;f1p=fp+fp*0.08;              N0p=round(Ff*f0p/Fs);     N1p=round(Ff*f1p/Fs);
f0rol=frol-frol*0.08; f1rol=frol+frol*0.08; N0rol=round(Ff*f0rol/Fs); N1rol=round(Ff*f1rol/Fs);
f0cep=fcep-fcep*0.08; f1cep=fcep+fcep*0.08; N0cep=round(Ff*f0cep/Fs); N1cep=round(Ff*f1cep/Fs);
f0vn=fvn-fvn*0.08; f1vn=fvn+fvn*0.08;       N0vn=round(Ff*f0vn/Fs);   N1vn=round(Ff*f1vn/Fs);
f0tk=ftk-ftk*0.08; f1tk=ftk+0.08;           N0tk=round(Ff*f0tk/Fs);   N1tk=round(Ff*f1tk/Fs);
f0nar=fnar-fnar*0.08; f1nar=fnar+fnar*0.08; N0nar=round(Ff*f0nar/Fs); N1nar=round(Ff*f1nar/Fs);

Nt0p=Fs/f0p;       Nt1p=Fs/f1p;
Nt0rol=Fs/f0rol;   Nt1rol=Fs/f1rol;
Nt0cep=Fs/f0cep;   Nt1cep=Fs/f1cep;
Nt0vn=Fs/f0vn;     Nt1vn=Fs/f1cep;
Nt0tk=Fs/f0tk;     Nt1tk=Fs/f1tk;
Nt0nar=Fs/f0nar;   Nt1nar=Fs/f1nar;

AA=DspData*Mvp; %Чтение данных в вектор АА с масштабированием
sA=size(AA);    %Длинна вектора
t0=sA/Fs;       %Длительность сигнала
t=0:1/Fs:(t0-1/Fs); %Дискретное времЯ равное t0 секунд
me=mean(AA);    %ОпределЯем среднее значение
A0=AA-me;       %Приводим данные к нулю
f=Fs*(0:1000)/Ff; %Нормируем частоту
%ПРОЕКТИРУЕМ ФВЧ
Wp=6*2/Fs;                    %Частота пропусканиЯ 5 Гц
Ws=2*2/Fs; Rs=3;Rp=1;         %Частота задерживаниЯ 2Гц, неравномерность в полосах
[n,Wn]=buttord(Wp,Ws,Rp,Rs);   %Расчет порЯдка фильтра Баттерворда
[bh,ah]=butter(n,Wn,'high');      %Расчет фильтра
[H,w]=freqz(bh,ah,1024);            
A01=filter(bh,ah,A0); %Фильтруем исходный сигнал ФВЧ

%ППОЕКТИРУЕМ ФНЧ
Wp=13500*2/Fs;                    %Частота пропусканиЯ 13500 Гц
Ws=13600*2/Fs; Rs=4;Rp=1;         %Частота задерживаниЯ 13600Гц, неравномерность в полосах
[n,Wn]=buttord(Wp,Ws,Rp,Rs);   %Расчет порЯдка фильтра Баттерворда
[bl,al]=butter(n,Wn);              %Расчет фильтра
%График АЧХ фильтра
[H1,w]=freqz(bl,al,1024);
%figure (1)%Создаем отдельное графческое окно
%subplot (2,1,1), plot(w*Fs/(2*pi),abs(H1),'--b'),grid
%title('АЧХ ФВЧ')
%xlabel('f, Гц')
%ylabel('А, дБ')
%subplot (2,1,2), plot(w*Fs/(2*pi),abs(H),'-.b'),grid
%title('АЧХ ФНЧ')
%xlabel('f, Гц')
%ylabel('А, дБ')
Xa=filter(bl,al,A01); %Фильтруем исходный сигнал ФНЧ

%Проектируем режекторный фильтр длЯ удалениЯ электромагнитных помех сети 50 Гц и 100 Гц
%Wn=[48*2/Fs 52*2/Fs 98*2/Fs 102*2/Fs];
%b=fir1(1000,Wn,'DC-1','scale');
%figure(1)
%[HH,w]=freqz(b,1,1024);
%subplot (3,1,2),plot(w*Fs/(2*pi),abs(HH)),grid
%title('Режекторный фильтр на 50 и 100 Гц')
%xlabel('f, Гц')
%ylabel('А, дБ')
%[Xa,gg]=latcfilt(b,Xa); %Фильтруем исходный сигнал режекторным фильтром
%t=t';
%figure (2)%Создаем отдельное графческое окно 
%plot(t(1:1000),Xa(1:1000),'k'),grid,
%title('Сигнал с датчика вибрации после фильтрации')
%xlabel('t, сек')
%ylabel('А, м/с^2')

%sound(Xa,30000)   %Прослушать звук 

%Интегрируем сигнал
Xv0=cumtrapz(Xa); %ВычислЯем интеграл
m2=mean(Xv0);
Xv=Xv0-m2;      %Нормируем сигнал
Xs0=cumtrapz(Xv); %Расчитываем вибросмещение (интеграл от виброскорости)
m3=mean(Xs0);
Xs=Xs0-m3;

%БЛОК РАССЧЕТА СКЗ
%Проектируем режекторный фильтр 10-1000 Гц
Wn=[5*2/Fs 1050*2/Fs];
b=fir1(3000,Wn,'DC-0','scale');
[Xvskz,gg]=latcfilt(b,Xa); %Фильтруем исходный сигнал 
skz0=sqrt((sum(Xvskz.^2))/120000);

%Преобразование Фурье длЯ виброускорениЯ
Ya=fft(Xa,Ff);  
Fyya=Ya.*conj(Ya)/Ff;    %Амплитуды спектра
f1=Fs*(0:7500)/Ff;   %Частоты спектра
%Преобразование Фурье длЯ виброскорости
Yv=fft(Xvskz,Ff);  
Fyyv=Yv.*conj(Yv)/Ff;    %Амплитуды спектра
f=Fs*(0:540)/Ff;   %Частоты спектра

figure,%Создаем отдельное графческое окно
abYa=abs(Ya)*2/Ff;
plot (f1, abYa(1:7501),'k'),grid
title('Спектр виброускорениЯ')
xlabel('f, Гц')
ylabel('Sa, мм/с^2')
%figure,%Создаем отдельное графческое окно
%plot (f1, Fyya(1:7501),'k'),grid
%title('Спектр мощности виброускорениЯ')
%xlabel('f, Гц')
%ylabel('Sa')
%figure,%Создаем отдельное графческое окно
%plot (f, Fyyv(1:541),'k'),grid
%title('Спектр мощности виброскорости')
%xlabel('f, Гц')
%ylabel('Sv')
figure,%Создаем отдельное графческое окно
abYv=abs(Yv)*2/Ff;
plot (f, abYv(1:541),'k'),grid
title('Спектр виброскорости')
xlabel('f, Гц')
ylabel('Sv, мм/с')

cpektro=input('Строить спектрограмму? 1-да; 0-нет :');
if cpektro>0;
figure      %Создаем отдельное графческое окно
specgram(Xv,[],Fs)%Построение спектрограммы с частотой дисктетизации Fs
colorbar
title ('Спектрограмма вибросигнала') %Название спектрограммы
%colormap gray %При желании ч/б вариант
end
vev=input('Строить вейвлет-преобразование? 1-да; 0-нет :');
if vev>0;   
figure
subplot(2,1,1), plot (t(10000:10100),Xa(10000:10100)),grid,
title ('Форма сигнала a(t)')
subplot(2,1,2), c=cwt(Xa(10000:10100),1:64,'sym4','plot');
colorbar
title ('Bейвлет-преобразование a(t)') %Название
colormap(pink(64)) %При желании ч/б вариант
end

YYa=rceps(Xa);%Вещественный кепстр

%БЛОК ВЫДЕЛЕНИЯ ОГИБАЮЩЕЙ 
disp('БЛОК ВЫДЕЛЕНИЯ ОГИБАЮЩЕЙ')
W1og=11000; W2og=13200;
pfog=input('Изменить полосу фильтрации (по умолчанию 11 - 13 кГц)? 1-да; 0-нет :');
if pfog>0;
W1og=input('Введите нижную частоту пропусканиЯ W1=');
W2og=input('Введите верхнию частоту пропусканиЯ W2=');
end
%Проектируем полосовой фильтр 
Wp=[W1og*2/Fs W2og*2/Fs];
Ws=[(W1og-0.1*W1og)*2/Fs (W2og+0.1*W2og)*2/Fs]; Rp=0.2;Rs=10;
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
Xog=filter(b,a,Xa); %Фильтруем исходный сигнал режекторным фильтром
XaHi0=abs(hilbert(Xog)); %ВыделЯем огибающую с помощью преобразованиЯ Гильберта 
%Проектируем ФВЧ
Wp=5*2/Fs;                    %Частота пропусканиЯ 5 Гц
Ws=3*2/Fs; Rs=3;Rp=1;         %Частота задерживаниЯ 3Гц, неравномерность в полосах
[n,Wn]=buttord(Wp,Ws,Rp,Rs);   %Расчет порЯдка фильтра Баттерворда
[b,a]=butter(n,Wn,'high');      %Расчет фильтра
XaHi=filter(b,a,XaHi0); %Фильтруем исходный сигнал ФВЧ
YaHi=fft(XaHi,Ff); %Преобразование Фурье
FyyaHi=YaHi.*conj(YaHi)/Ff;    %Амплитуды спектра мощности
abYHi=abs(YaHi)*2/Ff;
f2=Fs*(0:63)/Ff;   %Частоты спектра

figure,%Создаем отдельное графческое окно
%plot (f, FyyaHi(1:151),'k'),grid
%title('Спектр мощности огибающей')
%xlabel('f, Гц')
%ylabel('А, (м/с2)^2*с/Гц')
plot (f2, abYHi(1:64),'k'),grid
title('Спектр огибающей')
xlabel('f, Гц')
%ylabel('А, м/с^2')
%figure,%Создаем отдельное графческое окно
%plot (t(2:800), YYa(2:800),'k'),grid
%title('Кепстр')
%xlabel('t, сек')
%ylabel('С(t)')

        %БЛОК РАСЧЕТА ПИК-ФАКТОРА СИГНАЛА В ПОЛОСЕ ЧАСТОТ W1pf-W2pf кГц
disp('БЛОК РАСЧЕТА ПИК-ФАКТОРА СИГНАЛА В ПОЛОСЕ ЧАСТОТ')
W1pf=9000;
W2pf=11000;
pf=input('Изменить полосу фильтрации (по умолчанию 9 - 11 кГц)? 1-да; 0-нет :');
if pf>0;
W1pf=input('Введите нижную частоту пропусканиЯ W1pf=');
W2pf=input('Введите верхнию частоту пропусканиЯ W2pf=');
end
%Проектируем полосовой фильтр 
Wp=[W1pf*2/Fs W2pf*2/Fs];
Ws=[(W1pf-0.05*W1pf)*2/Fs (W2pf+0.05*W2pf)*2/Fs]; Rp=0.2;Rs=10;
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
[H,w]=freqz(b,a,1024);
Xpf=filter(b,a,Xa); %Фильтруем исходный сигнал режекторным фильтром
%figure (3)
%plot(w*Fs/(2*pi),abs(H)),grid
%Разбиваем сигнал на отрезкi to 1499
Xa01=buffer(Xpf,1499);
PF1=0;
for I=1:79
Xpf0=Xa01(:,I);
skz1=sqrt(sum(Xpf0.^2)/1499);
PF=max(abs(Xpf0))/skz1;
PF1=PF1+PF;
end
PF=PF1/79;
%РАСЧЕТ ЗНАЧЕНИЙ ДИАГНОСТИЧЕСКИХ ПРИЗНАКОВ
disp('ЗНАЧЕНИЯ ДИАГНОСТИЧЕСКИХ ПРИЗНАКОВ')
%1. CKZ
skz0
%2. ЗначениЯ амплитуд виброскорости, нормиров. виброскор и амплитуд огибающей
disp('АМПЛ.СКОР МОЩН.CКОР  АМПЛ.ОГИБ.')
AvpS=[max(abYv(N0p:N1p)) max(Fyyv(N0p:N1p)) max(FyyaHi(N0p:N1p))]
AvrolS=[max(abYv(N0rol:N1rol)) max(Fyyv(N0rol:N1rol)) max(FyyaHi(N0rol:N1rol))]
AvcepS=[max(abYv(N0cep:N1cep)) max(Fyyv(N0cep:N1cep)) max(FyyaHi(N0cep:N1cep))]
AvvnS=[max(abYv(N0vn:N1vn)) max(Fyyv(N0vn:N1vn)) max(FyyaHi(N0vn:N1vn))]
AvtkS=[max(abYv(N0tk:N1tk)) max(Fyyv(N0tk:N1tk)) max(FyyaHi(N0tk:N1tk))]
AvnarS=[max(abYv(N0nar:N1nar)) max(Fyyv(N0nar:N1nar)) max(FyyaHi(N0nar:N1nar))]
cpAv=(AvrolS+AvcepS+AvvnS+AvtkS+AvnarS)/5
%3. ЗначениЯ амплитуд кепстра
disp('АМПЛИТУДЫ КЕПСТРА')
%Ctp=max(YYa(Nt0p:Nt1p))
%Ctrol=max(YYa(Nt0rol:Nt1rol))
%Ctvcep=max(YYa(Nt0cep:Nt1cep))
%Ctvvn=max(YYa(Nt0vn:Nt1vn))
%Ctvtk=max(YYa(Nt0tk:Nt1tk))
%Ctvnar=max(YYa(Nt0nar:Nt1nar))
%4. Значение ПИК-ФАКТОРА
disp('ПИК-ФАКТОР')
PF
% ПОСТАНОВКА ДИАГНОЗА
if PF<=3.5
    Dtpf='Диагноз ПФ: Хорошее состоЯние подшипника';
elseif 3.5<PF<=4.5
    Dtpf='Диагноз ПФ: Предаварийное состоЯние подшипника';
    elseif PF>4.5
    Dtpf='Диагноз ПФ: Аварийное состоЯние подшипника';
end
Dtpf