[WIP] Windowed transformation and Bearing Fault example

rul_pm/dataset/MFPT.py

@@ -0,0 +1,90 @@
+import io
+import numpy as np
+import pandas as pd
+from rul_pm import DATA_PATH
+from rul_pm.dataset.lives_dataset import AbstractLivesDataset
+from scipy.io import loadmat
+from pathlib import Path
+
+DATA_PATH = Path("/home/luciano/fuentes/lru_gcd_bearing/rul_pm/dataset/data/")
+MFPT_PATH = DATA_PATH / "MFPT_Fault_Data_Sets"
+
+
+
+ROLLER_DIAMETER = 0.235
+PITCH_DIAMETER = 1.245
+NUMBER_OF_ELEMENTS = 8
+CONTACT_ANGLE = 0
+SHAFT_RATE = 25
+
+def BPFO():
+    return 0.5*(NUMBER_OF_ELEMENTS * SHAFT_RATE)*(1- (ROLLER_DIAMETER/PITCH_DIAMETER * np.cos(CONTACT_ANGLE)))
+
+def BPFI():
+    return 0.5*(NUMBER_OF_ELEMENTS * SHAFT_RATE)*(1 + (ROLLER_DIAMETER/PITCH_DIAMETER * np.cos(CONTACT_ANGLE)))
+
+def FTF():
+    return 0.5*(SHAFT_RATE)*(1 - (ROLLER_DIAMETER/PITCH_DIAMETER * np.cos(CONTACT_ANGLE)))
+
+def BSF():
+    return 0.5*((PITCH_DIAMETER*SHAFT_RATE)/ROLLER_DIAMETER)*(1 - (ROLLER_DIAMETER/PITCH_DIAMETER * np.cos(CONTACT_ANGLE))**2)
+
+def data_load(filename, label):
+    """
+    This function is mainly used to generate test data and training data.
+    filename:Data location
+    """
+    data = loadmat(filename)["bearing"]
+    data = {name: np.squeeze(data[0][0][i]) for i, (name, _) in enumerate(data.dtype.descr)}
+    df =  pd.DataFrame(data)
+
+    df['label'] = int(label)
+    df['time'] = np.arange(0, df.shape[0]) / df.sr
+    df['filename'] = filename.stem
+
+    return df
+
+CLASS_NORMAL = 0
+CLASS_OUTER_RACE_FAULT = 1
+CLASS_INNER_RACE_FAULT = 2
+
+folders = [
+    ("1 - Three Baseline Conditions", CLASS_NORMAL),
+    ("2 - Three Outer Race Fault Conditions", CLASS_OUTER_RACE_FAULT),
+    ("3 - Seven More Outer Race Fault Conditions", CLASS_OUTER_RACE_FAULT),
+    ("4 - Seven Inner Race Fault Conditions", CLASS_INNER_RACE_FAULT)
+]
+
+
+class MFPTDataset(AbstractLivesDataset):
+    def __init__(self):
+        self.lives = []
+        mat_files = list((MFPT_PATH / folders[0][0]).glob("*.mat"))
+        mat_files = sorted(mat_files, key=lambda x: x.stem)
+        for file in mat_files:
+            self.lives.append(data_load(file, label=CLASS_NORMAL))
+
+        for folder, label in folders[1:]:
+            mat_files = list((MFPT_PATH / folder).glob("*.mat"))
+            mat_files = sorted(mat_files, key=lambda x: x.stem)
+            for file in mat_files:
+                self.lives.append(data_load(file, label=label))
+
+
+    def get_life(self, i):
+        """
+
+        Returns
+        -------
+        pd.DataFrame
+            DataFrame with the data of the life i
+        """
+        return self.lives[i]
+
+    @property
+    def nlives(self):
+        return len(self.lives)
+
+
+
+

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/GetBearFreqRatio.m

@@ -0,0 +1,37 @@
+function d = GetBearFreqRatio(rd,pd,ca,ne,type,side)
+% Cage passing frequency:	 f/2 (1-d/e cos (?) )
+% Inner race passing frequency	bf/2 (1 + d/e cos (?) )
+% Outer race passing frequency	bf/2 (1 - d/e cos (?) )
+% Ball passing frequency	ef/d (1-(d/e)2 cos2(?) )
+% Where 
+% f is the driving frequency 
+% b is the number of rolling elements
+% d is the ball bearing diameter 
+% e is the bearing pitch diameter and
+% ? is the bearing contact angle.
+
+%Get the bearing Freq. ratio
+%rd:    roller diameter
+%pd:    pitch diameter
+%ca:    contact angle in degrees
+%ne:    number of elements
+%type   enum[cage, roller, outer race, inner race]
+%side: inner or outer race fixed.  inner = 1, outer = 2;
+%Eric Bechhoefer, April 10, 2009 for PHM Conference
+%Ref: Timken Bearings Catalog 2008
+rdpd = rd/pd;
+cs = cos(ca*pi/180);
+pdrd = 1/rdpd;
+if type == 1,       %cage freq ratio
+    if side == 1,
+        d = 0.5*(1-rdpd*cs);
+    else
+        d = 0.5*(1+rdpd*cs);
+    end
+elseif type == 2,   %roller freq ratio
+    d = (1-(rdpd*cs)^2)*pdrd;
+elseif type == 3,   %outer race freq ratio
+    d = 0.5*(1-rdpd*cs)*ne;
+else                %inner race freq 
+    d = 0.5*(1+rdpd*cs)*ne;
+end

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/NiceBearing.m

