# Connexions

You are here: Home » Content » Purdue Digital Signal Processing Labs (ECE 438) » Lab 7c - Power Spectrum Estimation

### Lenses

What is a lens?

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

#### Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
• NSF Partnership

This collection is included inLens: NSF Partnership in Signal Processing
By: Sidney Burrus

Click the "NSF Partnership" link to see all content affiliated with them.

Click the tag icon to display tags associated with this content.

• Featured Content

This collection is included inLens: Connexions Featured Content
By: Connexions

Click the "Featured Content" link to see all content affiliated with them.

Click the tag icon to display tags associated with this content.

#### Also in these lenses

• UniqU content

This collection is included inLens: UniqU's lens
By: UniqU, LLC

Click the "UniqU content" link to see all content selected in this lens.

• Lens for Engineering

This module and collection are included inLens: Lens for Engineering
By: Sidney Burrus

Click the "Lens for Engineering" link to see all content selected in this lens.

### Recently Viewed

This feature requires Javascript to be enabled.

### Tags

(What is a tag?)

These tags come from the endorsement, affiliation, and other lenses that include this content.

Inside Collection (Course):

Course by: Charles A. Bouman. E-mail the author

# Lab 7c - Power Spectrum Estimation

Module by: Charles A. Bouman. E-mail the author

Questions or comments concerning this laboratory should be directed to Prof. Charles A. Bouman, School of Electrical and Computer Engineering, Purdue University, West Lafayette IN 47907; (765) 494-0340; bouman@ecn.purdue.edu

## Introduction

In the first and second weeks of this experiment, we introduced methods of statistically characterizing random processes. The sample autocorrelation and cross correlation are examples of “time domain” characterizations of random signals. In many applications, it can also be useful to get the frequency domain characteristics of a random process. Examples include detection of sinusoidal signals in noise, speech recognition and coding, and range estimation in radar systems.

In this week, we will introduce methods to estimate the power spectrum of a random signal given a finite number of observations. We will examine the effectiveness of the periodogram for spectrum estimation, and introduce the spectrogram to characterize a nonstationary random processes.

## Power Spectrum Estimation

In this section, you will estimate the power spectrum of a stationary discrete random process. The power spectrum is defined as

S x x ( ω ) = lim N E 1 N n = 0 N - 1 x ( n ) e - j ω n 2 S x x ( ω ) = lim N E 1 N n = 0 N - 1 x ( n ) e - j ω n 2
(1)

There are 4 steps for calculating a power spectrum:

1. Select a window of length NN and generate a finite sequence x(0),x(1),,x(N-1)x(0),x(1),,x(N-1).
2. Calculate the DTFT of the windowed sequence x(n)x(n), (n=0,1,,N-1)(n=0,1,,N-1), square the magnitude of the DTFT and divide it by the length of the sequence.
3. Take the expectation with respect to xx.
4. Let the length of the window go to infinity.

In real applications, we can only approximate the power spectrum. Two methods are introduced in this section. They are the periodogram and the averaged periodogram.

## Periodogram

Click here for help on the fft command and the random command.

The periodogram is a simple and common method for estimating a power spectrum. Given a finite duration discrete random sequence x(n)x(n), (n=0,1,,N-1)(n=0,1,,N-1), the periodogram is defined as

P x x ( ω ) = 1 N X ( ω ) 2 = 1 N n = 0 N - 1 x ( n ) e - j ω n 2 P x x ( ω ) = 1 N X ( ω ) 2 = 1 N n = 0 N - 1 x ( n ) e - j ω n 2
(2)

where X(ω)X(ω) is the Discrete Time Fourier Transform (DTFT) of x(n)x(n).

The periodogram Pxx(ω)Pxx(ω) can be computed using the Discrete Fourier Transformation (DFT), which in turn can be efficiently computed by the Fast Fourier Transformation (FFT). If x(n)x(n) is of length NN, you can compute an NN-point DFT.

P x x ( ω k ) = 1 N n = 0 N - 1 x ( n ) e - j ω k n 2 P x x ( ω k ) = 1 N n = 0 N - 1 x ( n ) e - j ω k n 2
(3)
ω k = 2 π k N , k = 0 , 1 , , N - 1 ω k = 2 π k N , k = 0 , 1 , , N - 1
(4)

Using Equation 3 and Equation 4, write a Matlab function called Pgram to calculate the periodogram. The syntax for this function should be

[P,w] = Pgram(x)

where xx is a discrete random sequence of length NN. The outputs of this command are PP, the samples of the periodogram, and ww, the corresponding frequencies of the samples. Both PP and ww should be vectors of length NN.

Now, let xx be a Gaussian (Normal) random variable with mean 0 and variance 1. Use Matlab function random or randn to generate 1024 i.i.d. samples of xx, and denote them as x0x0, x1x1, ..., x1023x1023. Then filter the samples of xx with the filter which obeys the following difference equation

y ( n ) = 0 . 9 y ( n - 1 ) + 0 . 3 x ( n ) + 0 . 24 x ( n - 1 ) . y ( n ) = 0 . 9 y ( n - 1 ) + 0 . 3 x ( n ) + 0 . 24 x ( n - 1 ) .
(5)

Denote the output of the filter as y0,y1,,y1023y0,y1,,y1023.

### Note:

The samples of x(n)x(n) and y(n)y(n) will be used in the following sections. Please save them.

Use your Pgram function to estimate the power spectrum of y(n)y(n), Pyy(ωk)Pyy(ωk). Plot Pyy(ωk)Pyy(ωk) vs. ωkωk.

