Determination of Standard Deviation and Means square Error

Question

john amoo otoo 2022-11-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1861903-determination-of-standard-deviation-and-means-square-error

评论： Walter Roberson 2022-12-7

Can any one look in to the error I am getting from extracting the standard deviation and Mean Square Error of Matlab DTW of extracted features for fault signal and Steady state signals

14 个评论
显示 12更早的评论隐藏 12更早的评论

john amoo otoo 2022-11-30

移动：John D'Errico 2022-12-5

% GETMSWTFEAT Gets the Multiscale Wavelet Transform features, these

% include: Energy, Variance, Standard Deviation, and Waveform Length

% feat = getmswtfeat(x,winsize,wininc,SF)

% ------------------------------------------------------------------

% The signals in x are divided into multiple windows of size

% "winsize" and the windows are spaced "wininc" apart.

% Inputs

% ------

% signals: columns of signals

% winsize: window size (length of x)

% wininc: spacing of the windows (winsize)

% SF: sampling frequency

%

% Outputs

% -------

% =========================================================================

function feature_out = getmswtfeat(signals,winsize,wininc,SF)

if nargin < 4

if nargin < 3

if nargin < 2

error('A sliding window approach requires the window size (winsize) as input')

end

error('A sliding window approach requires the window increment (wininc) as input')

end

error('Please provide the sampling frequency of this signal')

end

%% The number of decomposition levels

decomOption = 1;

if decomOption==1

J=8; % Number of decomposition levels set manually here

elseif decomOption==2

J=wmaxlev(winsize,'sym2'); % Number of decomposition levels set based on window size and wavelet family

else

J=(log(SF/2)/log(2))-1; % Number of decomposition levels set based on sampling frequency (SF)

end

%% make sure you have some parameters pre-defined

% specify the number of samples

datasize = size(signals,1);

% based on the number of samples, winsize, and wininc, how many windows we

% will have? this is "numwin"

numwin = floor((datasize - winsize)/wininc)+1;

% how many signals (electrodes) are we processing

Nsignals = size(signals,2);

% how many features we plan to extract

%%

NF = 8;

% predefine zeros matrix to allocate memory for output features

%feature_out = zeros(numwin,(J+1)*NF*Nsignals);

feature_out = zeros(numwin,(J+1)*Nsignals);

for dims =1:Nsignals

x=signals(:,dims);

%% Chop the signal according to a sliding window approach

% allocate memory

feat = zeros(winsize,numwin);

st = 1;

en = winsize;

for i = 1:numwin

feat(1:winsize,i) = x(st:en,:)-mean(x(st:en,:));

st = st + wininc;

en = en + wininc;

end

%% Multisignal one-dimensional wavelet transform decomposition

dec = mdwtdec('col',feat,J,'sym2');

% Proceed with Multisignal 1-D decomposition energy distribution

if isequal(dec.dirDec,'c')

dim = 1;

end

[cfs,longs] = wdec2cl(dec,'all');

level = length(longs)-2;

if dim==1

cfs = cfs';

longs = longs';

end

numOfSIGs = size(cfs,1);

num_CFS_TOT = size(cfs,2);

absCFS = abs(cfs);

absCFS0 = (cfs);

cfs_POW2 = absCFS.^2;

Energy = sum(cfs_POW2,2);

percentENER = zeros(size(cfs_POW2));

notZER = (Energy>0);

percentENER(notZER,:) = cfs_POW2(notZER,:);%./Energy(notZER,ones(1,num_CFS_TOT));

%% or try this version below and tell us which one is the best on your data

% percentENER(notZER,:) = cfs_POW2(notZER,:);

%% Pre-define and allocate memory

tab_ENER = zeros(numOfSIGs,level+1);

tab_VAR = zeros(numOfSIGs,level+1);

tab_STD = zeros(numOfSIGs,level+1);

tab_WL = zeros(numOfSIGs,level+1);

tab_entropy = zeros(numOfSIGs,level+1);

%% Feature extraction section

st = 1;

for k=1:level+1

nbCFS = longs(k);

en = st+nbCFS-1;

%tab_ENER(:,k) = sum(percentENER(:,st:en),2);% energy per waveform

%tab_VAR(:,k) = var(percentENER(:,st:en),0,2); % variance of coefficients

%tab_STD(:,k) = std(percentENER(:,st:en),[],2); % standard deviation of coefficients

