GaussianWaves

Plot FFT using Matlab – FFT of sine wave & cosine wave

Key focus: Learn how to plot FFT of sine wave and cosine wave using Matlab. Understand FFTshift. Plot one-sided, double-sided and normalized spectrum.

Introduction

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT). Often we are confronted with the need to generate simple, standard signals (sine, cosine, Gaussian pulse, squarewave, isolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Matlab and how to represent them in frequency domain using FFT. If you are inclined towards Python programming, visit here.

This article is part of the following books
● Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
● Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Sine Wave

In order to generate a sine wave in Matlab, the first step is to fix the frequency of the sine wave. For example, I intend to generate a f=10 Hz sine wave whose minimum and maximum amplitudes are and respectively. Now that you have determined the frequency of the sinewave, the next step is to determine the sampling rate. Matlab is a software that processes everything in digital. In order to generate/plot a smooth sine wave, the sampling rate must be far higher than the prescribed minimum required sampling rate which is at least twice the frequency – as per Nyquist Shannon Theorem. A oversampling factor of is chosen here – this is to plot a smooth continuous-like sine wave (If this is not the requirement, reduce the oversampling factor to desired level). Thus the sampling rate becomes . If a phase shift is desired for the sine wave, specify it too.

f=10; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
phase = 1/3*pi; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=sin(2*pi*f*t+phase); %replace with cos if a cosine wave is desired
plot(t,x);
title(['Sine Wave f=', num2str(f), 'Hz']);
xlabel('Time(s)');
ylabel('Amplitude');

Representing in Frequency Domain

Representing the given signal in frequency domain is done via Fast Fourier Transform (FFT) which implements Discrete Fourier Transform (DFT) in an efficient manner. Usually, power spectrum is desired for analysis in frequency domain. In a power spectrum, power of each frequency component of the given signal is plotted against their respective frequency. The command computes the -point DFT. The number of points – – in the DFT computation is taken as power of (2) for facilitating efficient computation with FFT. A value of is chosen here. It can also be chosen as next power of 2 of the length of the signal.

Different representations of FFT:

Since FFT is just a numeric computation of -point DFT, there are many ways to plot the result.

1. Plotting raw values of DFT:

The x-axis runs from to – representing sample values. Since the DFT values are complex, the magnitude of the DFT is plotted on the y-axis. From this plot we cannot identify the frequency of the sinusoid that was generated.

NFFT=1024; %NFFT-point DFT      
X=fft(x,NFFT); %compute DFT using FFT        
nVals=0:NFFT-1; %DFT Sample points       
plot(nVals,abs(X));      
title('Double Sided FFT - without FFTShift');        
xlabel('Sample points (N-point DFT)')        
ylabel('DFT Values');

2. FFT plot – plotting raw values against Normalized Frequency axis:

In the next version of plot, the frequency axis (x-axis) is normalized to unity. Just divide the sample index on the x-axis by the length of the FFT. This normalizes the x-axis with respect to the sampling rate . Still, we cannot figure out the frequency of the sinusoid from the plot.

NFFT=1024; %NFFT-point DFT      
X=fft(x,NFFT); %compute DFT using FFT        
nVals=(0:NFFT-1)/NFFT; %Normalized DFT Sample points         
plot(nVals,abs(X));      
title('Double Sided FFT - without FFTShift');        
xlabel('Normalized Frequency')       
ylabel('DFT Values');

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

As you know, in the frequency domain, the values take up both positive and negative frequency axis. In order to plot the DFT values on a frequency axis with both positive and negative values, the DFT value at sample index has to be centered at the middle of the array. This is done by using function in Matlab. The x-axis runs from to where the end points are the normalized ‘folding frequencies’ with respect to the sampling rate .

NFFT=1024; %NFFT-point DFT      
X=fftshift(fft(x,NFFT)); %compute DFT using FFT      
fVals=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points        
plot(fVals,abs(X));      
title('Double Sided FFT - with FFTShift');       
xlabel('Normalized Frequency')       
ylabel('DFT Values');

4. FFT plot – Absolute frequency on the x-axis Vs Magnitude on Y-axis:

Here, the normalized frequency axis is just multiplied by the sampling rate. From the plot below we can ascertain that the absolute value of FFT peaks at and . Thus the frequency of the generated sinusoid is . The small side-lobes next to the peak values at and are due to spectral leakage.

NFFT=1024;      
X=fftshift(fft(x,NFFT));         
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,abs(X),'b');      
title('Double Sided FFT - with FFTShift');       
xlabel('Frequency (Hz)')         
ylabel('|DFT Values|');

5. Power Spectrum – Absolute frequency on the x-axis Vs Power on Y-axis:

The following is the most important representation of FFT. It plots the power of each frequency component on the y-axis and the frequency on the x-axis. The power can be plotted in linear scale or in log scale. The power of each frequency component is calculated as

Where is the frequency domain representation of the signal . In Matlab, the power has to be calculated with proper scaling terms (since the length of the signal and transform length of FFT may differ from case to case).

