%Abd Kadir Bin Mahamad.
%This code takes data received from a bearing and converts it from the time
%domain to the frequency domain in order to determine bearing defect size.
%This code also includes the HFRT to determine where the bearing fails.
% Figures generated
% Figure 1 = Raw Data
% Figure 2 = Raw data spectrum (fft)
% Figure 3 = Envelope spectrum
clear;
%close all;
%constants for the system for paris's equation
Co=0.13;
n=1.1;
%The following variables dimention the bearing used in research 
rpm=1750; %[rev/minute] bearing rotation speed
f=rpm/60; %[rev/minute] bearing rotation speed
ID=0.9843; %[inches] inside diameter
OD=2.0472; %[inches] outside diameter
no_balls=19; %[balls] number of balls in the bearing
ball_d=0.3126; %[inches] ball diameter
alphadeg=13.13; %[degrees]
alpha=alphadeg*(pi/180);%[radians] taper angle
dm=1.537; %[inches] pitch diameter
%determine frequencies of ball pass inner race, ball pass outer race,
%and rollers. Defect is on outer race for this research.
f_bpir = no_balls*f/2*(1+(ball_d/dm)*cos(alpha));
f_bpor = no_balls*f/2*(1-(ball_d/dm)*cos(alpha));
f_r = no_balls*(f*dm)/2*(1-((ball_d/dm)*cos(alpha))^2);
f_c =f/2*(1-(ball_d/dm)*cos(alpha));
%Load data and find pertinent parameters
[filename, pathname,X] = ...
     uigetfile({'*.mat';'*.*'},'File Selector');

filename = fullfile(pathname, filename);
data = load (filename, '-mat');

%rawdata = load ( file );
%Data1 = load([path02 file02 '.mat']);
%nodes_number = nodes(length(nodes));
%load '12feabnormal.mat';
datas = struct2cell(data);
[p q]=size(datas)
if p>=5
[Y1,Y2,Y3,Y4,Y5] = deal(datas{:});
elseif p>=4
    [Y1,Y2,Y3,Y4] = deal(datas{:});
elseif p>=3
    [Y1,Y2,Y3] = deal(datas{:});
else
    [Y1,Y2] = deal(datas{:});
end
vib_datapoints=length(Y1); 
%number of vibration data points
vib=Y1(1:vib_datapoints,1); %pulls out the correct column of data
vib_time=1.0; %[seconds] length of time data was sampled
vib_fs=30000; %[samples/second]sampling frequency
% This is a plot of the raw data

figure(1);
plot(0:1/vib_datapoints:(1-1/vib_datapoints),vib);
meanvib=mean(vib);
ms=std(vib);
mRMSp=rms(vib);
maxvib=max(vib);
minvib=min(vib);
mskew=skewness(vib);
mkurtos=kurtosis(vib);
mcrestfactor=maxvib/mRMSp;
mCLF=maxvib/(sum(sqrt(abs(vib)))/vib_datapoints).^2;
mSF=mRMSp/(sum(abs(vib))/vib_datapoints);
mIF=maxvib/(sum(abs(vib))/vib_datapoints);

% mR=corrcoef(vib);
% 
% mcov=cov(vib);
% median=median(vib);
% mvar=var(vib);

title(['Raw data'])
xlabel('Time (s)');
ylabel('Acceleration (g)');
% Next step is to find fft/power spectrum
X=fft(vib)/vib_datapoints;
df=vib_fs/vib_datapoints;
%This is a plot of the Power Spectrum (fft(rawdata))

figure(2)
freqvalue=df:df:vib_fs/2-df;
%realplot=freqvalue,abs(X(2:vib_datapoints/2));
plot(freqvalue,abs(X(2:vib_datapoints/2)));
maxps=max(abs(X(2:vib_datapoints/2)));
mfreqcentre=sum(freqvalue'.*abs(X(2:vib_datapoints/2))/sum(abs(X(2:vib_datapoints/2))));
mRMSfreq=sqrt(sum(freqvalue'.^2.*abs(X(2:vib_datapoints/2))/sum(abs(X(2:vib_datapoints/2)))));
bezafreq=freqvalue'-mfreqcentre';
mstdfreq=sqrt(sum(bezafreq.^2.*abs(X(2:vib_datapoints/2))/sum(abs(X(2:vib_datapoints/2)))));
title(['Power Spectrum'])
xlabel('Frequency (Hz)');
ylabel('Power');
%Third step of program is to find the envelope spectrum of the data
%This shows the defect location.
%Resonant frequency of the system is around 5500 Hz
Wn_vib=[4000 6000]*2/vib_fs;
cutoff=400;
cutoff_vib=cutoff*2/vib_fs; %cutoff frequency for bandpassed filtration of data
[top,bottom]=butter(5,Wn_vib); %use a butterworth filter for bandpassing
[t,b]=butter(5,cutoff_vib);
xx=filter(top,bottom,vib);
X=fft(vib)/vib_datapoints;
xxx = abs(xx);
xxx=filter(t,b,xxx);
XXX=fft(xxx)/vib_datapoints; %envelope spectrum output

figure(3)
plot(df:df:cutoff-df,abs(XXX(2:floor(cutoff/df))));
maxen=max(abs(XXX(2:floor(cutoff/df))));
title('Envelope Spectrum');
xlabel('Frequency (Hz)');
ylabel('Power');
% Cutting down Power Spectrum plot to encompass only the energy
%around the resonant frequency of the system.
low_lim=10;
up_lim=4000;
x_vect=df:df:vib_fs/2-df;
y_vect=abs(X(2:vib_datapoints/2));
lxv=length(x_vect);
index_low=round(low_lim/df);
index_high=round(up_lim/df);
diff=index_high-index_low;
for ii=1:1:diff;
x_res(ii)=x_vect(index_low+ii-1);
y_res(ii)=y_vect(index_low+ii-1);
end
format long g;
%maxenergy =max(x_res);
%uisave

% [c1A,c1D]=dwt(vib,'db4');
% [c2A,c2D]=dwt(c1A,'db4');
% [c3A,c3D]=dwt(c2A,'db4');
% [c4A,c4D]=dwt(c3A,'db4');
% mcoed1p=mean(c1D); % calculating mean value of coefficient d1 for primary
% mcoed2p=mean(c2D); % for d2
% mcoed3p=mean(c3D);
% mcoed4p=mean(c4D);
% mcoea4pa=mean(c4A); % for a4

save resultnor0Y3 m*