%tab_WL(:,k) = sum(abs(diff(percentENER(:,st:en)')))'; % waveform length

%percentENER(:,st:en) = percentENER(:,st:en)./repmat(sum(percentENER(:,st:en),2),1,size(percentENER(:,st:en),2));

prob = percentENER(:,st:en);%./repmat(sum(percentENER(:,st:en),2),1,longs(k)) + eps;

tab_entropy(:,k) = -sum(prob.*log(prob),2);%./size(percentENER(:,st:en),2);

st = en + 1;

end

%feature_out(:,(1:(NF*(J+1)))+(dims-1)*((J+1)*NF)) =log([tab_ENER tab_VAR tab_STD tab_WL tab_entropy]);

feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_entropy;

feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_STD;

feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_VAR;

feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_WL;

end

john amoo otoo 2022-11-30

移动：John D'Errico 2022-12-5

load('John_FaultData0km0ohm.mat')

%% does this signal has any NaNs, if so remove

SteadyStateNoneFaultState = rmmissing(SteadyStateNoneFaultState);

DataPG0km0hmsFaultData = rmmissing(DataPG0km0hmsFaultData);

SSTime = SteadyStateNoneFaultState.Time;

SSNFS = SteadyStateNoneFaultState.SteadyStateNoneFaultState;

DPTime = DataPG0km0hmsFaultData.Time;

DPFD = DataPG0km0hmsFaultData.DataPG0km0hmsFaultData;

%% Define the filtering and inspect

n = 8; %% defines window length

w = [-ones(n,1); ones(n,1)];

SSNFS = filter(w, n, SSNFS);

DPFD = filter(w, n, DPFD);

%% Normalize signals

med_training = prctile(SSNFS,50);

iqr_training = iqr(SSNFS);

SSNFS = (SSNFS-med_training)./iqr_training;

med_fault = prctile(DPFD,50);

iqr_fault = iqr(DPFD);

DPFD = (DPFD-med_fault)./iqr_fault;

% Plot & Observe the data

subplot(2,1,1)

plot(DPTime, DPFD)

title('filtered DataPG0km5hmsFaultData')

subplot(2,1,2)

plot(SSTime, SSNFS)

title('filtered SteadyStateNoneFaultState')

%% Let's observe the FFT power spectrum for differences

feat_fault = getmswtfeat(DPFD,32,16,100000);

feat_Good = getmswtfeat(SSNFS,32,16,100000);

john amoo otoo 2022-12-6

编辑：Walter Roberson 2022-12-6

在 MATLAB Online 中打开

%% Feature extraction section
st = 1;
for k=1:level+1
    nbCFS                = longs(k);
    en                   = st+nbCFS-1;
    %tab_ENER(:,k)        = sum(percentENER(:,st:en),2);% energy per waveform
    %tab_VAR(:,k)         = var(percentENER(:,st:en),0,2); % variance of coefficients
    %tab_STD(:,k)         = std(percentENER(:,st:en),[],2); % standard deviation of coefficients
    %tab_WL(:,k)          = sum(abs(diff(percentENER(:,st:en)')))'; % waveform length
    %percentENER(:,st:en) = percentENER(:,st:en)./repmat(sum(percentENER(:,st:en),2),1,size(percentENER(:,st:en),2));
    prob                 = percentENER(:,st:en);%./repmat(sum(percentENER(:,st:en),2),1,longs(k)) + eps;
    tab_entropy(:,k)     = -sum(prob.*log(prob),2);%./size(percentENER(:,st:en),2);
    tab_VAR(:,k)         = var(percentENER(:,st:en),0,2);% variance of coefficients
    tab_STD(:,k)         = std(percentENER(:,st:en),[],2);% standard deviation of coefficients
    tab_WL(:,k)          = sum(abs(diff(percentENER(:,st:en)')))'; % waveform length
    st = en + 1;
end
%feature_out(:,(1:(NF*(J+1)))+(dims-1)*((J+1)*NF)) =log([tab_ENER tab_VAR tab_STD tab_WL tab_entropy]);
feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_entropy;
feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_STD;
feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_VAR;
feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_WL;
end

john amoo otoo 2022-12-7

移动：Walter Roberson 2022-12-7

在 MATLAB Online 中打开

% GETMSWTFEAT Gets the Multiscale Wavelet Transform features, these
% include: Energy, Variance, Standard Deviation, and Waveform Length
% feat = getmswtfeat(x,winsize,wininc,SF)
% ------------------------------------------------------------------
% The signals in x are divided into multiple windows of size
% "winsize" and the windows are spaced "wininc" apart.
% Inputs
% ------
%    signals:	columns of signals
%    winsize:	window size (length of x)
%    wininc:	spacing of the windows (winsize)
%    SF:        sampling frequency
%
% Outputs
% -------
% =========================================================================
% REFERENCE: MATLAB CODE: Multi Scale Wavelet Decomposition: Dr. Rami Khushaba
% Email: Rami.Khushaba@sydney.edu.au
% URL: www.rami-khushaba.com (Matlab Code Section)
function feature_out = getmswtfeat(signals,winsize,wininc,SF)
if nargin < 4
    if nargin < 3
        if nargin < 2
            error('A sliding window approach requires the window size (winsize) as input')
        end
        error('A sliding window approach requires the window increment (wininc) as input')
    end
    error('Please provide the sampling frequency of this signal')
end
%% The number of decomposition levels 
decomOption = 1;
if decomOption==1
    J=8; % Number of decomposition levels set manually here
elseif decomOption==2
    J=wmaxlev(winsize,'db2'); % Number of decomposition levels set based on window size and wavelet family
else
    J=(log(SF/2)/log(2))-1; % Number of decomposition levels set based on sampling frequency (SF)
end
%% make sure you have some parameters pre-defined
% specify the number of samples
datasize = size(signals,1);
% based on the number of samples, winsize, and wininc, how many windows we
% will have? this is "numwin"
numwin = floor((datasize - winsize)/wininc)+1;
% how many signals (electrodes) are we processing
Nsignals = size(signals,2);
% how many features we plan to extract
%% 
NF = 8;
% predefine zeros matrix to allocate memory for output features
%feature_out = zeros(numwin,(J+1)*NF*Nsignals);
feature_out = zeros(numwin,(J+1)*Nsignals);
for dims =1:Nsignals
    x=signals(:,dims);
    %% Chop the signal according to a sliding window approach
    % allocate memory
    feat = zeros(winsize,numwin);
    st = 1;
    en = winsize;
    for i = 1:numwin
        feat(1:winsize,i) = x(st:en,:)-mean(x(st:en,:));
        st = st + wininc;
        en = en + wininc;
    end
    
    %% Multisignal one-dimensional wavelet transform decomposition
    dec = mdwtdec('col',feat,J,'sym8');
    % Proceed with Multisignal 1-D decomposition energy distribution
    
    if isequal(dec.dirDec,'c')
        dim = 1;
    end
    [cfs,longs] = wdec2cl(dec,'all');
    level = length(longs)-2;
    
    if dim==1
        cfs = cfs';
        longs = longs';
    end
    numOfSIGs             = size(cfs,1);
    num_CFS_TOT           = size(cfs,2);
    absCFS                = abs(cfs);
    absCFS0               = (cfs);
    cfs_POW2              = absCFS.^2;
    Energy                = sum(cfs_POW2,2);
    percentENER           = zeros(size(cfs_POW2));
    notZER                = (Energy>0);
    percentENER(notZER,:) = cfs_POW2(notZER,:);%./Energy(notZER,ones(1,num_CFS_TOT));
    
    %% or try this version below and tell us which  one is the best on your data
    % percentENER(notZER,:) = cfs_POW2(notZER,:);
    
    %% Pre-define and allocate memory
    tab_ENER    = zeros(numOfSIGs,level+1);
    tab_VAR     = zeros(numOfSIGs,level+1);
    tab_STD     = zeros(numOfSIGs,level+1);
    tab_WL      = zeros(numOfSIGs,level+1);
    tab_entropy = zeros(numOfSIGs,level+1);
    
    %% Feature extraction section
    st = 1;
    for k=1:level+1
        nbCFS                = longs(k);
        en                   = st+nbCFS-1;
        %tab_ENER(:,k)        = sum(percentENER(:,st:en),2);% energy per waveform
        %tab_VAR(:,k)         = var(percentENER(:,st:en),0,2); % variance of coefficients
        %tab_STD(:,k)         = std(percentENER(:,st:en),[],2); % standard deviation of coefficients
        %tab_WL(:,k)          = sum(abs(diff(percentENER(:,st:en)')))'; % waveform length
        %percentENER(:,st:en) = percentENER(:,st:en)./repmat(sum(percentENER(:,st:en),2),1,size(percentENER(:,st:en),2));
        prob                 = percentENER(:,st:en);%./repmat(sum(percentENER(:,st:en),2),1,longs(k)) + eps;
        tab_entropy(:,k)     = -sum(prob.*log(prob),2);%./size(percentENER(:,st:en),2);
        tab_VAR(:,k)         = var(percentENER(:,st:en),0,2);% variance of coefficients
        tab_STD(:,k)         = std(percentENER(:,st:en),[],2);% standard deviation of coefficients
        tab_WL(:,k)          = sum(abs(diff(percentENER(:,st:en)')))'; % waveform length
        st = en + 1;
    end
    %feature_out(:,(1:(NF*(J+1)))+(dims-1)*((J+1)*NF)) =log([tab_ENER tab_VAR tab_STD tab_WL tab_entropy]);
    feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_entropy;
    feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_STD;
    feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_VAR;
    feature_out(:,(1:((J+1)))+(dims-1)*(J+1)) =tab_WL;
end