Module

Medical Signal
Analysis

2022
—
Statistical Characterization of
Signal using MATLAB
—
Osmalina Nur Rahma, S.T., [Link].
 Module 2
Statistical Aspects

Introduction
Practical signal processing frequently involves statistical aspects. If you are using a sensor to
measure, say temperature, or light, or pressure, or anything else, you usually get a signal with noise
in it, and then there is a problem of noise removal. The almost immediate idea could be to apply
averaging; however, our advice is first trying to know better about the noise you have. There are
many other contexts where the data you get suffer from interference, lack of precision, variations
along time, etc. For example, suppose you want to measure the period of a pendulum using a watch:
the scientific procedure is to repeat the measurements (the values obtained will be different for
each measurement), get a data set, and then statistically process this set.
The best you can do in this example is to get a good estimate, in statistical terms. In this chapter
some aspects of probability and statistics, particularly relevant for signal processing, are selected.
First several kinds of probability density distributions are considered, and then parameters to
characterize random signals are introduced.

1. Random Variable
A random variable is a quantity whose value is not fixed but depends somehow on chance.
Typically, the value of a random variable may consist of a fixed part and a random component due
to uncertainty or disturbance. Other types of random variables take their values as a result of the
outcome of a random experiment.

2. Random Signals and Probability Density Distributions

The chief objective of this section is to introduce probability density distributions and
functions, selecting three illustrative cases: the uniform, the normal and the log-normal
distributions. The normal distribution is, in particular, a very important case.
Suppose there is a continuous random variable y(t), the distribution function Fy(v) of this
variable is the following:
𝐹𝑦 (v) = P (y(t) ≤ v), −∞ < v <∞
where P() is the probability.

The probability density function of y(t) is:

A well-known example of probability distribution
function, the so-called normal distribution, has a bell
shaped probability density function as shown in Fig. 1. In
this figure a shaded zone has been painted corresponding
to an interval [a, b] of the values that y(t) can have. The
probability of y(t) value to fall into this interval is given
by the area of the shaded zone.
The abbreviation “PDF” will be used in this book to
denote “Probability Density Function”.
Fig. 1 Probability Density Function
2.1. Random Signal with Uniform PD
A random signal taking equiprobable values in successive instants has a uniform PDF. For
example, the sequence of values that would be obtained recording the final angles (0◦.. 360◦) where
a roulette wheel stops along several runs in gambling days. Fig. 2.1 shows a random signal with
uniform PDF. It has been obtained using the rand() function provided by MATLAB. Notice that
the values are from 0 to 1.
This signal can be easily modified by adding a constant and/or multiplying by a constant: the
result will also have a uniform PDF.

Fig. 2.1 A random signal with uniform PDF

Program 2.1 Random signal with uniform PDF

% Random signal with uniform PDF

clc;clear;
fs=100; %sampling frequency in Hz
t=0:1/fs:5; %time intervals set for 5 seconds
N=length(t); %number of data points
y=rand(N,1); %random signal data set
plot(t,y,'-k'); %plots figure
xlabel('seconds');
title('random signal with uniform PDF');

An interesting way to check the quality of the MATLAB random variable generation functions
is by plotting a histogram of the signal values along time. For this purpose relatively large signal
data sets should be generated. MATLAB provides the hist() function to obtain the histogram.

Program 2.2 Histogram of a random signal with uniform PDF

% Histogram of a random signal with uniform PDF

clc;clear;
fs=100; %sampling frequency in Hz
 t=0:1/fs:100; %time intervals set for 100 seconds
N=length(t); %number of data points
y=rand(N,1); %random signal data set
v=0:0.02:1; %value intervals set
hist(y,v); colormap(cool); %plots histogram
xlabel('values');
title('Histogram of random signal with uniform PDF');

% Uniform PDF
v=0:0.01:1; %values set
ypdf=unifpdf(v,0,1); %uniform PDF
figure; plot(v,ypdf,'k'); hold on; %plots figure
axis([-0.5 1.5 0 1.1]);
xlabel('values'); title('uniform PDF');
plot([0 0],[0 1],'--k');
plot([1 1],[0 1],'--k');

