How to incorporate multiple IF conditions to generate noisy data?

Hello, I wanted to simulate a signal which has noise in the form of spikes (both positive and negative). If a single IF condition is inserted, one can easily generate positive spikes of desired height. Suppose we also wish to insert an additional condition that if n (i) <0, we get a negative spike of magnitude -0.5 in the noise. I have tried several approaches but the vector size of noise increases if I include a separate IF condition. What would be an appropriate way to incorporate both outcomes of IF in one vector with the same size as "t"? Thanks.

t=0:0.01:20;
x=normpdf(t, 10, 0.1); % Gaussian peak
n=rand(1,length(t));
noise=[];
for i=1:length(t)
    if(n(i)>=1)
        noise=[noise 1.5];
    else
        noise=[noise 0];
    end
end
x=x+noise;

[Mathieu NOE]

hello

there is already a function that is doing what you are looking for

see pulstran

pulstran Pulse train generator.

pulstran generates pulse trains from either continuous functions or sampled prototype pulse

% Example 2: Generate a periodic Gaussian pulse signal at 10 kHz,

% with 50% bandwidth. The pulse repetition frequency is 1 kHz,

% sample rate is 50 kHz, and pulse train length is 10msec. The

% repetition amplitude should attenuate by 0.8 each time. Uses

% a function handle to refer to the generator function.

T = 0 : 1/50E3 : 10E-3;

D = [0 : 1/1E3 : 10E-3 ; 0.8.^(0:10)]';

Y = pulstran(T,D,@gauspuls,10E3,.5); plot(T,Y)

[FW]

Thanks, I checked the pulstran function. I wanted to generate a "single point" spike, as you can see in the code. But that code only generates a positive spike. I was wondering how to generate negative and positive single point spikes. This is simulate a liquid flow cell where a flowing bubble generates an unwanted spike.

[Mathieu NOE]

so I could not make your code to work , so I did it my own way

maybe it helps you customize the shape of your signal according to your needs

samples = 1000;

peaks_width = 10; % in samples 
data = 1 + 0.15*rand(1,samples); % signal baseline = a constant + some random noise 

% positive peaks (multiple)
time_ind = 100:200:900;
data(time_ind) = 11;

% negative peaks (multiple)
time_ind = 150:200:950;
data(time_ind) = -9;

% first pass of sliding avg : creates square pulses
out = myslidingavg(data, peaks_width);
% correct output amplitude to get the same max amplitude of peaks
cor_factor = max(abs(detrend(data)))./max(abs(detrend(out)));
out = mean(out) + (out-mean(out)).*cor_factor;

% second pass of sliding avg : creates triangular pulses
out2 = myslidingavg(out, peaks_width);
% correct output amplitude to get the same max amplitude of peaks
cor_factor = max(abs(detrend(out)))./max(abs(detrend(out2)));
out2 = mean(out2) + (out2-mean(out2)).*cor_factor;

% third pass of sliding avg : creates parabolic pulses
out3 = myslidingavg(out2, peaks_width);
% correct output amplitude to get the same max amplitude of peaks
cor_factor = max(abs(detrend(out2)))./max(abs(detrend(out3)));
out3 = mean(out3) + (out3-mean(out3)).*cor_factor;

figure(1),
plot((1:samples),data,'b',(1:samples),out,'r',(1:samples),out2,'m',(1:samples),out3,'c');

function out = myslidingavg(in, N)
%   OUTPUT_ARRAY = MYSLIDINGAVG(INPUT_ARRAY, N)
%
%  The function 'slidingavg' implements a one-dimensional filtering, applying a sliding window to a sequence. Such filtering replaces the center value in
%  the window with the average value of all the points within the window. When the sliding window is exceeding the lower or upper boundaries of the input
%  vector INPUT_ARRAY, the average is computed among the available points. Indicating with nx the length of the the input sequence, we note that for values
%  of N larger or equal to 2*(nx - 1), each value of the output data array are identical and equal to mean(in).
%
%  *  The input argument INPUT_ARRAY is the numerical data array to be processed.
%  *  The input argument N  is the number of neighboring data points to average over for each point of IN.
%
%  *  The output argument OUTPUT_ARRAY is the output data array.

if (isempty(in)) | (N<=0)                                              % If the input array is empty or N is non-positive,
 disp(sprintf('SlidingAvg: (Error) empty input data or N null.'));     % an error is reported to the standard output and the
 return;                                                               % execution of the routine is stopped.
end % if

if (N==1)                                                              % If the number of neighbouring points over which the sliding 
 out = in;                                                             % average will be performed is '1', then no average actually occur and
 return;                                                               % OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the routine
end % if                                                               % is stopped.

nx   = length(in);             % The length of the input data structure is acquired to later evaluate the 'mean' over the appropriate boundaries.

if (N>=(2*(nx-1)))                % If the number of neighbouring points over which the sliding 
 out = mean(in)*ones(size(in));   % average will be performed is large enough, then the average actually covers all the points
 return;                          % of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be gained by such an approach.
end % if                          % The execution of the routine is stopped.

out = zeros(size(in));         % In all the other situations, the initialization of the output data structure is performed.

if rem(N,2)~=1                 % When N is even, then we proceed in taking the half of it:
 m = N/2;                      % m = N     /  2.
else                           % Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we proceed taking the half of it:
 m = (N-1)/2;                  % m = (N-1) /  2.