@@ -0,0 +1,16 @@
+function [rd,pd,ca,ne,side] = NiceBearing
+
+%rd:    roller diameter
+%pd:    pitch diameter
+%ca:    contact angle in degrees
+%ne:    number of elements
+%type   enum[cage, roller, outer race, inner race]
+%side: inner or outer race fixed.  inner = 1, outer = 2;
+%Eric Bechhoefer, April 10, 2009 for PHM Conference
+rd = .235;
+ca = 0;
+ne = 8;
+pd = 1.245;
+
+side = 2;
+

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/bearingAnalysis.m

@@ -0,0 +1,52 @@
+function cis = bearingAnalysis(v,zct)
+
+cis = zeros(1,6);
+
+sr = 97656;
+faultRates = [0.42 2.87 9.25 6.72 1];%cage, ball, inner, outer, 1/rev
+idx = find(zct>6,1); %only interested in first six seconds of data.
+rate = mean(1./diff(zct(1:idx))/2); 
+cis(6) = rate;
+faultFreq = faultRates * rate;
+
+ [env,dty] = envelope1(v,1/sr,128,8500,11500);
+
+ [spec, freq] = psde(env, 4096,1/dty, 2048);
+
+ plot(freq,spec);
+ ax = axis();
+ hold on
+ for i = 1:5,
+     f = faultFreq(i);
+     plot([f f],ax(3:4),'LineWidth',2)
+     cis(i) = BrngEng(spec,freq(2),f);
+ end
+ hold off
+ legend('env','cage','ball','outer','inner','1/rev')
+ xlabel('Hz')
+ ylabel('Gs')
+ axis([0 200 ax(3:4)])
+
+ function e = BrngEng(P,dFrq,hz)
+
+hLw = hz*.97;
+hHi = hz*1.02;
+
+bLw = floor(hLw/dFrq);
+if bLw == 0, 
+    bLw = 1;
+end
+bHi = ceil(hHi/dFrq);
+rng = bLw:bHi;
+[e,idx] = max(P(rng));
+
+% n = length(P);
+% freq = (0:(n-1))*dFrq;
+% plot(freq,P,(rng(idx)-1)*dFrq,e,'x','LineWidth',2)
+% xlabel('Frequency (Hz)')
+% ylabel('Gs (Hz)')
+% ax = axis;
+% h = line([hz hz],ax(3:4)); 
+% set(h,'LineWidth',2);
+% set(h,'Color',[1 0 0]);
+% axis([0 500 0 ax(4)])

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/bearingScript.m

@@ -0,0 +1,31 @@
+function bearingScript(bearing)
+%runs a the beairng analysis. This is from the MFPT data set
+if nargin == 0
+    load InnerRaceFault_vload_7.mat
+end
+
+
+[rd,pd,ca,ne,side] = NiceBearing; %get the dimensions of the nice bearing
+faultRates = ones(5,1); %[cage, ball, outer, inner, shaft];
+for i = 1:4
+    faultRates(i) = GetBearFreqRatio(rd,pd,ca,ne,i,side);%calculate the bearing fault rate
+end
+faultFreq = faultRates * bearing.rate;
+
+% [env,dty] = envelope1(bearing.gs,1/bearing.sr,128,2000,4000);
+[env,dty] = envelope2(bearing.gs,1/bearing.sr,2000,4000);
+
+[spec, freq] = psde(env, 8192,1/dty, 4096);
+
+plot(freq,spec);
+ax = axis();
+hold on
+for i = 1:5
+    f = faultFreq(i);
+    plot([f f],ax(3:4),'LineWidth',2)
+end
+hold off
+legend('env','cage','ball','outer','inner','1/rev')
+xlabel('Hz')
+ylabel('Gs')
+axis([0 200 ax(3:4)])

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/envelope1.m