2.2. Random Signal with Normal (Gaussian) PDF

As said before, the normal distribution is very important, both from theoretical and practical
points of view). Most noise and perturbation models employed in systems or automatic control
theory are of Gaussian nature. The practical reason is provided by the central limit theorem, which
in words says that if a phenomenon is the accumulation of many small additive random effects, it
tends to a normal distribution. For instance, the number of travels per day of an elevator. Fig. 2.2
shows a random signal with normal PDF. It has been obtained using the randn() function
provided byMATLAB (notice the slight name difference compared to rand(), which corresponds
to uniform PDF).
The values of the signal in Fig. 2.2 have positive and negative values. The figure has been
obtained with Program 2.3, using the randn() function, which generates a signal with mean zero,
variance one and a standard deviation one.

Figure 2.2 A random signal with normal PDF

Program 2.3 Random signal with normal PDF

% Random signal with normal PDF

clc;clear;
fs=100; %sampling frequency in Hz
t=0:1/fs:5; %time intervals set for 5 seconds
N=length(t); %number of data points
 y=randn(N,1); %random signal data set
plot(t,y,'-k'); %plots figure
xlabel('seconds');
title('random signal with normal PDF');

Program 2.4 Histogram of a random signal with normal PDF

% Histogram of a random signal with normal PDF

clc;clear;
fs=100; %sampling frequency in Hz
t=0:1/fs:100; %time intervals set for 100 seconds
N=length(t); %number of data points
y=randn(N,1); %random signal data set
v=-4:0.1:4; %value intervals set
hist(y,v); colormap(cool); %plots histogram
xlabel('values');
title('Histogram of random signal with normal PDF');

% Normal PDF
v=-3:0.01:3; %values set
mu=0; sigma=1; %random variable parameters
ypdf=normpdf(v,mu,sigma); %normal PDF
figure; plot(v,ypdf,'k'); hold on; %plots figure
axis([-3 3 0 0.5]);
xlabel('values'); title('normal PDF');

2.3. Random Signal with Long-Normal PDF

A random variable y is log-normally distributed if log(y) has a normal distribution. The log-
normal distribution is related the multiplicative product of many small independent factors. It is
observed for instance in environment, microbiology, human medicine, social sciences, or
economics contexts. For example, the case of latent periods (time from infection to first
symptoms) of infectious diseases.
Figure 2.3 shows a random signal with log-normal PDF. It has been obtained, using the
lognrnd() function, with the Program 2.5. A mean zero and a standard deviation one has been
specified inside the parenthesis of lognrnd(); other mean and standard deviation values can be
explored. Notice that signal values are always positive.

Figure 2.3 A random signal with log-normal PDF

 Program 2.5 Random signal with log-normal PDF

% Random signal with log-normal PDF

clc;clear;
fs=100; %sampling frequency in Hz
t=0:1/fs:5; %time intervals set for 5 seconds
N=length(t); %number of data points
mu=0; sigma=1; %random signal parameters
y=lognrnd(mu,sigma,N,1); %random signal data set
plot(t,y,'-k'); %plots figure
xlabel('seconds');
title('random signal with log-normal PDF');

Program 2.6 Histogram of a random signal with log-normal PDF

% Histogram of a random signal with log-normal PDF

clc;clear;
fs=100; %sampling frequency in Hz
t=0:1/fs:100; %time intervals set for 100 seconds
N=length(t); %number of data points
mu=0; sigma=1; %random signal parameters
y=lognrnd(mu,sigma,N,1); %random signal data set
v=0:0.1:12; %value intervals set
hist(y,v); colormap(cool); %plots histogram
axis([0 8 0 700])
xlabel('values');
title('Histogram of random signal with log-normal PDF');

% Log-normal PDF
v=-3:0.01:6; %values set
mu=0; sigma=1; %random variable parameters
ypdf=lognpdf(v,mu,sigma); %log-normal PDF
figure; plot(v,ypdf,'k'); hold on; %plots figure
axis([0 6 0 0.7]);
xlabel('values'); title('log-normal PDF');