Next, estimate the power spectrum of yy, Pyy(ωk)Pyy(ωk) using 1/41/4 of the samples of yy. Do this only using samples y0y0, y1y1, , y255y255. Plot Pyy(ωk)Pyy(ωk) vs. ωkωk.

### INLAB REPORT

1. Hand in your labeled plots and your Pgram code.
2. Compare the two plots. The first plot uses 4 times as many samples as the second one. Is the first one a better estimation than the second one? Does the first give you a smoother estimation?
3. Judging from the results, when the number of samples of a discrete random variable becomes larger, will the estimated power spectrum be smoother?

## Averaged Periodogram

The periodogram is a simple method, but it does not yield very good results. To get a better estimate of the power spectrum, we will introduce Bartlett's method, also known as the averaged periodogram. This method has three steps. Suppose we have a length-N sequence x(n)x(n).

1. Subdivide x(n)x(n) into KK nonoverlapping segments of length MM. Denote the ithith segment as xi(n)xi(n).
xi(n)=x(n+iM),i=0,1,,K-1,n=0,1,,M-1xi(n)=x(n+iM),i=0,1,,K-1,n=0,1,,M-1
(6)
2. For each segment xi(n)xi(n), compute its periodogram
Pxixi(i)(ωk)=1Mn=0M-1xi(n)e-jωkn2Pxixi(i)(ωk)=1Mn=0M-1xi(n)e-jωkn2
(7)
ωk=2πkMωk=2πkM
(8)
where k=0,1,,M-1k=0,1,,M-1 and i=0,1,,K-1i=0,1,,K-1.
3. Average the periodograms over all KK segments to obtain the averaged periodogram, PxxA(k)PxxA(k),
PxxA(ωk)=1Ki=0K-1Pxixi(i)(ωk)PxxA(ωk)=1Ki=0K-1Pxixi(i)(ωk)
(9)
ωk=2πkMωk=2πkM
(10)
where k=0,1,,M-1k=0,1,,M-1.

Write a Matlab function called AvPgram to calculate the averaged periodogram, using the above steps. The syntax for this function should be

[P,w] = AvPgram(x, K)

where xx is a discrete random sequence of length NN and the KK is the number of nonoverlapping segments. The outputs of this command are PP, the samples of the averaged periodogram, and ww, the corresponding frequencies of the samples. Both PP and ww should be vectors of length M where N=KM. You may use your Pgram function.

### Hint:

The command A=reshape(x,M,K) will orient length MM segments of the vector xx into KK columns of the matrix AA.

Use your Matlab function AvPgram to estimate the power spectrum of y(n)y(n) which was generated in the previous section. Use all 1024 samples of y(n)y(n), and let K=16K=16. Plot PP vs. ww.

### INLAB REPORT

1. Submit your plot and your AvPgram code.
2. Compare the power spectrum that you estimated using the averaged periodogram with the one you calculated in the previous section using the standard periodogram. What differences do you observe? Which do you prefer?

## Power Spectrum and LTI Systems

Consider a linear time-invariant system with frequency response H(ejω)H(ejω), where Sxx(ω)Sxx(ω) is the power spectrum of the input signal, and Syy(ω)Syy(ω) is the power spectrum of the output signal. It can be shown that these quantities are related by

S y y ( ω ) = H ( e j ω ) 2 S x x ( ω ) . S y y ( ω ) = H ( e j ω ) 2 S x x ( ω ) .
(11)

In the "Periodogram" section, the sequence y(n)y(n) was generated by filtering an i.i.d. Gaussian (mean=0, variance=1) sequence x(n)x(n), using the filter in Equation 5. By hand, calculate the power spectrum Sxx(ω)Sxx(ω) of x(n)x(n), the frequency response of the filter, H(ejω)H(ejω), and the power spectrum Syy(ω)Syy(ω) of y(n)y(n).

### Hint:

In computing Sxx(ω)Sxx(ω), use the fact that a b 2 = a b * . a b 2 = a b * .
Plot Syy(ω)Syy(ω) vs. ωω, and compare it with the plots from the "Periodogram" and "Averaged Periodogram" sections. What do you observe?

### INLAB REPORT

2. Submit your analytical calculations for Sxx(ω)Sxx(ω), H(ejω)H(ejω), Syy(ω)Syy(ω).
3. Compare the theoretical power spectrum of Syy(ω)Syy(ω) with the power spectrum estimated using the standard periodogram. What can you conclude?
4. Compare the theoretical power spectrum of Syy(ω)Syy(ω) with the power spectrum estimated using the averaged periodogram. What can you conclude?

## Estimate the Power Spectrum of a Speech Signal

Download the file speech.au for this section. For help on the following Matlab topics select the corresponding link: how to load and play audio signals and specgram function.

The methods used in the last two sections can only be applied to stationary random processes. However, most signals in nature are not stationary. For a nonstationary random process, one way to analyze it is to subdivide the signal into segments (which may be overlapping) and treat each segment as a stationary process. Then we can calculate the power spectrum of each segment. This yields what we call a spectrogram.

While it is debatable whether or not a speech signal is actually random, in many applications it is necessary to model it as being so. In this section, you are going to use the Matlab command specgram to calculate the spectrogram of a speech signal. Read the help for the specgram function. Find out what the command does, and how to calculate and draw a spectrogram.

Draw the spectrogram of the speech signal in speech.au. When using the specgram command with no output arguments, the absolute value of the spectrogram will be plotted. Therefore you can use

speech=auread('speech.au');

to read the speech signal and use

specgram(speech);

to draw the spectrogram.

### INLAB REPORT

1. Hand in your spectrogram plot.
2. Describe the information that the spectrogram is giving you.

## Content actions

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

#### Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

#### Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks