**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# MatlabScript for adding Gaussian Noise

4 views (last 30 days)

Show older comments

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

Walter Roberson
on 30 Nov 2022

john amoo otoo
on 30 Nov 2022

Walter Roberson
on 30 Nov 2022

Ah but your code uses level with a lower case L but the mat has Level with an uppercase L

john amoo otoo
on 30 Nov 2022

Let me check

john amoo otoo
on 30 Nov 2022

Walt,

I have updated the script and I now have an error at feat_fault.The error says I need to add an end to a function. I have added the end and I still have an error

%% 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

%% Chop the signal according to a sliding window approach

% allocate memory

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,'db8'); % 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

signals = Data3Phase0hmsFaultData;

for dims =1:Nsignals

x=signals(:,dims);

winsize = 32;

numwin = 16;

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,'db8');

% 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

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);

end

>> Noisysig

Error: File: Noisysig.m Line: 101 Column: 1

All functions in a script must be closed with an 'end'.

>>

Walter Roberson
on 30 Nov 2022

for dims =1:Nsignals

signals = Data3Phase0hmsFaultData;

for dims =1:Nsignals

Notice that the second for loop is the same as the first one. You do not have an end corresponding to the first version of the loop.

john amoo otoo
on 1 Dec 2022

Walt. thanks this script seems to work without any errors but I do not get any output results in the work spave and plots

%% 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;

%% Chop the signal according to a sliding window approach

% allocate memory

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);

signals = Data3Phase0hmsFaultData;

for dims =1:Nsignals

x=signals(:,dims);

end

winsize = 32;

numwin = 16;

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,'db8');

% 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

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

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

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

% Plot & Observe the data

subplot(2,1,1)

plot(noisySig)

title('noisySignal')

end

john amoo otoo
on 1 Dec 2022

Thank you Walt. Code/script runs perfect now

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
on 1 Dec 2022

What difficulty are you encountering with the code you have that you commented out?

### Answers (0)

### See Also

### Categories

### Tags

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)