3. Mean, Median, Mode, Variance, Standard Deviation

Fig. 2.4 shows a skewed PDF where the mean, the median and the mode values are marked
(Program 2.7). In symmetrical PDFs these three values would be coincident.
Consider the random variable y, with 𝐹𝑦 (v) as PDF. The mean μ of the variable y is the
expected value of y (Eq. 2.1). It is also called the average value of y. Using a mass analogy, it may
be regarded as the center of mass of the distribution.

(Eq. 2.1)

The moments about the mean, or central moments, for the variable y are given by:
(Eq. 2.2)
(μ denotes the mean of y)

A median y0 of the variable y is any point that divides the mass of the distribution into two
equal parts, that is:
 (Eq. 2.3)
A point vi such that:
(Eq. 2.4)

(where ε is an arbitrarily small positive quantity) is called a mode of y.

A mode is a value of y corresponding to a peak of the PDF. When the PDF has only one peak,
the distribution is said to be unimodal.

Fig. 2.4 Mean, median and mode marked on a PDF

In measurement tasks, depending on the PDF of the signal being obtained it would be
advisable to consider the mean, or the median, or the mode or modes, as the value of interest. In
particular, while the mean of a variable y may not exist, the median will exist.

Program 2.7 A skewed PDF with mean, median and mode with normal PDF

% A skewed PDF with mean, median and mode with normal PDF
clc;clear;
% Normal PDF
v=-3:0.01:3; %values set
mu=0; sigma=1; %random variable parameters
ypdf=normpdf(v,mu,sigma); %normal PDF
figure; plot(v,ypdf,'k'); hold on; %plots figure
axis([-3 3 0 0.5]);
xlabel('values'); title('normal PDF');

hold on;
mu=mean(y); %mean of y
vo=median(y); %median of y
mo=mode(y); %modus of y
variance=var(y) %variance of y
std_dev=std(y); %std deviation of y
[pky,pki]=max(ypdf); %peak of the PDF
plot([mu mu],[0 0.33],'--k'); %mean
plot([vo vo],[0 0.37],':k'); %median
plot([v(pki) v(pki)],[0 pky],'-.k'); %mode
 The variance of the signal y(t) is the second moment of y(t) about the mean. The positive
square root of the variance is the standard deviation σ. The variance is related with how large
is the range of values of y(t). The random signal y(t) is called strict-sense stationary if all its
statistical properties are invariant to a shift of the time origin.
MATLAB offers the following functions: mean(), median(), var(), std() (for standard
deviation).

Task 1
1. Do the Program 2.1 to 2.7. Report the plot and compare the result!
2. Check the given signal whether it is uniform PDF, normal PDF or long-normal PDF. Also
measure its statistical aspect (Mean, median, mode, variance and std deviation)!

4. Correlation, Covariance and Shifted Correlations

The basic measurements of mean, variance, standard deviation, and rms generally do not
capture the important features of a signal. For example, if we have digital versions of the EEG
signals, we can easily compute their mean, variance, standard deviation, and rms value, but these
would not tell us much about the signals or the neural processes that created them. More insight
might be gained by comparing these signals with some reference signal(s).
Comparing a signal with a reference signal, or perhaps a group or “family” of reference
signals, is an oft used tool in signal analysis. Such reference signals, or signal families, tend to be
much less complicated than the signal such as a sine wave or family of sine waves.
A quantitative comparison can tell you how much your complicated signal is like a simpler,
easier to understand reference signal or family. If enough comparisons are made with a well-
chosen reference family, these comparisons can actually provide an alternative representation of
the signal. Sometimes this new representation of the signal is more informative or enlightening
than the original.

4.1. Correlation
Correlation between different signals is illustrated in Figure 2.5, which shows various pairs
of waveforms and the correlation between them. The lack of correlation between two sinusoids
shows that correlation does not always measure general similarity. Mathematically they are as un-
alike as possible, even though they have similar behavioral patterns.
The linear correlation between two digital functions or signals can be obtained using the
Pearson correlation coefficient defined as:

