Matlab Script that I can add a Gaussian Noise to a Discrete Wavelet Transform processed fault signal

Please does any one has a matlab script where I can add a Gausssian Noise to a Discrete Wavelete Transform processed fault signal. This will mimic a lightning disturbance and switching disturbance during fault connditions and normal steady state. Forf the protective relay to distinguish between this lightning disturbance and fault conditions

Answer 1

[William Rose]

@John Amoo-Otoo,

If you have a fault signal x(t), processed by a DWT (or not), add Gaussian noise as shown.

s=1;                    %standard deviation of the added noise
xn=x+s*randn(size(x));  %x=signal without added noise

The fact that the signal has been processed by a DWT is not really relevant, unless there is more information that you want to share. Good luck with your work on lightning and faults.

[John Amoo-Otoo]

Thank you William

[William Rose]

@John Amoo-Otoo, you're welcome.

[John Amoo-Otoo]

William, quick question, if the Gaussian Noise is 10db, 20db, 30db, 40 db. How does it fit in this equation for the script?

s=1; %standard deviation of the added noise

xn=x+s*randn(size(x)); %x=signal without added noise

[William Rose]

@John Amoo-Otoo,

[edit: fix typos; change "y" to "x", to be consistent with my earlier comment and code]

To get noise with different amplitudes, you adjust s, the standard deviation.

In general, s=10^(a/20), where a is the noise amplitude in dB. Therefore

a=[10,20,30,40];    %noise in dB
s=10.^(a/20)        %noise amplitude factor

s = 1×4

3.1623 10.0000 31.6228 100.0000

dB is usally defined in relative terms. For example, the signal to noise ratio (SNR) is usually expressed in dB. If the dB values which you quote are intended as SNRs, then the values of s computed above should be divided by the signal amplitude. If you do that, then the SNR will be 10, 20, 30, or 40 dB. The only problem with this advice is that the amplitude of the signal may not be totally obvious. If the signal is a sinusoid, x=A*cos(wt), then the amplitude is A. If the signal is a complicated real world signal, then what is the noise-free amplitude? At this point some judgement will be required. You could go with A=std(x), or A=(max(x)-min(x))/2, or other ideas, where x is the noise-free signal.

[John Amoo-Otoo]

Thank you William

[John Amoo-Otoo]

William,

Please attached is the error that I got when I added the small script of adding Gaussian Noise to the original signal

[William Rose]

@John Amoo-Otoo,

Line 115 of function getmswfeat.m is

prob = percentENER(:,st:en);

This line should assign a subset of array percentENER() to the array prob.

The error indicates that the range st:en extends beyond the number of columns in array percentENER. The error indicates that percentENER has only 52 columns. One can infer from this error rmessage that st:en includes one or more values greater than 52. Therefore it is impossible to access columns st:en, which is an error.

[John Amoo-Otoo]

@william, so please how do I correct it to havemore than 52 columns

[William Rose]

@John Amoo-Otoo, to answer that I would need to run your script. processdata.m is tha script name that appears on the PDF you posted. YOu would need to post the script and any data files and other funciton files called.

I may not be able to answer at this time, because I am using an old computer with Matlab 2016, while my usual computer is waiting for a fix. There is a good chance that your script or function files will call some Matlab function that did not exist in 2016. Then I would not be able to run your code.

[John Amoo-Otoo]

William, Please attach the script and data

[William Rose]

@John Amoo-Otoo,

As I predicted in my previous comment, I cannot run processdata.m because it calls rmmissing(), which is not implemented in Matlab 2016a.

[John Amoo-Otoo]

Okay

[William Rose]

@John Amoo-Otoo,

I worked around the calls to rmmissing(). Then I get the same error you reported earlier*:

Index exceeds matrix dimensions.

Error in getmswtfeat (line 115)

prob = percentENER(:,st:en);

Error in processdataWR (line 37)

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

The error above occurs because percentENER is 43x52, but st=689. Line 115 tells Matlab to get data from columns 689 and following of a matrix which has only 52 columns. Of course this is an error.

