MatlabScript for adding Gaussian Noise

4 views (last 30 days)
john amoo otoo
john amoo otoo on 30 Nov 2022
Commented: Walter Roberson on 1 Dec 2022
Please I have been able to fix my script but ran into an error. Attached is the script for adding a Gaussian Noise to a signal and extracting the features
%% clear
load('John_3PhaseFault0ohm.mat')
%% does this signals has any NaNs, if so remove
SteadyStateNoneFaultState = rmmissing(SteadyStateNoneFaultState);
Data3Phase0hmsFaultData = rmmissing(Data3Phase0hmsFaultData);
SSTime = SteadyStateNoneFaultState.Time;
SSNFS = SteadyStateNoneFaultState.SteadyStateNoneFaultState;
DPTime = Data3Phase0hmsFaultData.Time;
DPFD = Data3Phase0hmsFaultData.Data3Phase0hmsFaultData;
%% Feature extraction section
for k=1:level+1
signals = Data3Phase0hmsFaultData;
reqSNR = [15]; %noise in dB
sigEner = norm(signals(:,k))^2; % energy of the signals
noiseEner = sigEner/(10^(reqSNR/10)); % energy of noise to be added
noiseVar = noiseEner/(length(signals)); % variance of noise to be added
ST = sqrt(noiseVar); % std. deviation of noise to be added
noise = ST*randn(size(signals)); % noise
noisySig = signals+noise; % noisy signals
end
% Plot & Observe the data
subplot(2,1,1)
plot(noisySig);
title('Noise Signal')
subplot(2,1,2)
plot(noise);
title('Noise')
%% Let's observe the FFT power spectrum for differences
feat_fault = getmswtfeat(noisySig,32,16,100000);
Error
>> Noisysig
Unrecognized function or variable 'level'.
Error in Noisysig (line 11)
for k=1:level+1
>>
  10 Comments
john amoo otoo
john amoo otoo on 1 Dec 2022
Edited: Walter Roberson on 1 Dec 2022
Walt
Walter Roberson,please could you look into this. I want to add a scipt to tabulate standard deviation, mean square error and variance.Anyhelp on this
% 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,'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
Walter Roberson
Walter Roberson on 1 Dec 2022
What difficulty are you encountering with the code you have that you commented out?

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!