(Eq. 2.5)

where rxy is a common symbol in signal analysis for the correlation between x and y. (For
the Pearson correlation coefficient, the symbol rxy is also used.) The variables x and y could
represent any two waveforms. Again x and y are the means of x and y.
Fig. 2.5 Three pairs of signals and the correlation between them as given by the Pearson correlation coefficient defined in Eq.
2.5. The high correlation between the sine wave and triangular wave (center) correctly expresses the similarity between them,
but the zero correlation between the sinusoids in the upper plot does not reflect their general similarity.

Example 2.1

Find the correlation between a sine wave and a square compare to the correlation beetween
a cosine wave and a square. All waveforms have amplitudes of 1.0 V (peak-to-peak) and periods
of 4.0 s.

Program 2.8 Correlation between two signals

% Correlation between two signals

clc;clear;
N = 500; % Number of points
Tt = 4.0 % Desired total time
f = 0.25; % Wave frequency in Hz
fs = N/Tt; % Calculate sampling frequency
t = (0:N-1)/fs; % Time vector from 0 (approx.) to 4 sec
x = sin(2*pi*f*t); % 0.25 Hz sine wave
y = cos(2*pi*f*t); % 0.25 Hz cosine wave
z = [ones(1,N/2) -ones(1,N/2)]; % 0.25 Hz square wave
rxz = mean(x.*z); % Correlation (Eq. 2.5) x and z
rxy = mean(x.*y); % Correlation x and y
plot(t,x,'r')
hold on
plot(t,y,'y')
plot(t,z,'g')
disp([rxz rxy]) % Output correlation
 Based on Program 2.8, we could find correlation between a sine and square wave but none
between a cosine and square wave. This is another example of how correlation does not always
represent [Link] find out how to get around this problem using a technique called
“cross-correlation” (see sub 4.3).
However, any correlation method, the Pearson correlation or the unnormalized versions,
could test if two signals are orthogonal. Orthogonal signals and functions can be very useful
signal processing tools. In common usage, “orthogonal” means perpendicular: if two lines are
orthogonal, they are perpendicular. If the basic signals are orthogonal, you can determine each
one independently without worrying about the other signals in the collection. Orthogonality
simplifies many calculations where multiple signals are involved. Some analyses could not be
done, at least not practically, using nonorthogonal signals. Orthogonality is not limited to just
two signals. Whole families exist where each signal is orthogonal to all other members in the
family. Such families of orthogonal signals are called “orthogonal sets.”

4.2. Covariance
The covariance is the correlation normalized by N given by Eq. 2.5, whereas the matrix of
correlations is the Pearson correlation normalized as in Eq. 2.6.

(Eq. 2.6)

Covariance computes the variance that is shared between two (or more) signals.
Covariance is usually defined in discrete notation as:

(Eq. 2.7)

Example 2.2

Determine if 1.0 Hz sine and cosine waves are orthogonal and if a 2.0 Hz cosine wave is
orthogonal to a 1.0 Hz cosine wave. Of course, we already know the answer, but this example
shows how well the cov and corrcoef routines perform. The 1 and 2 Hz cosine waves are called
“harmonically related,” as they have frequencies that are multiples. Also determine if a 1.0 Hz
sawtooth is orthogonal to the sinusoidal waveforms. Make the peak-to-peak amplitude of the four
signals ±1.0.

Program 2.9 Correlation and covariance between two signals

% Application of the covariance matrices to sinusoids that are orthogonal

%and a sawtooth
clear;clc;
N=1000; % Number of points
Tt=2; % desired total time
fs=N/Tt; % Calculate sampling frequency
t=(0:N-1)/fs; % Time vector from 0 (approx.) to 2 sec
x(:,1)=cos(2*pi*t)'; % Generate a 1 Hz cosine
x(:,2)=sin(2*pi*t)'; % Generate a 1 Hz sine
x(:,3)=cos(4*pi*t)'; % Generate a 2 Hz cosine
x(:,4)=sawtooth(2*pi*t)'; % Generate a 1 Hz sawtooth
S=cov(x) % Print covariance matrix
Rxx=corrcoef(x) % and correlation matrix
 Example 2.2 uses covariance and correlation analysis to determine if sines and cosines of