@@ -0,0 +1,28 @@
+function [env,dty] = envelope1(data,dt,nfilt,lowf,highf)
+% [env,dty] = envelope1(data,dt,nfilt,lowf,highf);
+%Inputs:
+% data      :data vector, time domain
+% dt        :sampling time interval
+% nfilt     :size of filter for low pass filtering
+% lowf      :low frequency limit of bandpass filter
+% highf     :high frequency limit of bandpass filter
+%Outputs:
+% env       :Envelope of data
+% dty       :decimated sample rate
+
+c = (highf + lowf)/2;
+if 1/dt/2 < c
+    y = [];
+    dt = [];
+else
+
+    ndata = length(data);
+    z = data(:) .*exp(-2 * pi * 1i * c * dt * (0:ndata-1)');
+    bw = highf - lowf;
+    b = fir1(nfilt,bw*dt);
+    x = filter(b,1,z);
+    r = fix(1/(bw*2*dt));
+    env = abs( x(nfilt+1:r:ndata) );      % crop first nfilt elements
+
+    dty = dt*r;
+end

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/envelope2.m

@@ -0,0 +1,27 @@
+function [env,dty] = envelope2(data,dt,lowf,highf)
+% [env,dty] = envelope1(data,dt,nfilt,lowf,highf);
+%Inputs:
+% data      :data vector, time domain
+% dt        :sampling time interval
+% lowf      :low frequency limit of bandpass filter
+% highf     :high frequency limit of bandpass filter
+%Outputs:
+% env       :Envelope of data
+% dty       :decimated sample rate
+
+n = length(data);
+dfq = 1/dt/n;
+idxLow = floor(lowf/dfq);
+idxHi = ceil(highf/dfq);
+D = fft(data);
+idx = idxHi-idxLow + 1;
+
+D(1:idx) = D(idxLow:idxHi);
+D(idx+1:end) = 0;
+data = abs(ifft(D));
+bw = highf - lowf;
+r = fix(1/(bw*2*dt));
+env =  data(1:r:n);   % crop first nfilt elements
+
+dty = dt*r;
+

rul_pm/dataset/data/MFPT_Fault_Data_Sets/5 - Analyses/psde.m

@@ -0,0 +1,56 @@
+function [Spec, freq] = psde(x, winln,Fs, noverlap)
+%[Spec, freq] = psde(x, winln,Fs, noverlap);
+%Energy power spectrual density
+%using Welch's method using Hann window
+% x:        time domain data
+% windln:   window length, e.g. 2048,
+% Fs;       sampleing frequency
+% noverlap: length of the overlp
+n = winln;
+m = n/2;
+window = .5*(1 - cos(2*pi*(1:m)'/(n+1))); %Hann window
+window = [window; window(end:-1:1)];
+window = window(:);
+
+nfft = length(window);
+
+n = length(x);		    % Number of data points
+nwind = length(window); % length of window
+if n < nwind            % zero-pad x if it has length less than the window length
+    x(n+1:nwind)=0;  
+    n=nwind;
+end
+% Make sure x is a column vector; do this AFTER the zero-padding 
+% in case x is a scalar.
+x = x(:);		
+
+k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
+index = 1:nwind;
+
+KMU = k*sum(window)^2;% alt. Normlzng scale factor ==> peaks are about right
+
+Spec = zeros(nfft,1); 
+for i=1:k
+    xw = window.*detrend(x(index));
+    index = index + (nwind - noverlap);
+    xx = fft(xw,nfft);
+    Xx = abs(xx).^2;
+    Spec = Spec + Xx;
+end
+
+% Select first half
+if ~any(any(imag(x)~=0)),   % if x is not complex
+    if rem(nfft,2),    % nfft odd
+        select = (1:(nfft+1)/2)';
+    else
+        select = (1:nfft/2)';
+    end
+    Spec = Spec(select);
+else
+    select = (1:nfft)';
+end
+freq = (select - 1)*Fs/nfft;
+
+
+Spec = sqrt(Spec*(4/KMU));   % normalize
+

rul_pm/graphics/kurtogram.py

@@ -0,0 +1,36 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def plot_kurtogram(Kwav, freq, Level_w, ax=None):
+    """
+    Plots the kurtogram.
+
+    Parameters
+    -----------
+    Kwav: np.array
+        kurtogram
+    freq_w: np.array
+        frequency vector
+    Level_w: np.array
+        vector of levels
+    ax: Optional[Axis]
+
+    """
+    if ax is None:
+        _, ax = plt.subplots(figsize=(9, 5))
+
+    aspect = (1/(Kwav.shape[0] / Kwav.shape[1]))
+
+    im = ax.imshow(np.clip(Kwav, 0, np.inf), aspect=aspect, cmap="viridis")
+    x_ticks = [x for x in ax.get_xticks() if x >= 0 and x < len(freq)]
+    ax.set_xticks(x_ticks)
+
+    ax.set_xticklabels([np.round(freq[int(i)] / 1000, 2) for i in x_ticks])
+    ax.set_xlabel("Frequency [kHz]")
+    ax.set_ylabel("Level")
+    ax.set_yticks(list(range(len(Level_w))))
+    ax.set_yticklabels([np.round(j, 2) for j in Level_w])
+    ax.grid(False)
+    ax.figure.colorbar(im)
+    return ax