NFFT=1024;
L=length(x);         
X=fftshift(fft(x,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,Px,'b');      
title('Power Spectral Density');         
xlabel('Frequency (Hz)')         
ylabel('Power');

If you wish to verify the total power of the signal from time domain and frequency domain plots, follow this link.
Plotting the power spectral density (PSD) plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.

NFFT=1024;      
L=length(x);         
X=fftshift(fft(x,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,10*log10(Px),'b');        
title('Power Spectral Density');         
xlabel('Frequency (Hz)')         
ylabel('Power');

6. Power Spectrum – One-Sided frequencies

In this type of plot, the negative frequency part of x-axis is omitted. Only the FFT values corresponding to to sample points of -point DFT are plotted. Correspondingly, the normalized frequency axis runs between to . The absolute frequency (x-axis) runs from to .

L=length(x);        
NFFT=1024;       
X=fft(x,NFFT);       
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(0:NFFT/2-1)/NFFT;      
plot(fVals,Px(1:NFFT/2),'b','LineSmoothing','on','LineWidth',1);         
title('One Sided Power Spectral Density');       
xlabel('Frequency (Hz)')         
ylabel('PSD');

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

For further reading

[1] Power spectral density – MIT opencourse ware↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
□ Sinusoidal signals
□ Square wave
□ Rectangular pulse
□ Gaussian pulse
□ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
□ Real and complex DFT
□ Fast Fourier Transform (FFT)
□ Interpreting the FFT results
□ FFTShift
□ IFFTShift
● Obtaining magnitude and phase information from FFT
□ Discrete-time domain representation
□ Representing the signal in frequency domain using FFT
□ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
□ Energy of a signal
□ Power of a signal
□ Classification of signals
□ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
□ Polynomial functions
□ Representing single variable polynomial functions
□ Multiplication of polynomials and linear convolution
□ Toeplitz matrix and convolution
● Methods to compute convolution
□ Method 1: Brute-force method
□ Method 2: Using Toeplitz matrix
□ Method 3: Using FFT to compute convolution
□ Miscellaneous methods
● Analytic signal and its applications
□ Analytic signal and Fourier transform
□ Extracting instantaneous amplitude, phase, frequency
□ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
□ Design specification
□ General considerations in design

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Generate multiple sequences of correlated random variables

In the previous post, a method for generating two sequences of correlated random variables was discussed. Generation of multiple sequences of correlated random variables, given a correlation matrix is discussed here.

Correlation Matrix

Correlation matrix defines correlation among N variables. It is a symmetric matrix with the element equal to the correlation coefficient between the and the variable. The diagonal elements (correlations of variables with themselves) are always equal to 1.

Sample problem:

Let’s say we would like to generate three sets of random sequences X,Y,Z with the following correlation relationships.

Correlation co-efficient between X and Y is 0.5
Correlation co-efficient between X and Z is 0.3
Obviously the variable X correlates with itself 100% – i.e, correlation-coefficient is 1

Putting all these relationships in a compact matrix form, gives the correlation matrix. We take arbitrary correlation value (0.3) for the relationship between Y and Z.

Now, the task is to generate three sets of random numbers X,Y and Z that follows the relationship above. The problem can be addressed in many ways. Two most common methods finding the solution are

Cholesky Decomposition method
Spectral Decomposition ( also called Eigen Vector Decomposition) method

The Cholesky Decomposition method is discussed here.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

Generating Correlated random number using Cholesky Decomposition:

Cholesky decomposition is the matrix equivalent of taking square root operation on a given matrix. As with any scalar values, positive square root is only possible if the given number is a positive (Imaginary roots do exist otherwise). Similarly, if a matrix need to be decomposed into square-root equivalent, the matrix need to be positive definite.

The method discussed here, seeks to decompose the given correlation matrix using Cholesky decomposition.

where U and L are upper and lower triangular matrices. We will consider Upper triangular matrix here. Equivalently, lower triangular matrix can also be used, in which case the order of output needs to be reversed.

For this decomposition to work, the correlation matrix should be positive definite. The correlated random sequences (where X,Y,Z are column vectors) that follow the above relationship can be generated by multiplying the uncorrelated random numbers R with U .

Steps to follow:

Generate three sequences of uncorrelated random numbers R – each drawn from a normal distribution. For this case, the R matrix will be of size where k is the number of samples we wish to generate and we allocate the k samples in three columns, where the columns indicate the place holder for each variable X, Y and Z. Multiply this matrix with the Cholesky decomposed upper triangular version of the correlation matrix.

Python code

import numpy as np
from scipy.linalg import cholesky
from scipy.stats import pearsonr #to calculate correlation coefficient

#for plotting and visualization
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('fivethirtyeight')
import seaborn as sns

C = np.array([[1, -0.5, 0.3], 
              [-0.5, 1, 0.2],
              [0.3, 0.2, 1]]) #Construct correlation matrix
U = cholesky(C) #Cholesky decomposition
R = np.random.randn(10000,3) #Three uncorrelated sequences
Rc = R @ U #Array of correlated random sequences

#compute and display correlation coeff from generated sequences
def pearsonCorr(x, y, **kws): 
    (r, _) = pearsonr(x, y) #returns Pearson’s correlation coefficient, 2-tailed p-value)
    ax = plt.gca()
    ax.annotate("r = {:.2f} ".format(r),xy=(.7, .9), xycoords=ax.transAxes)
    
#Visualization
df = pd.DataFrame(data=Rc, columns=['X','Y','Z'])
graph = sns.pairplot(df)
graph.map(pearsonCorr)

Figure 1: Pairplot of correlated random variables generated using Cholesky decomposition (Python)

Matlab code

x=[  1  0.5 0.3; 0.5  1  0.3; 0.3 0.3  1 ;]; %Correlation matrix
U=chol(x); %Cholesky decomposition 

R=randn(10000,3); %Random data in three columns each for X,Y and Z
Rc=R*U; %Correlated matrix Rc=[X Y Z]

%Verify Correlation-Coeffs of generated vectors
coeffMatrixXX=corrcoef(Rc(:,1),Rc(:,1));
coeffMatrixXY=corrcoef(Rc(:,1),Rc(:,2));
coeffMatrixXZ=corrcoef(Rc(:,1),Rc(:,3));

%Extract the required correlation coefficients
coeffXX=coeffMatrixXX(1,2) %Correlation Coeff for XX;
coeffXY=coeffMatrixXY(1,2) %Correlation Coeff for XY;
coeffXZ=coeffMatrixXZ(1,2) %Correlation Coeff for XZ;

%Scatterplots
subplot(3,1,1)
plot(Rc(:,1),Rc(:,1),'b.')
title(['Scatterd Plot - X and X calculated \rho=' num2str(coeffXX)])
xlabel('X')
ylabel('X')

subplot(3,1,2)
plot(Rc(:,1),Rc(:,2),'r.')
title(['Scatterd Plot - X and Y calculated \rho=' num2str(coeffXY)])
xlabel('X')
ylabel('Y')

subplot(3,1,3)
plot(Rc(:,1),Rc(:,3),'m.')
title(['Scatterd Plot - X and Z calculated \rho=' num2str(coeffXZ)])
xlabel('X')
ylabel('Z')

Scattered plots to verify the simulated data

Figure 2: Scatter plots showing the correlation trend between the three random variables (Matlab)

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Random Variables - Simulating Probabilistic Systems
● Introduction
● Plotting the estimated PDF
● Univariate random variables
□ Uniform random variable
□ Bernoulli random variable
□ Binomial random variable
□ Exponential random variable
□ Poisson process
□ Gaussian random variable
□ Chi-squared random variable
□ Non-central Chi-Squared random variable
□ Chi distributed random variable
□ Rayleigh random variable
□ Ricean random variable
□ Nakagami-m distributed random variable
● Central limit theorem - a demonstration
● Generating correlated random variables
□ Generating two sequences of correlated random variables
□ Generating multiple sequences of correlated random variables using Cholesky decomposition
● Generating correlated Gaussian sequences
□ Spectral factorization method
□ Auto-Regressive (AR) model

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Generate two correlated random sequences

Note: There is a rating embedded within this post, please visit this post to rate it.

This article discusses the method of generating two correlated random sequences using Matlab. If you are looking for the method on generating multiple sequences of correlated random numbers, I urge you to go here.

Generating two vectors of correlated random numbers, given the correlation coefficient , is implemented in two steps. The first step is to generate two uncorrelated random sequences from an underlying distribution. Normally distributed random sequences are considered here.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

Step 1: Generate two uncorrelated Gaussian distributed random sequences

x1=randn(1,100); %Normal random numbers sequence 1
x2=randn(1,100); %Normal random numbers sequence 2
subplot(1,2,1); plot(x1,x2,'r*');
title('Uncorrelated RVs X_1 and X_2');
xlabel('X_1'); ylabel('X_2');

Step 2: Generate correlated random sequence z

In the second step, the required correlated sequence is generated as

rho=0.9;
z=rho*x1+sqrt(1-rhoˆ2)*x2;%transformation
subplot(1,2,2); plot(x1,z,'r*');
title(['Correlated RVs X_1 and Z , \rho=',num2str(rho)]);
xlabel('X_1'); ylabel('Z');

The resulting sequence Z will have correlation with respect to

Results plotted below.

*Figure : Scatter plots – Correlated random variables and on right*

Continue reading this article on the method to generate multiple vectors of correlated random numbers.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Sampling in Matlab and downsampling an audio file

Generating a continuous signal and sampling it at a given rate is demonstrated here. In simulations, we may require to generate a continuous time signal and convert it to discrete domain by appropriate sampling.

For baseband signal, the sampling is straight forward. By Nyquist Shannon sampling theorem, for faithful reproduction of a continuous signal in discrete domain, one has to sample the signal at a rate higher than at-least twice the maximum frequency contained in the signal (actually, it is twice the one-sided bandwidth occupied by a real signal. For a baseband signal bandwidth ( to ) and maximum frequency in a given band are equivalent).

Matlab or any other simulation softwares process everything in digital i.e, discrete in time. This is because, the signals are represented as discrete samples in computer memory. Therefore, we cannot generate a real continuous-time signal on it, rather we can generate a “continuous-like” signal by using a very very high sampling rate. When plotted, such signals look like a continuous signal.

Let’s generate a simple continuous-like sinusoidal signal with frequency . In order to make it appear as a continuous signal when plotting, a sampling rate of is used.

fs=500e3; %Very high sampling rate 500 kHz
f=10e3; %Frequency of sinusoid
nCyl=5; %generate five cycles of sinusoid
t=0:1/fs:nCyl*1/f; %time index
x=cos(2*pi*f*t);

plot(t,x)
title('Continuous sinusoidal signal');
xlabel('Time(s)');
ylabel('Amplitude');

*Figure 1: Continuous time signal in Matlab*

Pretending the above generated signal as a “continuous” signal, we would like to convert the signal to discrete-time equivalent by sampling. By Nyquist Shannon Theorem, the signal has to be sampled at at-least . Let’s sample the signal at and then at for illustration.

fs1=30e3; %30kHz sampling rate
t1=0:1/fs1:nCyl*1/f; %time index
x1=cos(2*pi*f*t1);

fs2=50e3; %50kHz sampling rate
t2=0:1/fs2:nCyl*1/f; %time index
x2=cos(2*pi*f*t2);

subplot(2,1,1);
plot(t,x);
hold on;
stem(t1,x1);
subplot(2,1,2);
plot(t,x);
hold on;
stem(t2,x2);

*Figure 2: Sampling a Continuous time signal in Matlab*

Manipulating audio files in Matlab

Matlab’s standard installation comes with a set of audio files. The audio files,that can be considered as one-dimensional vectors, can be inspected and played using xpsound command. With this command, we can visualize the audio files in three ways

● Time series (data-vector as function of time)
● Power spectral density (distribution of frequency content)
● Spectrogram (frequency content as function of time)

The output of the xpsound command plotting time-series plot of a sample audio file looks like this

We can also load and plot the time-series plot using inbuilt Matlab commands as follows

>> load('gong') %load the variables for the 'gong' audio file, this loads the sample frequency and the sample values
>> Fs  %sampling frequency

Fs =

        8192

>> y(1:10) %first 10 sample values in the file

ans =

   -0.0027
   -0.0045
   -0.0074
   -0.0110
   -0.0128
   -0.0173
   -0.0223
   -0.0223
   -0.0200
   -0.0092
>> t=0:1/Fs:length(y)/Fs-1/Fs; %time index
>> figure;plot(t,y);xlabel('Time (s)'),ylabel('y')

*Figure 4: ‘Gong’ sound visualized as time-series in Matlab*

Example: DownSampling

In signal processing, downsampling is the process of throwing away samples without applying any low-pass filtering. Mathematically, downsampling by a factor of implies, starting from the very first sample we throw away every $M-1$ samples (i.e, keep every -th sample.

For example, if is a vector of input samples, downsampling by implies

Going back to the previous example of ‘gong’ audio vector loaded in the Matlab variable space, the downsampling operation can be coded as follows.

>>M=2 % downsample by 2
>>y_down = y(1:M:end); % keep every M-th sample

Note: Downsampling↗ is not same as decimation. Decimation↗ implies reducing the sampling rate of a signal by applying low-pass filtering as the first step and then followed by decimation.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Reference

[1] Julius O smith III, “Spectral audio signal processing”, Center for Computer Research in Music and Acoustics, Stanford.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Power Delay Profile

Power delay profile gives the signal power received on each multipath as a function of the propagation delays of the respective multipaths.

Power delay profile (PDP)

A multipath channel can be characterized in multiple ways for deterministic modeling and power delay profile (PDP) is one such measure. In a typical PDP plot, the signal power on each multipath is plotted against their respective propagation delays.

In a typical PDP plot, the signal power () of each multipath is plotted against their respective propagation delays (). A sample power delay profile plot, shown in Figure 1, indicates how a transmitted pulse gets received at the receiver with different signal strengths as it travels through a multipath channel with different propagation delays. PDP is usually supplied as a table of values, obtained from empirical data and it serves as a guidance to system design. Nevertheless, it is not an accurate representation of the real environment in which the mobile is destined to operate at.

Figure 1: A typical discrete power delay profile plot for a multipath channel with 3 paths

The PDP, when expressed as an intensity function , gives the signal intensity received over a multipath channel as a function of propagation delays. The PDP plots, like the one shown in Figure 1, can be obtained as the spatial average of the complex channel impulse response as

RMS delay spread and mean delay

The RMS delay spread and mean delay are two most important parameters that characterize a frequency selective channel. They are derived from power delay profile. The delay spread of a multipath channel at any time instant, is a measure of duration of time over which most of the symbol energy from the transmitter arrives the receiver.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

In the wide-sense stationary uncorrelated scattering (WSSUS) channel models [1], the delays of received waves arriving at a receive antenna are treated as uncorrelated. Therefore, for the WSSUS model, the underlying complex process is assumed as zero-mean Gaussian random proces and hence the RMS value calculated from the normalized PDP corresponds to standard deviation of PDP distribution.

Figure 2: Relation between scattering function, power delay profile, Doppler power spectrum, spaced frequency
correlation function and spaced time correlation function

For continuous PDP (as in Figure 2), the RMS delay spread () can be calculated as

where, the mean delay is given by

For discrete PDP (as in Figure 1), the RMS delay spread () can be calculated as

where, is the power of the path, is the delay of the path and the mean delay is given by

Knowledge of the delay spread is essential in system design for determining the trade-off between the symbol rate of the system and the complexity of the equalizers at the receiver. The ratio of RMS delay spread () and symbol time duration () quantifies the strength of intersymbol interference (ISI). Typically, when the symbol time period is greater than 10 times the RMS delay spread, no ISI equalizer is needed in the receiver. The RMS delay spread obtained from the PDP must be compared with the symbol duration to arrive at this conclusion.

Frequency selective and non-selective channels

With the power delay profile, one can classify a multipath channel into frequency selective or frequency non-selective category. The derived parameter, namely, the maximum excess delay together with the symbol time of each transmitted symbol, can be used to classify the channel into frequency selective or non-selective channel.

PDP can be used to estimate the average power of a multipath channel, measured from the first signal that strikes the receiver to the last signal whose power level is above certain threshold. This threshold is chosen based on receiver design specification and is dependent on receiver sensitivity and noise floor at the receiver.

Maximum excess delay, also called maximum delay spread, denoted as (), is the relative time difference between the first signal component arriving at the receiver to the last component whose power level is above some threshold. Maximum delay spread () and the symbol time period () can be used to classify a channel into frequency selective or non-selective category. This classification can also be done using coherence bandwidth (a derived parameter from spaced frequency correlation function which in turn is the frequency domain representation of power delay profile).

Maximum excess delay is also an important parameter in mobile positioning algorithm. The accuracy of such algorithm depends on how well the maximum excess delay parameter conforms with measurement results from actual environment. When a mobile channel is modeled as a FIR filter (tapped delay line implementation), as in CODIT channel model [2], the number of taps of the FIR filter is determined by the product of maximum excess delay and the system sampling rate. The cyclic prefix in a OFDM system is typically determined by the maximum excess delay or by the RMS delay spread of that environment [3].

A channel is classified as frequency selective, if the maximum excess delay is greater than the symbol time period, i.e, . This introduces intersymbol interference into the signal that is being transmitted, thereby distorting it. This occurs since the signal components (whose powers are above either a threshold or the maximum excess delay), due to multipath, extend beyond the symbol time. Intersymbol interference can be mitigated at the receiver by an equalizer.

In a frequency selective channel, the channel output can be expressed as the convolution of input signal and the channel impulse response , plus some noise .

On the other hand, if the maximum excess delay is less than the symbol time period, i.e, , the channel is classified as frequency non-selective or zero-mean channel. Here, all the scattered signal components (whose powers are above either a specified threshold or the maximum excess delay) due to the multipath, arrive at the receiver within the symbol time. This will not introduce any ISI, but the received signal is distorted due to inherent channel effects like SNR condition. Equalizers in the receiver are not needed. A time varying non-frequency selective channel is obtained by assuming the impulse response . Thus the output of the channel can be expressed as

Therefore, for a frequency non-selective channel, the channel output can be expressed simply as product of time varying channel response and the input signal. If the channel impulse response is a deterministic constant, i.e, time invariant, then the non-frequency selective channel is expressed as follows by assuming

This is the simplest situation that can occur. In addition to that, if the noise in the above equation is white Gaussian noise, the channel is called additive white Gaussian noise (AWGN) channel.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] P. A. Bello, Characterization of randomly time-variant linear channels, IEEE Trans. Comm. Syst., vol. 11, no. 4, pp. 360–393, Dec. 1963.↗
[2] Andermo, P.G and Larsson, G., Code division testbed, CODIT, Universal Personal communications, 1993. Personal Communications, Gateway to the 21st Century. Conference Record., 2nd International Conference on , vol.1, no., pp.397,401 vol.1, 12-15 Oct 1993.↗
[3] Huseyin Arslan, Cognitive Radio, Software Deﬁned Radio, and Adaptive Wireless Systems, pp. 238, 2007, Dordrecht, Netherlands, Springer.↗

Topics in this chapter

Small-scale Models for Multipath Effects
● Introduction
● Statistical characteristics of multipath channels
□ Mutipath channel models
□ Scattering function
□ Power delay profile
□ Doppler power spectrum
□ Classification of small-scale fading
● Rayleigh and Rice processes
□ Probability density function of amplitude
□ Probability density function of frequency
● Modeling frequency flat channel
● Modeling frequency selective channel
□ Method of equal distances (MED) to model specified power delay profiles
□ Simulating a frequency selective channel using TDL model

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Multipath channel models: scattering function

Understand various characteristics of a wireless channel through multipath channel models. Discuss Wide Sense Stationary channel, uncorrelated scattering channel, wide sense stationary uncorrelated scattering channel models and scattering function.

Introduction

Wireless channel is of time-varying nature in which the parameters randomly change with respect to time. Wireless channel is very harsh when compared to AWGN channel model which is often considered for simulation and modeling. Understanding the various characteristics of a wireless channel and understanding their physical significance is of paramount importance. In these series of articles, I seek to expound various statistical characteristics of a multipath wireless channel by giving more importance to the concept than the mathematical derivation.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

Complex baseband mutipath channel model:

In a multipath channel, multiple copies of a signal travel different paths with different propagation delays and are received at the receiver at different phase angles and strengths. These rays add constructively or destructively at the receiver front end, thereby giving rise to rapid fluctuations in the channel. The multipath channel can be viewed as a linear time variant system where the parameters change randomly with respect to time. The channel impulse response is a two dimensional random variable – that is a function of two parameters – instantaneous time and the propagation delay . The channel is expressed as a set of random complex gains at a given time and the propagation delay . The output of the channel can be expressed as the convolution of the complex channel impulse response and the input

If the complex channel gains are typically drawn from a complex Gaussian distribution, then at any given time , the absolute value of the impulse response is Rayleigh distributed (if the mean of the distribution or Rician distributed (if the mean of the distribution . These two scenarios model the presence or absence of a Line of Sight (LOS) path between the transmitter and the receiver.

Here, the values for the channel impulse response are samples of a random process that is defined with respect to time and the multipath delay τ. That is, for each combination of and τ, a randomly drawn value is assigned for the channel impulse response. As with any other random process, we can calculate the general autocorrelation function as

Given the generic autocorrelation function above, following assumptions can be made to restrict the channel model to the following specific set of categories

Wide Sense Stationary channel model
Uncorrelated Scattering channel model
Wide Sense Stationary Uncorrelated Scattering channel model

Wide Sense Stationary (WSS) channel model

In this channel model, the impulse response of the channel is considered Wide Sense Stationary (WSS) , that is the channel impulse response is independent of time . In other words, the autocorrelation function is independent of time instant $t$ and it depends on the difference between the time instants where and . The autocorrelation function for WSS channel model is expressed as

Uncorrelated Scattering (US) channel model

Here, the individual scattered components arriving at the receiver front end (at different propagation delays) are assumed to be uncorrelated. Thus the autocorrelation function can be expressed as

Wide Sense Stationary Uncorrelated Scattering (WSSUS) channel model

The WSSUS channel model combines the aspects of WSS and US channel model that are discussed above. Here, the channel is considered as Wide Sense Stationary and the scattering components arriving at the receiver are assumed to be uncorrelated. Combining both the worlds, the autocorrelation function is

Scattering function

The autocorrelation function of the WSSUS channel model can be represented in frequency domain by taking Fourier transform with respect one or both variables – difference in time and the propagation delay . Of the two forms, the Fourier transform on the variable gives specific insight to channel properties in terms of propagation delay and the Doppler Frequency simultaneously. The Fourier transform of the above two-dimensional autocorrelation function on the variable is called scattering function and is given by

Fourier transform of relative time is Doppler Frequency. Thus, the scattering function is a function of two variables – Dopper Frequency and the propagation delay . It gives the average output power of the channel as a function of Doppler Frequency and the propagation delay .

Figure 1: Relation between scattering function, power delay profile, Doppler power spectrum, spaced frequency
correlation function and spaced time correlation function (click to know more)

Two important relationships can be derived from the scattering function – Power Delay Profile (PDP) and Doppler Power Spectrum. Both of them affect the performance of a given wireless channel. Power Delay Profile is a function of propagation delay and the Doppler Power Spectrum is a function of Doppler Frequency.

Power Delay Profile gives the signal intensity received over a multipath channel as a function of propagation delays. It is obtained as the spatial average of the complex baseband channel impulse response as

Power Delay Profile can also be obtained from scattering function, by integrating it over the entire frequency range (removing the dependence on Doppler frequency).

Similarly, the Doppler Power Spectrum can be obtained by integrating the scattering function over the entire range of propagation delays.

Fourier Transform of Power Delay Profile and Inverse Fourier Transform of Doppler Power Spectrum:

Power Delay Profile is a function of time which can be transformed to frequency domain by taking Fourier Transform. Fourier Transform of Power Delay Profile is called spaced-frequency correlation function. Spaced-frequency correlation function describes the spreading of a signal in frequency domain. This gives rise to the importance channel parameter – Coherence Bandwidth.

Similarly, the Doppler Power Spectrum describes the output power of the multipath channel in frequency domain. The Doppler Power Spectrum when translated to time-domain by means of inverse Fourier transform is called spaced-time correlation function. Spaced-time correlation function describes the correlation between different scattered signals received at two different times as a function of the difference in the received time. It gives rise to the concept of Coherence Time.

Next Topic : Power Delay Profile

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

For further reading

[1] L. Bernadó, T. Zemen, F. Tufvesson, A. F. Molisch and C. F. Mecklenbräuker, “The (in-) validity of the WSSUS assumption in vehicular radio channels,” 2012 IEEE 23rd International Symposium on Personal, Indoor and Mobile Radio Communications – (PIMRC), Sydney, NSW, 2012, pp. 1757-1762, doi: 10.1109/PIMRC.2012.6362634.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

BLUE estimator

Why BLUE :

We have discussed Minimum Variance Unbiased Estimator (MVUE) in one of the previous articles. Following points should be considered when applying MVUE to an estimation problem

MVUE is the optimal estimator
Finding a MVUE requires full knowledge of PDF (Probability Density Function) of the underlying process.
Even if the PDF is known, finding an MVUE is not guaranteed.
If PDF is unknown, it is impossible find an MVUE using techniques like Cramer Rao Lower Bound (CRLB)
In practice, knowledge of PDF of the underlying process is actually unknown.

Considering all the points above, the best possible solution is to resort to finding a sub-optimal estimator. When we resort to find a sub-optimal estimator

We may not be sure how much performance we have lost – Since we will not able to find the MVUE estimator for bench marking (due to non-availability of underlying PDF of the process).
We can live with it, if the variance of the sub-optimal estimator is well with in specification limits

Common Approach for finding sub-optimal Estimator:

Restrict the estimator to be linear in data
Find the linear estimator that is unbiased and has minimum variance
This leads to Best Linear Unbiased Estimator (BLUE)
To find a BLUE estimator, full knowledge of PDF is not needed. Just the first two moments (mean and variance) of the PDF is sufficient for finding the BLUE

Definition of BLUE:

Consider a data set whose parameterized PDF depends on the unknown parameter . As the BLUE restricts the estimator to be linear in data, the estimate of the parameter can be written as linear combination of data samples with some weights

Here is a vector of constants whose value we seek to find in order to meet the design specifications. Thus, the entire estimation problem boils down to finding the vector of constants – . The above equation may lead to multiple solutions for the vector . However, we need to choose those set of values of , that provides estimates that are unbiased and has minimum variance.

Thus seeking the set of values for for finding a BLUE estimator that provides minimum variance, must satisfy the following two constraints

The estimator must be linear in data
Estimate must be unbiased

Constraint 1: Linearity Constraint:

Linearity constraint was already given above. Just repeated here for convenience.

Constraint 2: Constraint for unbiased estimates:

For the estimate to be considered unbiased, the expectation (mean) of the estimate must be equal to the true value of the estimate.

Thus,

Combining both the constraints (1) and (2) or (3),

Now, the million dollor question is : “When can we meet both the constraints ? “. We can meet both the constraints only when the observation is linear. That is is of the form , where is the unknown parameter that we wish to estimate.

Consider a data model, as shown below, where the observed samples are in linear form with respect to the parameter to be estimated.

Here , is zero mean process noise , whose PDF can take any form (Uniform, Gaussian, Colored etc., ). The mean of the above equation is given by

Substuiting (6) in (4) ,

Looking at the last set of equality,

The above equality can be satisfied only if

Given this condition is met, the next step is to minimize the variance of the estimate. Minimizing the variance of the estimate,

Finding BLUE:

As discussed above, in order to find a BLUE estimator for a given set of data, two constraints – linearity & unbiased estimates – must be satisfied and the variance of the estimate should be minimum. Thus the goal is to minimize the variance of which is subject to the constraint . This is a typical Lagrangian Multiplier↗ problem, which can be considered as minimizing the following equation with respect to (Remember !!! this is what we would like to find ).

Minimizing with respect to is equivalent to setting the first derivative of w.r.t to zero.

Substituting (12) in (9)

Finally, from (12) and (13), the co-effs of the BLUE estimator (vector of constants that weights the data samples) is given by

The BLUE estimate and the variance of the estimates are as follows

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Linear Models – Least Squares Estimator (LSE)

Key focus: Understand step by step, the least squares estimator for parameter estimation. Hands-on example to fit a curve using least squares estimation

Background:

The various estimation concepts/techniques like Maximum Likelihood Estimation (MLE), Minimum Variance Unbiased Estimation (MVUE), Best Linear Unbiased Estimator (BLUE) – all falling under the umbrella of classical estimation – require assumptions/knowledge on second order statistics (covariance) before the estimation technique can be applied. Linear estimators, discussed here, do not require any statistical model to begin with. It only requires a signal model in linear form.

Linear models are ubiquitously used in various fields for studying the relationship between two or more variables. Linear models include regression analysis models, ANalysis Of VAriance (ANOVA) models, variance component models etc. Here, one variable is considered as a dependent (response) variable which can be expressed as a linear combination of one or more independent (explanatory) variables.

Studying the dependence between variables is fundamental to linear models. For applying the concepts to real application, following procedure is required

Problem identification
Model selection
Statistical performance analysis
Criticism of the model based on statistical analysis
Conclusions and recommendations

Following text seeks to elaborate on linear models when applied to parameter estimation using Ordinary Least Squares (OLS).

Linear Regression Model

A regression model relates a dependent (response) variable y to a set of k independent explanatory variables {x₁, x₂ ,…, x_k} using a function. When the relationship is not exact, an error term e is introduced.

If the function f is not a linear function, the above model is referred as Non-Linear Regression Model. If f is linear, equation (1) is expressed as linear combination of independent variables x_k weighted by unknown vector parameters θ = {θ₁, θ₂,…, θ_k } that we wish to estimate.

Equation (2) is referred as Linear Regression model. When N such observations are made

where,
y_i – response variable
x_i – independent variables – known expressed as observed matrix X with rank k
θ_i – set of parameters to be estimated
e – disturbances/measurement errors – modeled as noise vector with PDF N(0, σ² I)

It is convenient to express all the variables in matrix form when N observations are made.

Denoting equation (3) using (4),

Except for X which is a matrix, all other variables are column/row vectors.

Ordinary Least Squares Estimation (OLS)

In OLS – all errors are considered equal as opposed to Weighted Least Squares where some errors are considered significant than others.

If is a k ⨉ 1 vector of estimates of θ, then the estimated model can be written as

Thus the error vector e can be computed from the observed data matrix y and the estimated as

Here, the errors are assumed to be following multivariate normal distribution with zero mean and standard deviation σ².

To determine the least squares estimator, we write the sum of squares of the residuals (as a function of ) as

The least squares estimator is obtained by minimizing . In order to get the estimate that gives the least square error, differentiate with respect to and equate to zero.

Thus, the least squared estimate of θ is given by

where the operator T denotes Hermitian Transpose (conjugate transpose).

Summary of computations

Step 1: Choice of variables. Choose the variable to be explained (y) and the explanatory variables { x₁, x₂ ,…, x_k } where x₁ is often considered a constant (optional) that always takes the value 1 – this is to incorporate a DC component in the model.
Step 2: Collect data. Collect n observations of y and for a set of known values of { x₁, x₂ ,…, x_k }. Example: { x₁, x₂ ,…, x_k } is the pilot data in OFDM using which we would like to estimate the channel impulse response θ and y is the received vector of samples. Store the observed data y in an – n⨉1 vector and the data on the explanatory variables in the n⨉k matrix X.
Step 3: Compute the estimates. Compute the least squares estimates by the formula

The superscript T indicates Hermitian Transpose (conjugate transpose) operation.

Key Points

We do not need a probabilistic assumption but only a deterministic signal model.
It has a broader range of applications.
Least squares is unbiased.
Estimating the disturbance variance (k variables to estimate and n observations are available).
To keep the variance low, the number of observations must be greater than the number of variables to estimate.
The observation matrix X should have maximum rank – this leads to independent rows and columns which always happens with real data. This will make sure (X^TX) is invertible.
Least Squares Estimator can be used in block processing mode with overlapping segments – similar to Welch’s method of PSD estimation.
Useful in time-frequency analysis.
Adaptive filters are utilized for non-stationary applications.

LSE applied to curve fitting

Matlab snippet for implementing Least Estimate to fit a curve is given below.

x = -5:.1:5; % set of x- values - known explanatory variables
y = 5.3 + 1.2* x; % Straight line without noise
e=randn(size(y));
y = y + e; % adding random noise to get observed variable - 
%Linear model - Y=Xa+e where a - parameters to be estimated

X = [ ones(length(x),1) x']; %first column treated aas all ones since x_1=1
y = y'; %column vector for proper dimension during multiplication
a = inv(X'*X)*X'*y  % Least Squares Estimator - equivalent code X\y
h=plot ( x , y , 'o'); %original data
hold on;
plot( x , a(1)+ a(2)*x , 'r-' ); %Fitted line
legend('observed samples',['y=' num2str(a(1)) '+' num2str(a(2)) 'x']) 
title('Least Squares Estimate for Curve Fitting');
xlabel('X values');
ylabel('Y values');

Simulation Results

*Figure 1: Least Squares Estimate for Curve Fitting*

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

AutoCorrelation (Correlogram) and persistence – Time series analysis

The agenda for the subsequent series of articles is to introduce the idea of autocorrelation, AutoCorrelation Function (ACF), Partial AutoCorrelation Function (PACF) , using ACF and PACF in system identification.

Introduction

Given time series data (stock market data, sunspot numbers over a period of years, signal samples received over a communication channel etc.,), successive values in the time series often correlate with each other. This series correlation is termed persistence or inertia and it leads to increased power in the lower frequencies of the frequency spectrum. Persistence can drastically reduces the degrees of freedom in time series modeling (AR, MA , ARMA models). In the test for statistical significance, presence of persistence complicates the test as it reduces the number of independent observations.

Autocorrelation function

Correlation of a time series with its own past and future values- is called autocorrelation. It is also referred as “lagged or series correlation”. Positive autocorrelation is an indication of a specific form of persistence, the tendency of a system to remain in the same state from one observation to the next (example: continuous runs of 0’s or 1’s). If a time series exhibits correlation, the future values of the samples probabilistic-ally depend on the current & past samples. Thus the existence of autocorrelation can be exploited in prediction as well as modeling time series. Autocorrelation can be accessed using the following tools

● Time series plot
● Lagged scatterplot
● AutoCorrelation Function (ACF)

Generating a sample time series

For the purpose of illustration, let’s begin by generating two time series data using Auto-Regressive AR(1) process. AR(1) process relates the current sample x[n] of the output of an LTI system, its immediate past sample x[n-1] and the white noise term w[n].

A generic AR(1) system is given by

Here and are the model parameters which we will tweak to generate different set of time series data and is a constant which will be set to zero . Thus the model can be equivalently written as

Let’s generate two time series data from the above model.

Model 1: a₀=0, a₁=1

The “Filter” function in Matlab will be utilized to generate the output process x[n]. The filter function, in its basic form – X=filter(B,A,W), takes three inputs. The vectors B and A denote the numerator and denominator co-efficients (model parameters here) of the transfer function of the LTI system in standard difference equation form, W is the white noise vector to the LTI filter and the output of filter is X.

The transfer function of model 1 is therefore,

Where B=1 and A=[1 -1] and the input W is a white noise – which can be generate using randn function. Therefore, the above model can be implemented with the command x=filter(1,[1 -1,randn(1000,1)) generate 1000 samples of x[n]

A=[1 -1];  %model co-effs
% generating using numerator/denominator form with noise
x1 = filter(1,A,randn(1000,1));
plot(x1,’b’);

Model 2: a₀=1, a₁=0.5

Transfer function of this model is

Where B=1 and A=[1 -0.5] and the input W is a white noise – which can be generated using randn function. Therefore, the above model can be implemented with the command x=filter(1,[1 -0.5], randn(1000,1)) to generate 1000 samples of x[n]

A=[1 -0.5];  %model co-effs
% generating using numerator/denominator form with noise
x2 = filter(1,A,randn(1000,1));
plot(x2,’r’);

Time-series plot of two models – where one model shows persistence and the other does not

In the plot above, the output from model 1 exhibits persistence or positive correlation – positive deviations from mean tend to be followed by positive deviations for some duration and the negative deviations from mean tend to be followed by negative deviations for sometime. When the positive deviations are followed by negative deviations or vice-versa, it is a characteristic of negative correlation. Positive correlations are strong indications of long runs of several consecutive observations above or below mean. Negative correlations indicate low incidence of such runs. The output of the model 2 always jumps around the mean value and there is no consistent departure from the mean – no persistence (no positive correlation). The interpretation of time series plots for clues on persistence is a subjective matter and is left for trained eyes. However, it can be considered as a preliminary analysis.

Persistence – an indication of non-stationarity:

For time series analysis, it is imperative to work with stationary process. Many of the formulated theorems in statistical signal processing assume a series to be stationary (atleast in weak sense). Processes whose Probability Density Functions do not change with time are termed stationary (sub classifications include strict sense stationarity (SSS), weak sense stationarity (WSS) etc.,). For analysis, the joint probability distribution must remain unchanged should there be any shift in the time series. Time series with persistence – changing mean with time – are non-stationary – therefore many theorems in signal processing will not apply as such.

Plotting the histogram of the two series (see next figure) , we can immediately identify that the data generated by model 1 is non-stationary – histogram varies between selected portion of the signal. Whereas, the histogram of the output from model 2 is pretty much same – therefore, this is a stationary signal and is suitable for further analysis.

Persistence – non-stationary and stationary signal

Lagged Scatter Plots

Autocorrelation trend can also be ascertained by lagged scatter plots. In lagged scatter plots, the samples of time series are plotted against one another with one lag at a time. A Strong positive autocorrelation will show of as a linear positive slope for the particular lag value. If the scatter plot is random, it indicates no-correlation for the particular lag.

figure;
x12 = x1(1:end-1);
x12 = x1(1:end-1);
x21 = x1(2:end);
x13 = x1(1:end-2);
x31 = x1(3:end);
x14 = x1(1:end-3);
x41 = x1(4:end);
x15 = x1(1:end-4);
x51 = x1(5:end);
subplot(2,2,1)
plot(x12,x21,'b*');
xlabel('X_1'); ylabel('X_2');
subplot(2,2,2)
plot(x13,x31,'b*');
xlabel('X_1'); ylabel('X_3');
subplot(2,2,3)
plot(x14,x41,'b*');
xlabel('X_1'); ylabel('X_4');
subplot(2,2,4)
plot(x15,x51,'b*');
xlabel('X_1'); ylabel('X_5');

figure;
x12 = x2(1:end-1);
x12 = x2(1:end-1);
x21 = x2(2:end);
x13 = x2(1:end-2);
x31 = x2(3:end);
x14 = x2(1:end-3);
x41 = x2(4:end);
x15 = x2(1:end-4);
x51 = x2(5:end);
subplot(2,2,1)
plot(x12,x21,'b*');
xlabel('X_1'); ylabel('X_2');
subplot(2,2,2)
plot(x13,x31,'b*');
xlabel('X_1'); ylabel('X_3');
subplot(2,2,3)
plot(x14,x41,'b*');
xlabel('X_1'); ylabel('X_4');
subplot(2,2,4)
plot(x15,x51,'b*');
xlabel('X_1'); ylabel('X_5');

The scatter plot of the model 1 for the first four lags indicate strong positive correlation at all the four lag values. The scatter plot of model 2 indicates a slightly positive correlation for lag=1 and no correlation for remaining lags. This trend can be clearing seen if we plot the Auto Correlation Function (ACF).

Auto Correlation Function (ACF) or Correlogram

ACF plot summarizes the correlation of a time series at various lags. It plots the correlation co-efficient of the series lagged by 1 delay at a time in the sample plot. Plotting the ACF for the output from both the models with the code below.

[x1c,lags] = xcorr(x1,100,'coeff');
%Plotting only positive lag values - autocorrelation is symmetric
stem(lags(101:end),x1c(101:end));
[x2c,lags] = xcorr(x2,100,'coeff');
stem(lags(101:end),x2c(101:end))

The ACF plot of model 1 indicates strong persistence across all the lags. The ACF plot of model 2 indicates significant correlation only at lag 1 (and lag 0 will obviously correlate fully) which concurs with the lagged scatter plots.

Auto Correlation Function or correlogram

Correlogram has very few significant spikes at very small lags and cuts off drastically/dies down quickly for stationary series. Thus model 2 produces stationary series, where as model 1 does not. Also, model 2 is suitable for further time series analysis.

Continue reading on constructing an autocorrelation matrix…

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Articles in this series

[1]	An Introduction to Estimation Theory
[2]	Bias of an Estimator
[3]	Minimum Variance Unbiased Estimators (MVUE)
[4]	Maximum Likelihood Estimation
[5]	Maximum Likelihood Decoding
[6]	Probability and Random Process
[7]	Likelihood Function and Maximum Likelihood Estimation (MLE)
[8]	Score, Fisher Information and Estimator Sensitivity
[9]	Introduction to Cramer Rao Lower Bound (CRLB)
[10]	Cramer Rao Lower Bound for Scalar Parameter Estimation
[11]	Applying Cramer Rao Lower Bound (CRLB) to find a Minimum Variance Unbiased Estimator (MVUE)
[12]	Efficient Estimators and CRLB
[13]	Cramer Rao Lower Bound for Phase Estimation
[14]	Normalized CRLB - an alternate form of CRLB and its relation to estimator sensitivity
[15]	Cramer Rao Lower Bound (CRLB) for Vector Parameter Estimation
[16]	The Mean Square Error – Why do we use it for estimation problems
[17]	How to estimate unknown parameters using Ordinary Least Squares (OLS)
[18]	Essential Preliminary Matrix Algebra for Signal Processing
[19]	Why Cholesky Decomposition ? A sample case:
[20]	Tests for Positive Definiteness of a Matrix
[21]	Solving a Triangular Matrix using Forward & Backward Substitution
[22]	Cholesky Factorization - Matlab and Python
[23]	LTI system models for random signals – AR, MA and ARMA models
[24]	Comparing AR and ARMA model - minimization of squared error
[25]	Yule Walker Estimation
[26]	AutoCorrelation (Correlogram) and persistence – Time series analysis
[27]	Linear Models - Least Squares Estimator (LSE)
[28]	Best Linear Unbiased Estimator (BLUE)

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Yule Walker Estimation and simulation in Matlab

If a time series data is assumed to be following an Auto-Regressive (AR(N)) model of given form,

the natural tendency is to estimate the model parameters a₁,a₂,…,a_N.

Least squares method can be applied here to estimate the model parameters but the computations become cumbersome as the order N increases.

Fortunately, the AR model co-efficients can be solved for using Yule Walker Equations.

Yule-Walker Equations

Yule Walker equations relate auto-regressive model parameters to auto-covariance r_xx[k] of random process x[n]. It can be applied to find the AR(N) model parameters from a given time-series data x[n] that indeed can be assumed to be AR process (visually examining auto-correlation function (ACF) and partial auto-correlation function (PACF) may give clue on whether the data can be assumed as an AR or MA process with appropriate model order N).

Steps to find model parameters using Yule-Walker Equations:

Given x[n], estimate auto-covariance of the process r_xx[k]
Solve Yule Walker Equations to find a_k (k=0,1,…,N) and σ² from

To solve Yule-Walker equations, one has to know the model order. Although the Yule-Walker equations can be used to find the model parameters, it cannot give any insight into the model order N directly. Several methods exists to ascertain the model order.

ACF and PACF Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) may give insight into the model order of the underlying process
Akaike Information Criterion (AIC)^[1] It associates the number of model parameters and the goodness of fit. It also associates a penalty factor to avoid over-fitting (increasing model order always improves the fit – the fitting process has to be stopped somewhere)
Bayesian Information Criterion (BIC)^[2] Very similar to AIC
Cross Validation^[3] The given time series is broken into subsets. One subset of data is taken at a time and model parameters estimated. The estimated model parameters are cross validated with data from the remaining subsets.

Deriving Yule-Walker Equations

A generic AR model

can be written in compact form as

Note that the scaling factor x[n] term is a₀=1

Multiplying (2) by x[n-l]

One can readily identify the auto-correlation and cross-correlation terms as

The next step is to compute the identified cross-correlation r_wx[l] term and relate it to the auto-correlation term r_xx[l-k]

The term x[n-l] can also be obtained from equation (1) as

Note that data and noise are always uncorrelated, therefore: x[n-k-l]w[n]=0. Also the auto-correlation of noise is zero at all lags except at lag 0 where its value is equal to σ² (remember the flat Power Spectral Density of white noise and its auto-correlation). These two properties are used in the following steps. Restricting the lags only to positive values and zero ,

Substituting (6) in (4),

Here there are two cases to be solved – l>0 and l=0
For l>0 case, equation (7) becomes,

Clue: notice the lower limit of the summation which changed from k=0 to k=1 .

Now, equation (8) can be written in matrix form

This is the Yule-Walker Equations which comprises of a set of N linear equations and N unknown parameters.

Representing equation (9) in a compact form

The solutions can be solved by

Once we solve for , equivalently a_k, – the model parameters, the noise variance σ² can be found by applying the estimated values of a_k in equation (7) by setting l=0

Matlab’s “aryule” efficiently solves the “Yule-Walker” equations using “Levinson Algorithm”^[4][5]

Simulation:

Let’s generate an AR(3) process and pretend that we do not anything about the model parameters. We will take this as input data to Yule-Walker and check if it can estimate the model parameters properly

Generating the data from the AR(3) process given above

a=[1  0.9   0.6 0.5]; %model parameters
x=filter(1,a,randn(1000,1)); %data generation with gaussian noise - 1000 samples
plot(x);
title('Sample data from AR(3) process');
xlabel('Sample index');
ylabel('Value');

*Figure 1: Simulated data from AR(3) process*

Plot the periodogram (PSD plot) for reference. The PSD plots from the estimated model will be checked against this reference plot later.

figure(2);
periodogram(x); hold on; %plot the original frequency response of the data

Running the simulation for three model orders and checking which model order suits the best.

N=[2,3,4];

for i=1:3,
    [a,p] = aryule(x,N(i))
    [H,w]=freqz(sqrt(p),a);
    hp = plot(w/pi,20*log10(2*abs(H)/(2*pi))); %plotting in log scale
end
xlabel('Normalized frequency (\times \pi rad/sample)');
ylabel('One-sided PSD (dB/rad/sample)');
legend('PSD estimate of x','PSD of model output N=2','PSD of model output N=3','PSD of model output N=4');

*Figure 2: Power Spectral density of modeled data*

The estimated model parameters and the noise variances computed by the Yule-Walker system are given below. It can be ascertained that the estimated parameters are approximately same as that of what is expected. See how the error decreases as the model order N increases. The optimum model order in this case is N=3 since the error has not changed significantly when increasing the order and also the model parameter a₄ of the AR(4) process is not-significantly different from zero.

*Figure 3: Estimated model parameters and prediction errors for the given model model order*

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References

[1] Akaike H. (1998) Information Theory and an Extension of the Maximum Likelihood Principle. In: Parzen E., Tanabe K., Kitagawa G. (eds) Selected Papers of Hirotugu Akaike. Springer Series in Statistics (Perspectives in Statistics). Springer, New York, NY. https://doi.org/10.1007/978-1-4612-1694-0_15
[2] McQuarrie, A. D. R., and Tsai, C.-L., 1998. Regression and Time Series Model Selection. World Scientific.↗
[3] Arlot, Sylvain, and Alain Celisse. “A survey of cross-validation procedures for model selection.” Statistics surveys 4 (2010): 40-79.
[4] J. Durbin.The fitting of time series in models.Review of the International Statistical Institute, 28:233-243, 1960
[5] G. H. Golub and C. F. Van Loan, “Matrix Computations”, The Johns Hopkins University Press, third edition, 1996.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Cookie	Duration	Description
cookielawinfo-checbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checkbox-necessary	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookies is used to store the user consent for the cookies in the category "Necessary".
cookielawinfo-checkbox-performance	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Performance".
viewed_cookie_policy	11 months	The cookie is set by the GDPR Cookie Consent plugin and is used to store whether or not user has consented to the use of cookies. It does not store any personal data.

Introduction

Sine Wave

Representing in Frequency Domain

Different representations of FFT:

1. Plotting raw values of DFT:

2. FFT plot – plotting raw values against Normalized Frequency axis:

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

4. FFT plot – Absolute frequency on the x-axis Vs Magnitude on Y-axis:

5. Power Spectrum – Absolute frequency on the x-axis Vs Power on Y-axis:

6. Power Spectrum – One-Sided frequencies

For further reading

Topics in this chapter

Books by the author

Correlation Matrix

Sample problem:

Generating Correlated random number using Cholesky Decomposition:

Steps to follow:

Python code

Matlab code

Further reading

Topics in this chapter

Books by the author

Step 1: Generate two uncorrelated Gaussian distributed random sequences <img decoding="async" src="https://s0.wp.com/latex.php?latex=X_1%2C+X_2&amp;bg=ffffff&amp;fg=000&amp;s=0&amp;c=20201002" alt="X_1, X_2" width="46" height="16">

Step 2: Generate correlated random sequence z

Further reading

Topics in this chapter

Books by the author

Manipulating audio files in Matlab

Example: DownSampling

Reference

Books by the author

Power delay profile (PDP)

RMS delay spread and mean delay

Frequency selective and non-selective channels

References:

Topics in this chapter

Books by the author

Introduction

Complex baseband mutipath channel model:

Wide Sense Stationary (WSS) channel model

Uncorrelated Scattering (US) channel model

Wide Sense Stationary Uncorrelated Scattering (WSSUS) channel model

Scattering function

Fourier Transform of Power Delay Profile and Inverse Fourier Transform of Doppler Power Spectrum:

For further reading

Books by the author

Why BLUE :

Common Approach for finding sub-optimal Estimator:

Definition of BLUE:

Constraint 1: Linearity Constraint:

Constraint 2: Constraint for unbiased estimates:

Finding BLUE:

Similar topics

Books by the author

Background:

Linear Regression Model

Ordinary Least Squares Estimation (OLS)

Summary of computations

Key Points

LSE applied to curve fitting

Simulation Results

Related topics:

Books by the author

Introduction

Autocorrelation function

Generating a sample time series

Model 1: a0=0, a1=1

Model 2: a0=1, a1=0.5

Persistence – an indication of non-stationarity:

Lagged Scatter Plots

Auto Correlation Function (ACF) or Correlogram

Articles in this series

Books by the author

Yule-Walker Equations

Deriving Yule-Walker Equations

Simulation:

References

Related topics

Books by the author

Step 1: Generate two uncorrelated Gaussian distributed random sequences

Model 1: a₀=0, a₁=1

Model 2: a₀=1, a₁=0.5