different frequencies are orthogonal. Again, two orthogonal signals have zero correlation. Either
covariance or correlation could be used to determine if signals are orthogonal.

4.3. Shifted Correlations: Cross-Correlation

Many of the real-world signals we are called upon to analyze are quite complex. One way to
deal with these complex signals is to see how much they are like less complicated signals by using
our new-found correlation tools and comparing the EEG signal with some less complicated signals
such as sinusoids. Comparing with sinusoids has the advantage as it is easy to interpret the results;
sinusoids represent oscillatory behavior so, when we compare a signal with a sinusoid at a given
frequency, we are actually searching for oscillatory behavior at that frequency. But the lack of
correlation between sine and cosine at the same frequency becomes troubling. Fig. 2.6 illustrates
the problem finding the correlation between a cosine and a sine with a phase shift (i.e., a sinusoid).
The correlation depends on the phase shift of the sine wave. In Figure 2.12A the sine is not
shifted and there is zero correlation, but when the sine is shifted by 45 degrees (Fig. 2.6B) the
correlation coefficient is 0.71, and at a 90-degree shift, the correlation coefficient is 1.0 (Figure
2.6C). Figure 2.6D shows that the correlation as a function of sinusoidal shift is itself a sine wave
ranging between ±1.

Fig. 2.5 A) The correlation between a 2 Hz cosine reference (dashed) and an unshifted 2.0 Hz sine wave is
zero. (B) When the sine wave is time shifted by the equivalent of 45 degree, the Pearson correlation is 0.71. (C) When
the sine wave is time shifted equal to 90 degree the sine wave becomes a cosine wave and the correlation is 1.0. (D)
Plotting the correlation as a function of the time shift shows a sinusoidal variation. The peak value, 1.0, comes at
0.125 s, which happens to be a phase shift of 90 degrees. At a time shift of 0.25 s the sine wave is a sine wave again
(but inverted) so the correlation is back to zero. At a time shift of 0.375 s the sine wave (now shifted by 270 degree)
becomes an inverted cosine wave so it has a correlation of ±1.0 with a cosine wave.
 To restate the problem, when using a reference signal to search for a particular behavior we
could miss it depending on the time position (or phase) of the reference signal. Figure 2.6D shows
correlations as a function of time shift of the sine wave reference signal and shows that at some
shift (0.125 s in this particular situation) it is an exact match. Whenever we search a target signal
for a particular behavior using a reference signal, we could try correlation at a bunch of different
time shifts and take the maximum correlation as representing the true similarity between reference
and signal. If we are to apply this shifting/correlation approach, we need to decide how much to
shift the reference signal between correlations. We do not want to shift the reference by so much
that we pass over the shift that gives maximum correlation between the reference and signal.

Example 2.3

Find the correlation between the EEG signal given in [Link] and
sinusoids at 6.5 and 14 Hz. The EEG data is stored in file CE1_1.mat and the signal was sampled
at 100 Hz. Is the oscillatory behavior greater at 14 Hz sinusoid or 6.5 Hz? Also generate two plots
of the EEG signal superimposed on each of the two sinusoids shifted for best correlation. Scale
the plots to display all the signals nicely.

Program 2.10 Correlation two signals