The error occurs inside a funciton which you did not write (accfording to the documentaiton in the code). I suggest you contact the author of the function for assistance.

You call getmswtfeat twice: once on each of the fiinal two lines of your script. The comment preceding the calls is

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

But getmswtfeat does not return the power spectrum, so either the comment is incorrect (i.e. you are not really interesting in observing the FFT), or you are using a function that does not accomplish your goal.

You have added noise to the signals inside getmswtfeat, The signals which are analyzed and plotted in the main program are unaffected by this added noise, because 1. The noisy signals are not passed back from getmswtfeat to the main program, and 2. even if the signals were passed back, the calls to getmswtfeat are at the very end of the main program.

*The line numbers in the error message are different than what you reported, because of the work-around code which I added.

[John Amoo-Otoo]

Okay william.

I updated my script and this is what I got

Answer 2

[John Amoo-Otoo]

William, for a script like this where will the Gaussian Noise fit. I have also attached the data that needs to be processed

[William Rose]

@John Amoo-Otoo,

I downloaded the m file and data file which you posted. I see that the m file is a function, not a script. It is not runnable on its own. Can you provide a script that calls this funciton, to help me understand it? Thank you.

[William Rose]

@John Amoo-Otoo,

Question 1. Why do you want to add Gaussian noise to the signal? I am asking because you asked "where will the Gaussian Noise fit", and the answer depends on what your goal is. If you want to test the ability of an algorithm to distinguish lightning versus some other fault, when the signal is noisy, then add the noise to the signals before you process them. If you have a downstream algorithm (after the DWT) which you want to evaluate, then maybe you woudl add noise after the DWT. Therefore please explain in more detail whay you want to add noise to a signal.

Quesiton 2a, 2b. The DWT function which you posted (which appears to have been written by a third party) says

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

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

You must specify values for winsize and wininc when calling the function. Q.2a. What values do you plan to use for winsize and for wininc? Q.2b. What do you plan to do with the output from this function?

-------------

I examined the data in the file you sent, by doing the following:

clear
data=load('John_FaultData112km0ohm');
t1=table2array(data.DataPP112km0hmsFaultData(:,1));
x1=table2array(data.DataPP112km0hmsFaultData(:,2));
t2=table2array(data.SteadyStateNoneFaultState(1:end-1,1));
x2=table2array(data.SteadyStateNoneFaultState(1:end-1,2));
%Skip the last row of t2, x2, due to NaNs in that row
figure;
subplot(121); plot(t1,x1,'-r.');
xlabel('Time'); ylabel('Amplitude'); title('Transient 1'); grid on
subplot(122); plot(t2,x2,'-r.');
xlabel('Time'); ylabel('Amplitude'); title('Transient 2'); grid on

The data file includes two recordings. Record 1 is 0.7 seconds long. Record 2 is 30 seconds long. Each record shows a single transient event. Transient 1 occurs 0.2 seconds into the record. Transient 2 occurs at the start of record 2. Transient 1 recovers about 90% of the way to the pre-transient value. Transient 2 never recovers - it is a step change.

The transients differ in amplitude by a factor of 10^5. Therefore a simple threshold detector would easily distinguish the two. The non-return-to-baseline of transient 2 is another distinguishing feaure.

[John Amoo-Otoo]

Answer to Question 1: The reason why I need to add Gaussian Noise to the orginal signal before processing is to be able to distinguish between a lightning and switching signal(Non-fault) which is quite similar to a Gaussian Noise to Normal Steady state operation and and also a fault condition. Transient 1 is the Fault simulation and Transient 2 is the steady state(Non Fault conidtion)

[John Amoo-Otoo]

Values for winsize is 32 and increment is 16

[John Amoo-Otoo]

The output is to distinguish between a fault condition and Non faut condition. The entropy output(Highest and Minmimum entropy) wil tell me a fault condition and when a Protection starts and will use it as a treshold for the protection settings