end % if

for i=1:nx,                 % For each element (i-th) contained in the input numerical array, a check must be performed:
  dist2start = i-1;         %  index distance from current index to start index (1)
  dist2end = nx-i;          %  index distance from current index to end index (nx)
  if dist2start<m || dist2end<m        % if we are close to start / end of data, reduce the mean calculation on centered data vector reduced to available samples
      dd = min(dist2start,dist2end);    % min of the two distance (start or end)
  else
      dd = m;
  end % if
  out(i) = mean(in(i-dd:i+dd));     % mean of centered data , reduced to available samples at both ends of the data vector

end % for i

[FW]

Thanks for your efforts. normpdf may not be working because I have installed the statistics toolbox. It is just a single Gaussian peak. I could not see your plots. The error is Line: 67 Column: 12 All functions in a script must be closed with an 'end'.

[Mathieu NOE]

can you simply add one more end at the last line ? this was missing when I did the copy paste

it should end like this :

.....
    for i=1:nx,                 % For each element (i-th) contained in the input numerical array, a check must be performed:
  dist2start = i-1;         %  index distance from current index to start index (1)
  dist2end = nx-i;          %  index distance from current index to end index (nx)
  if dist2start<m || dist2end<m        % if we are close to start / end of data, reduce the mean calculation on centered data vector reduced to available samples
      dd = min(dist2start,dist2end);    % min of the two distance (start or end)
  else
      dd = m;
  end % if
  out(i) = mean(in(i-dd:i+dd));     % mean of centered data , reduced to available samples at both ends of the data vector

end % for i

end

[FW]

Okay, now I can see it. Thanks. This is somewhat close to what I would like, but it has to be random. Let me add a picture in the OP of how the spikes I was originally mentioning. I should have posted a figure first. Sorry about that. The red is the original Gaussian and I wanted to add those spikes. The current problem is that I could only generate positive spikes. I wanted some random negative spikes.

[Mathieu NOE]

ok

and the negative spike must follow the positive spike ? i mean the pulses must be pos / neg / pos / neg .... ?

both with random time distances ?

[FW]

No I was thinking about random negative spikes. Suppose we also wish to insert an additional condition that if n (i) <0, we get a negative spike of magnitude -0.5 in the noise.

[Mathieu NOE]

here as a first example , I created the positive peaks of a random time basis

samples = 1000;

peaks_width = 10; % in samples 
data = 1 + 0.15*rand(1,samples); % signal baseline = a constant + some random noise 

% positive peaks (multiple)
% time_ind = 100:200:900;
nb_of_peaks = 10;
amplitude_pos_peak = 10;
time_ind = randi([100 900],nb_of_peaks,1);
[time_ind,ind] = sort(time_ind); % sort in ascending order
data(time_ind) = mean(data)+amplitude_pos_peak;

% negative peaks (multiple)
time_ind = 150:200:950;
data(time_ind) = -9;

% first pass of sliding avg : creates square pulses
out = myslidingavg(data, peaks_width);
% correct output amplitude to get the same max amplitude of peaks
cor_factor = max(abs(detrend(data)))./max(abs(detrend(out)));
out = mean(out) + (out-mean(out)).*cor_factor;

% second pass of sliding avg : creates triangular pulses
out2 = myslidingavg(out, peaks_width);
% correct output amplitude to get the same max amplitude of peaks
cor_factor = max(abs(detrend(out)))./max(abs(detrend(out2)));
out2 = mean(out2) + (out2-mean(out2)).*cor_factor;

% third pass of sliding avg : creates parabolic pulses
out3 = myslidingavg(out2, peaks_width);
% correct output amplitude to get the same max amplitude of peaks
cor_factor = max(abs(detrend(out2)))./max(abs(detrend(out3)));
out3 = mean(out3) + (out3-mean(out3)).*cor_factor;

figure(1),
plot((1:samples),data,'b',(1:samples),out,'r',(1:samples),out2,'m',(1:samples),out3,'c');

[Mathieu NOE]

Finally , in its simplest form , this will generate timely random positive and negative pulses

the threshold value will modify the amount (density) of pulses (to increase the density , reduce the threshold)

a threshold of 1 is the upper limit but then the pulse density is zero , so you should stay below 1

samples = 1000;
signal = zeros(samples,1);
tmp = 2*(-0.5+rand(samples,1)); % a random signal with amplitude in -1 / +1 range 

% positive peaks (multiple)
threshold = 0.97;
time_ind = find(tmp>threshold);
signal(time_ind) = 5;

% negative peaks (multiple)
time_ind = find(tmp<-threshold);
signal(time_ind) = -5;

figure(1),
plot((1:samples),signal,'b');

Answer 1

[Adam Danz]

No need for a loop at all.

noise = 1.5*(n>1) + -.5*(n<0)

However, in your example, n is always between 0 and 1.

Demo

rng('default')
n = rand(1,8)*3-1

n = 1×8

1.4442 1.7174 -0.6190 1.7401 0.8971 -0.7074 -0.1645 0.6406

●

noise = 1.5*(n>1) + -.5*(n<0)

noise = 1×8

1.5000 1.5000 -0.5000 1.5000 0 -0.5000 -0.5000 0

●

[FW]

Thank you, this is the type of conditional I was looking. Simple as it can be.