close all;clear;
clc;
load CE1_1; % Get the EEG data (in variable 'eeg')
eeg=CE(:,1)';
N = length(eeg); % Number of EEG data points
fs = 512; % Sampling frequency of EEG data (given)
f = [6.5 14]; % Frequencies of reference signals
t = (0:N-1)/fs; % Time vector
figure;plot(t,eeg);
for k1 = 1:2
x = cos(2*pi*f(k1)*t); % Generate reference signal at desired frequency
figure;plot(t,x);
for k = 2:N % Perform N-1 correlations
y = [x(k:end), x(1:k-1)]; % Shift reference circularly
rxy(k) = mean(eeg.*y); % Correlation as a function of shift k
end
[corr(k1),shift(k1)] = max(rxy); % Find maximum correlation
end
%
% Plotting section
for k = 1:2
subplot(2,1,k); % Plot the two sine waves separately
x = cos(2*pi*f(k)*t); % Recreate the reference signal for plotting
y = [x(shift(k):end), x(1:shift(k)-1)]; % Shift reference
y = y*(max(eeg)/max(y))/4; % Scale the sinusoid for good viewing
plot(t,eeg); hold on; % Plot EEG
plot(t,y); % Plot shifted reference
% xlim([2 2.6]); % Scale time axis for better viewing
... labels and text ...
end
 In the Signal Processing Toolbox, MATLAB features a routine that performs crosscorrelation
called xcorr. The calling structure is:
[rxy lags] = xcorr(x,y,maxlags); % Perform crosscorrelation

where x and y are the target and reference signals (makes no difference which is which) and
maxlags is an optional argument indicating the maximum number of shifts (the default is 2 N - 1
where N is the length of the larger input signal). The routine uses zero padding to generate the
additional points. The cross-correlation found in rxy and lags is a vector the same length as rxy and
contains the corresponding lags. This is useful in finding the lag that corresponds to a given
correlation (such as the maximum correlation) or for plotting.
The xcorr routine has a lot of bells and whistles, including options such as various ways to
bias the correlation, but we can produce a similar routine by revising the code in Example 2.3. So
in the next example we modify the code in Example 2.3 and make a routine similar to xcorr. We
use that routine to find the correlation between sinusoidal reference signals and the EEG signal,
but rather than just compare at two sinusoids, let us make the comparison over a range of sinusoids
from 1 to 25 Hz in 0.5-Hz increments. We also compare our routine to MATLAB's xcorr for these
same reference signals.

Example 2.4

Find the correlation between the EEG signal given in [Link] and
sinusoids at 6.5 and 14 Hz. The EEG data is stored in file CE1_1.mat and the signal was sampled
at 100 Hz. Is the oscillatory behavior greater at 14 Hz sinusoid or 6.5 Hz? Also generate two plots
of the EEG signal superimposed on each of the two sinusoids shifted for best correlation. Scale
the plots to display all the signals nicely.

Program 2.11 Cross-Correlation over a range of sinusoids from 1 to 25 Hz in 0.5-Hz

close all;clear;
clc;
load CE1_1; % Get the EEG data (in variable 'eeg')
eeg=CE(:,1)';
N = length(eeg); % Number of EEG data points
fs = 512; % Sampling frequency of EEG data (given)
t = (0:N-1)/fs; % Time vector
f = (1:0.5:30); % Frequencies of reference signals

figure;plot(t,eeg);
for k = 1:length(f)
x = cos(2*pi*f(k)*t); % Generate reference signal at desired frequency
rxy = crosscorr(x,eeg); % Compute crosscorrelation
max_corr(k)= max(rxy); % Find maximum correlation
rxy_M = xcorr(x,eeg,'biased'); % Find MATLAB crosscorrelation
max_corr_M(k) = max(rxy_M); % MATLAB's max correlation
end
figure;plot(f,max_corr,'k'); hold on; % Plot crosscorr results
plot(f,max_corr_M,'*k'); % Plot MATLAB results as *-points
......labels.......
Task 2
1. Do the Program 2.8 to 2.11. Report the plot and compare the result!
2. Analyse the result of Program 2.10. Which one is giving the better correlation between two
sinusoids shifted?
3. Analyse the result of Program 2.11. Which frequency of sinusoids shows the maximum cross-
correlation between a reference sinusoid and the EEG signal?
4. Is the cross-correlation function of two random Gaussian variables (N = 500) itself a Gaussian
random variable? Show. (Make your evidence definitive!)