Description

ECE 310 Project 1 Sampling Rate Conversion

Overview
In discrete‐time signal processing, it is often necessary to convert a signal from one sampling rate to
another. A common example is the conversion from the sampling rate of a compact disk (CD); signal
(44.1 KHz) to that of a Digital Audio Tape (DAT) signal (48 KHz). Another example is the aud io standard
in High‐Definition Television (HDTV) transmission, where at least three sampling rates are supported
(32, 44.1, and 48 KHz). Although in principle we may convert the signal back to analog form and
resample at the desired rate, it is usually preferable to perform t he entire conversion digitally. This is
clue to many considerations including the fact that conversion to analog form often introduces noise in
the signal, and that digital signal processing can be much more cost‐effective and flexible.
This MATLAB project asks you to perform a sampling rate conversion on segments of audio signals. The
input audio signals are quantized to 8 bits and sampled with a sampling frequency of 11,025 Hz. You are
required to convert the signal to a sampling rate of 24,000 Hz in a computationally efficient manner.
Although one conceptual way of realizing this sampling rate conversion process is to upsample the
signal, lowpass filter, and downsample it, a more clever implementation can lead to an implementation
that is many times more efficient. To do this, you can exploit various aspects of class to optimize the
system, such as multistage filter implementation, filter design, and polyphase implementation. By the
end of the project, you hopefully will have a much better understanding of both the theoretical aspect
of the system as well as various issues in implementing a practical DSP system at a software level.
Project Goal
A sampling rate converter which produces an output signal with a sampling rate which is
M
L times the
original sampling rate can be specified as shown in Figure 1.
Figure 1: Sampling rate conversion system.
For an ideal sampling rate converter, the lowpass filter in Figure 1 is an ideal lowpass filter with cutoff
frequency ⎟
⎠
⎞ ⎜
⎝
⎛ = LM c
ππ
ω ,min The goal of this project is to design an efficient DSP algorithm that
implements the system in Figure 1 subject to the following constraints:
• The system performs the correct sampling rate conversion from 11,025 Hz to 24,000 Hz. In
particular, do not assume that 11,025 Hz ≈ 11,000 Hz.
x[ ] n y n[ ] Lowpass
filter
↑ L ↓ M
f1997
If the system in Figure 1 is an ideal sampling rate converter , when the input x[n] is a unit impulse
δ n],[ the output y[n] has a Fourier transform ( ) jω eY which corresponds to an ideal lowpass filter with
a cutoff frequency at ωc = (11,025/ 24,000)π. For an equivalent system which you are to implement,
when the input x[n] is a unit impulseδ n][ , the output y[n] must have a Fourier transform ( ) jω eY which is
an approximation of a lowpass filter, and meets the specifications shown below in Table 1.
Passband Cutoff (ωp)
000,24
025,11 π
Passband Ripple ±0.1 dB or less
Stopband Frequency (ωs) 1.2ωp
Stopband Attenuation 70 dB or more
Phase Constraints |max grpdelay ‐ min grpdelay| ≤ 720 in the passband
Table 1
If your sampling rate conversion system functions properly, you should meet the specifications in Table
1, and when you play the output audio signal at 24,000 Hz, it should sound the same as the input audio
signal played at 11,025 Hz.
You should get a rough estimate of the efficiency of your design by determining how much computation
was required to perform the sampling rate conversion on the Wagner.wav signal. The method which
you will use to count the number of operations will be described shortly.
You are to write a MATLAB function srconvert such that the command srconvert(in) takes the
input signal in with an associated sampling rate of 11,025 Hz, and returns an output signal at a
sampling rate of 24,000 Hz. (See the section on writing MATLAB functions.) Once your function
srconvert.m is finalized, run the command y=srconvert([1 zeros(1,3000)]);. This will
produce a vector y which contains the response of your system to a unit impulse. Then call
verify(y) to verify that your design meets the design specification.
The Files
The project zip file contains the audio files and MATLAB functions which you will need for this project.
To access the zip file, click on the link below where you found this file on the web site.
In testing your system, you may find it helpful to first use a unit impulse as input to your system and
then later try using real audio signals as inputs and listening to the outputs. To load a sound file into a
vector x, type x=wavread(′filename.wav′);. To write a sound file for the MATLAB vector x into
the current directory, type wavwrite(x,sampfreq, ′filename.wav′). You may test your
system using any of the audio signals, although you should use the Wagner.wav signal to benchmark
your system performance.
f1997
MATLAB Utilities
You may find the following MATLAB functions useful for your design.
examlpf(h, wp, vs) ‐ allows you to examine the passband ripple, the group delay, as well as the
stopband attenuation of the lowpass filter whose impulse response is h. It generates three
simultaneous plots. The first one zooms in on the passband with cutoff frequency wp. The
second plot measures the group delay in the passband, and the third plot is simply the
magnitude response over the entire frequency range (-π,π). The passband (wp) and stopband
(ws) cutoff frequencies are normalized by π/2. That is, a cutoff frequency at π/2 should be
entered as 0.5.
poly1(h, M) ‐ returns a matrix E whose i
th row corresponds to the i
th polyphase components of the
FIR filter h, obtained via type I decomposition with downsampling factor M.
fftfilt(h , x) ‐ convolves the signal x with the filter h using the FFTs. Beware that this command
might not always leads to a more efficient implementation than conv(h ,x), depending on
the length of the signals.
upsample(h ,L) ‐ returns an upsampled version of h by a factor L, without any interpolation.
downsample(h ,M) ‐ downsamples h by M.
Computation Counts
MATLAB does not automatically record the number of floating point operations (flops) in your algorithm, so
to measure your computational efficiency, you’ll have to determine the required number of multiplies and
adds yourself. You can make this process somewhat automated by inserting a counter to be updated each
time the filter function is called, but you should also do a hand calculation to double check.
Do not include in the flops count the computation involved in designing any filters you need.
Example MATLAB Function
To give you an example of how to write a MATLAB function, here is the code from the file poly1.m.
function E=poly1(h,M)
%
% Performs type I polyphase decomposition of h in M components.
% The ith row of E corresponds to the ith polyphase component.
% Assumes that the first point of h is index 0.
%
h = [h zeros(1, ceil(length(h)/M)*M-length(h))];
E = reshape(h, M, length(h)/M);
f1997
Tips
The old homework problem in the Appendix might be useful as a guideline to efficient implementation of an
interpolation filter. You may also find useful page 85 of the Vaidyanathan paper cited below.
P. P. Vaidyanathan, “Multirate digital filters, filter banks, polyphase networks, and applications: A
tutorial,” Proceedings of the IEEE, vol. 78, pp. 56‐92, Jan. 1990.
https://www.systems.caltech.edu/dsp/ppv/papers/ProcIEEEmultirateTUTExtra.pdf (Available 7/08)
Appendix – An Old Homework Problem:
x[ ] n ( ) c Let represent a signal which is obtained by sampling a continuous time audio signal x t
/ 44.1 kHz T
using an ideal C/D
converter shown in Figure A1. Assume that the sampling rate is 1 = .
Figure A1
We wish to design a 4x (4 times) oversampling digital FIR interpolation filter. One way to do this is to use the
single‐stage design of System 1 shown in Figure A2.
Figure A2: System 1
The filter h n[ ] with Fourier transform )( j H e ω is designed using the window method with a Kaiser window.
Suppose the filter h n[ ] is designed to meet Specification 1:
Specification 1:
1 1 ( )1 , j
j
H e ω
ω
2 H e( ) ,
δ δ ωω
δ ω ω π
− ≤ ≤+ ≤
≤ ≤ ≤
where 1 ω = 0.23π 2 , ω = 0.27π
4 , and δ 10− =
h n[ ]
.
(a) Using a Kaiser window, estimate the length of the filter which meets Specification 1.
(b) Using the filter length estimate for in part(a), estimate the number of multiplications per second
required in System 1 if the system is implemented in polyphase form using four polyphase components.
h n[ ]
An alternate way to design the 4x oversampling interpolation filter is to use the two‐stage design of System 2
shown in Figure A3.
x[ ] n
↑ 4
1 y n[ ]
h n[ ]
( ) x[ ] n c x t C/D
T
f1997
Stage 1 Stage 2
2 y n[ ]
f1997
)
Figure A3: System 2
The filters and with frequency responses 1 h n[ ] 2 h n[ ] 1( j H e ω and 2 ( j H e ) ω , respectively, are designed using
the window method with Kaiser windows. Suppose that these filters are designed to meet Specification 2.
Specification 2:
1 1
1
2 1
2 1
2
1
1
( )1 , 2
() , 2
( )1 ,
() ,
j
j
j
H e
H e
H e
H e
ω
ω
ω
δ δω jω ω
δ ω ωπ
δ δω ω
δ πω ω π
′ ′ ≤ ≤+ ≤
≤ ≤ ′
′ ′ ≤ ≤+ ≤
≤ −≤ ′
−
≤
≤
−
where δ δ ′ = +− 1 1,ω1 , ω2 , and δ are given previously.
(c) Show that System 2 is equivalent to System 1 with
1 2 He H e H e ( ) ) ()( j jj ω ω ω
h n[ ]
j
=
2 (d) Show that if the filters and h [ ] n are designed using Kaiser windows to meet Specification 2, then
the resulting equivalent system in System 1 will have a filter
1
H e )( ω which meets Specification 1. (You
may assume that the filters designed have magnitude responses which are monotonically decreasing in
the transition band from the passband to the stopband).
(e) Using Kaiser windows, estimate the lengths of the filters and h which meet Specification 2. 1 h n[ ] [ ] n
h n[ ] h n[ ]
2
(f) Using the filter length estimates for and in part (e), estimate the number of multiplications
per second required in System 2 if each stage is separately implemented in polyphase form using two
polyphase components. What are the computational savings in multiplications, if any, over System 1?
1 2
x[ ] n
↑ 2 1 h n[ ] 2 ↑ 2 h n[ ]

ECE 310 Project 2 Testing a “Folk Theorem”

One area in which discrete-time signal processing has enjoyed widespread use is in the field of speech processing. In this project we consider various issues relating to applying digital filtering to speech. Specifically, in this project, we consider the effect of all-pass filtering on speech.

Historically, DSP courses have put considerable emphasis on design techniques for both IIR and FIR filters, as discussed in Chapter 7 of the course text. A variety of filter design algorithms are now implemented in common software packages such as fdatool in Matlab. This makes it unnecessary for most practitioners to learn the details of the design algorithms; however, ceding the design function to a software package makes it more important to understand the properties of different types of optimal filters, the meanings of the design parameters, and some of the trade-offs between filter classes.

In the context of IIR filter design, we suggest that you read Sections 7.0 and 7.1 including Example 7.1.

Project Assignment

Your writeup should have a summary of how you went about this project (anywhere from one to five pages, not to exceed five pages). In addition to the five page maximum, you can include any appropriate Matlab plots. Do NOT turn in any Matlab code. The total number of pages including description and plots should not exceed ten pages.

The writeup will receive a grade which reflects our assessment of your level of understanding of the various issues involved and the effort that you put into it.

Playing Sound in Matlab

Throughout this project it will be necessary to play audio files which you will have manipulated in Matlab. The functions sound and soundsc (sound scaled) work well for doing this on most hardware platforms.

Introduction

A common “folk theorem” states that the ear is insensitive to phase, i.e. that for audio, phase distortion is inaudible. If that is correct, then processing audio with an all-pass filter should not result in perceived distortion. This project tests this conjecture. The all-pass filter that we consider is of the form

.

For this part of the project, the file projIA.mat will need to be first loaded into Matlab. To access the project zip file, click on the link below where you found this file on the web site. After downloading the zip file, copy its contents into your Matlab working directory and type load projIA. The b_k and a_k coefficients above are stored in the variables b and a respectively; b(k + 1)= b_k and a(k + 1)= a_k.

• Implement the all-pass filter for N = 1 and show plots of the impulse response (the first 100 samples), the magnitude of the frequency response, and the group delay.
• Plot the pole-zero diagram. You may find the functions roots and zplane helpful in doing so.
• Process the speech file, stored as the variable speech (sampled at rate f_s = 11025 Hz), using the command filter to run the filter implemented in (a). Listen to the result and comment specifically on whether there is audible distortion.
• Implement the filter for N = 50 and again show plots of the impulse response (the first 5000 samples), the magnitude of the frequency response, and the group delay. Note: This is the hard part of the exercise, you can’t skip it.
• Repeat (c) with N = 50. How would you subjectively describe the distortion and what aspect of the filter is its primary cause?

ECE 310 Project 3 Enhancing Speech by Removing Noise

One area in which discrete-time signal processing has enjoyed widespread use is in the field of speech processing. In this project, we consider IIR and FIR filter design in the context of enhancing speech corrupted by additive noise. While this relies somewhat on material Chapters 6 and 7, some aspects can be started earlier.

This project is an excellent introduction to the issues discussed in Hardware Implementation Considerations for Removing Noise.

Historically, DSP courses have put considerable emphasis on design techniques for both IIR and FIR filters, as discussed in Chapter 7 of the course text. A variety of filter design algorithms are now implemented in common software packages such as Matlab. This makes it unnecessary for most practitioners to learn the details of the design algorithms; however, ceding the design function to a software package makes it more important to understand the properties of different types of optimal filters, the meanings of the design parameters, and some of the trade-offs between filter classes.

In the context of IIR filter design, we suggest that you read Sections 7.0 and 7.1 including Example 7.1. Also, read the examples in Section 7.3.1, but focus on what Butterworth and Chebyshev filters are rather than on the details of the bilinear transformation. For the discussion of FIR filter design, the suggested reading is Sections 7.5, 7.6, and 7.7 through 7.7.1. We will not be going into the details of the issues raised in Section 7.7.1, but we would recommend at least reading that section to get a sense of what some of those issues are.

Project Assignment

Your writeup should have a summary of how you went about it (anywhere from one to five pages, not to exceed five pages). In addition to the five page maximum, you can include any appropriate Matlab plots. Do NOT turn in any Matlab code. The total number of pages including description and plots A should not exceed 10 pages.

The writeup will receive a grade which reflects our assessment of your level of understanding of the various issues involved and the effort that you put into it.

Playing Sound in Matlab

Throughout these projects it will be necessary to play audio files which you will have manipulated in Matlab. The functions sound and soundsc work well for doing this on most hardware platforms.

For this project, the file projIB.mat will need to be first loaded into Matlab. To access the project zip file, click on the link below where you found this file on the web site.  After downloading the zip file, copy its contents into your Matlab working directory and type load projIB.

Matlab’s filter commands

In this part of the project, it may be helpful to use Matlab’s order estimation functions (e.g. buttord, cheb1ord, …). A caveat in the use of these functions is that Matlab’s definition of ‘ripple’ differs between the IIR and FIR filter design functions.  Do the Matlab Ripple project before proceeding with this project.

Filtering Noisy Speech

This project is concerned with designing a low-pass filter for the removal of high-pass noise from a speech signal. The noisy signal is stored in the variable noisy and was sampled at 44100 (fs) Hz. It consists of a summation of a speech signal which was band-limited at 4 kHz using a low-pass filter with a very narrow transition band, and a noise signal which was filtered at 4 kHz using a high-pass filter with a very narrow transition band. To filter the noise we must design a discrete-time filter with the following parameters:

1. Passband edge: 2500 Hz.
2. Stopband edge: 4000 Hz.
3. Maximum gain in the passband: 40 dB.
4. Minimum gain in the passband: 37 dB.
5. Maximum gain in the stopband: 55 dB.

Because it may take a few iterations to get each filter right, we suggest you write a .m Matlab script for each filter.

Note that in this part of the project, you may end up with very high order filters since the specifications are rather severe. You can try to implement the IIR filters according to the specs using Matlab’s built in tools (such as the fdatool or the command line tools butter, cheby1, cheby2, and ellip—type “help signal” for more information). If you do so, you will notice that the resulting systems might not be stable. There are several reasons for this, the most important of which being coefficient quantization—even with floating point precision. This issue arises because the specifications are very tight and some of the filter types have all of their poles concentrated near z = 1.

Fortunately, there is a workaround. What you need to do is to group the poles and zeros for the desired filter in pairs (conjugate or not) to create smaller stable second order filters of the form N(z)/D(z) where N(z) and D(z) are at most second order polynomials in z. A cascade of such filters will produce the desired system, and Matlab will be able to analyze it. The drawback is that you will have to implement your own methods to generate the plots and simulate the system.

Design a DT filter of each of the following types based on the specifications given above:

Butterworth, Chebyshev Type I, Chebyshev Type II, Elliptic, Parks-McClellan, Kaiser.

For each of the designs,

• Determine the order of the filter.
• Determine the number of multiplication operations per input sample required to implement the filter. Be sure to consider the structure you assume.
• Plot the magnitude response (in dB) from w = 0 to w = p  using freqz. Plot a detail of the magnitude response, focusing on the passband ripple (linear scale). Plot the group delay (in samples) using grpdelay. (Use subplot to fit the three plots on the same page for each filter.)
• Plot the pole-zero diagram.
• Plot the impulse response using filter and stem for 100 samples. (Use subplot to fit the pole-zero diagram and the impulse response on the same page.)

Filter noisy using your de-noising filter. Listen to the filtered and original files. How do they compare?

ECE 310 Project 4 Spectrogram Analysis and Applications

The frequency characteristics of nonstationary signals vary as a function of time. This requires Fourier analysis of these signals to be localized and time-dependent, resulting in analysis methods known as time-frequency distributions. One elementary time-frequency distribution based on the FFT is the short-time Fourier Transform (STFT). The STFT applied to speech is often called a speech spectrogram. Spectrograms are well suited to analyze speech signals with their time-varying narrowband features.

In this project, you will:

• Calculate spectrograms of frequency modulated signals.
• Construct narrowband and wideband spectrograms of speech signals.
• Estimate a speech signal from a modified STFT.
• Modify the time scale of a speech signal.

The project zip file contains the audio files and MATLAB functions which you will need for this project. To access the zip file, click on the link below where you found this file on the web site.

Spectrogram Definition

The spectrogram of a signal x[n] is defined as:

where w[n] is the window used. The window will determine what portion of the signal is used for analysis and controls the frequency resolution of the spectrogram. The parameter n denotes the reference position of the window on the signal. Let the window be of length N_win and nonzero only in the interval 0 ≤ k ≤ (N_win – 1). The above equation reduces to:

Computing Spectrograms in MATLAB

It is easy to construct spectrograms in Matlab since the STFT involves only windowing and the FFT. This project will require extensive use of the Matlab function spectrogram. It is suggested you take some time to fully understand the details of the spectrogram function. This can be done by typing help spectrogram at the Matlab prompt. For a more graphical presentation of the help (including figures) try doc spectrogram. A few important points about the spectrogram function for this project:

• Note that if the spectrogram is defined as it is in the previous section, the first column of the spectrogram generated by matlab is really , not X₀(w).
• The spectrogram function with no output arguments will plot the results to the current figure.
• For real signals x[n], the spectrogram function only calculates or plots a portion of the frequency components (from 0 to p) since the remaining are “redundant”. How are they “redundant”?
• It is common in speech processing to use a “short” time window and to plot frequency on the y axis. When you have loaded the speech signals mentioned below, the command spectrogram(s1,hamming(256),’yaxis’) will produce a good example spectrogram plot. Please use the ‘yaxis’ parameter for all your figures.

The MATLAB function diary may be helpful to keep track of commands when you are experimenting. Create separate M-file(s) to be able to rerun and test parts of your code. Do not turn in a log of all the commands you tried, just the source code that was relevant to solve the different parts of the project. In addition to this source code, please turn in the required plots, keeping in mind that you are expected to generate reasonable graphics.

Frequency-Modulated Signals

1. Spectrograms display frequency variations as a function of time. In this section, we will attempt to develop some intuition about the STFT by analyzing signals with known frequency variations, specifically linear frequency modulated (FM) chirp signals. The mathematical form of a continuous-time linear FM chirp is:

Generate a discrete-time linear FM chirp signal x[n] using a sampling frequency of 5 MHz and letting the chirp rate be m = 4.0 X 10⁹. Assume that the continuous-time chirp lasts 200 ms. Generate a spectrogram for this signal X_n(w) using 256-point FFT’s, a 256-point triangular window and an overlap of 255 samples between sections. (The MATLAB function triang can be used to construct the triangular window.)

2. Notice the distinct ridge in the time-frequency plot of part (1). Let us attempt to relate this to the functional form of the linear FM chirp signal. If the frequency modulation of x[n] is rewritten as , one could interpret f₁(t) as the “instantaneous frequency” of x(t). Another way to construct an “instantaneous frequency” of x(t) is the time derivative of the phase: where x(t) = cos [f(t)]. Calculate both of these definitions of “instantaneous frequency” for the signal in part (1) and determine which definition corresponds to the slope of the ridge in the spectrogram.
3. Sketch the spectrogram if the chirp rate was changed to m = 1.0 X 10¹⁰ . Verify your intuition by defining another chirp signal x₂[n] with this chirp rate and computing its spectrogram. Compare this spectrogram to the one obtained in part (1) and discuss your conclusions.

Narrowband and Wideband Spectrograms

The window used to compute the STFT determines the frequency resolution of the analysis. In this section, you will examine the effect of the window on the spectrogram of a voiced speech signal. You will need two speech files in the project zip file: s1.mat and s5.mat. Load in these files by typing load sl and load s5 at the Matlab prompt. Assume these speech signals are sampled at 8 kHz. If you wish to listen to these files, the Matlab function soundsc can be used to autoscale and play signals.

The exercises below depend on the frequency structure of human speech. Voiced human speech is produced by buzzing from the vocal cords (glottal excitation) into the resonant frequencies (formants) of the vocal tract. At high frequency resolution, we can see the individual harmonics of the glottal excitation signal at multiples of its pitch frequency.  At low frequency resolution, the individual harmonics are blurred, but the emphasis applied by the vocal tract resonances is more obvious in the wide peaks. As the lips, teeth and tongue change the shape of the vocal tract, the resonant frequencies (formants) change the sound to correspond to various recognizable vowel sounds. Thus, the specific formant frequencies determine the vowel sound perceived by the listener.

4. One application of a spectrogram is to estimate the fundamental frequency as a function of time. This estimation can be performed using a narrowband spectrogram which has fine frequency resolution and poor temporal resolution. Determine the proper spectrogram parameters (FFT length, window length and overlap between sections) to construct narrow band spectrograms of mat and s5.mat. Use a triangular window and FFT lengths that are powers of two. Experiment with different parameter values to obtain spectrograms that enable you to estimate and sketch the fundamental frequency of both speech signals as a function of time.
5. Another application of a spectrogram is to estimate the formant frequencies as a function of time. This estimate can be performed using a wideband spectrogram which has fine temporal resolution and poor frequency resolution. Determine the proper spectrogram parameters (FFT length, window length and overlap between sections) to construct wideband spectrograms of mat and s5.mat. Use a triangular window and FFT lengths that are powers of two. Experiment with different parameter values to obtain spectrograms that enable you to estimate and sketch the formant frequencies of both speech signals as a function of time.

Modified Short-Time Fourier Transforms

6. The STFT of a signal has a lot of redundant information. Note that an arbitrary function Xn(w) may not be a valid STFT. This can be seen by noting that overlapping portions of the input are used to construct X_n(w) for different values of n. Taking the inverse DFT of each X_i(w) may not result in the same values of x[n] for different values of i, therefore an arbitrary function X_n(w) may not be a valid STFT.

There are different methods to estimate a signal from a modified STFT. In this section, you will estimate the signal by determining the signal y[n] which produces the smallest error for the following criterion:

Write a Matlab function to implement this estimation from a modified STFT. To simplify the programming, assume that the STFT that you are processing was created with 1024-point FFT’s and a 256-point rectangular window with segments overlapping by 128 points.

One way to approach this problem is to look at how a STFT is constructed with the parameters above. The first column of the STFT is constructed by taking the first segment of the signal in time and computing its FFT. In this case, that would be the 1024 point FFT of the first 256 points zero-padded to length 1024. The second column of the STFT is constructed by taking the next segment of the signal in time and computing its FFT. Since the segments overlap by 128 points, this would mean the 256 points from 128 to 383 are zero-padded to length 1024, and the 1024 point FFT is computed.

Now, if we want to estimate the signal from a modified STFT, we reverse the concepts above. Using the first column of the modified STFT, we can take an IFFT and the first 256 points should correspond to the time signal from 0 to 255. The IFFT of the second column of the modified STFT can be computed and the first 256 points will correspond to the time signal from 128 to 383. And we can repeat this for the remaining columns. Note that there are now two estimates of the time signal from 128 to 255 since the IFFT’s of the first two columns overlap in time. Because of the error criterion we are using (and the rectangular windows), we can just take the average of the estimates at each point. If we continue doing this, note that except for 128 points at the beginning and the end of y[n], there are two estimates at every point in time. This illustrates the redundant information in the STFT.

The only input variable to this function, besides the modified STFT, will be the number of samples of the spectrogram (note: this could also be determined from the modified spectrogram itself). Also assume that the output signal y[n] will be real and that modifications to the STFT are made in a manner such that IDFT of each slice of X_n(k) will also be real. Even if conjugate symmetry is satisfied, taking the inverse DFT with the Matlab  function ifft may result in a nonzero imaginary component due to precision error. In this case, take the real part of the ifft result to eliminate the precision error. Also, note that the MATLAB function spectrogram will only compute X_n[k] for 0 ≤ k ≤ 512 since the original sequence was real, so you will have to reconstruct X_n[k] for 513 ≤ k ≤ 1024 to be able to compute a 1024-point inverse DFT.

Test this function by using the unmodified spectrogram of the speech signal in vowels.mat as an input to your program and comparing the output to the original signal. Plot the difference between the input and output.

Changing the Rate of Speech

7. There are many scenarios where changing the rate of speech may be useful. One example is slowing down speech for beginners trying to learn a new language. Simple schemes that temporally interpolate or decimate the speech signal will change the speech rate but also alter characteristics of the speech such as the fundamental frequency and formant frequencies. The altered speech characteristics can make the speech signal unintelligible. Ideally, one would like to be able to modify the rate of speech while preserving the speech characteristics. The STFT can be used to accomplish this. In this section, you will modify a STFT and then use the function that you wrote in the previous section to obtain a speech signal with a different time scale than the original signal, but with the frequency characteristics of the speech preserved.

First, load in the speech signal from vowels.mat and compute the spectrogram using a 256-point rectangular window with an overlap of 128 points per section, 1024-point FFT’s, and assuming a sampling rate of 8 kHz. Then compress the time scale by a factor of 2 by throwing out every other slice in time, and use the function from the previous section to obtain a faster speech signal with the same frequency content. Plot the speech signal before and after time scale modification. Listen to the signals to verify that your processing is working correctly.

ECE 310 Project 5 Spectral Estimation

In this project you will explore and compare various parametric and non-parametric techniques for power density spectrum (PDS) estimation. The goal is to approximate the PDS with non-parametric  direct methods, such as the periodogram, periodogram averaging, and the non-parametric  indirect Blackman-Tukey method. You will also investigate parametric all-pole modeling for PDS estimation. In the process, you will have a chance to sort through many issues involving finite length signals and the DFT.

Suggested Reading:  Sections 10.6, 10.7, and 11.0 to 11.5.

You are given a signal  that is the result of passing white Gaussian noise with variance one ( ) through a system .

 

Figure 1: Generating

The goal is to estimate the PDS of , . Recall that:

Thus estimating  is equivalent to estimating . Your goal is to construct estimates of  from the samples of the colored noise .

The data for this project is contained in the file pj2data.mat.

Copy this file into the local directory where you are working. This file contains two vectors. The vector y is a 512-point vector representing the discrete time signal . The vector Hejw2 is a 512-point DFT representing samples of the magnitude-squared response . It is the desired response that you are trying to estimate from . Hejw2 is given to you as a baseline for error calculation.

You can load the data into your MATLAB environment using the load command. Just cd to the directory where you copied pj2data.mat and at the prompt type:

load pj2data

Your MATLAB environment will now have two vectors y and Hejw2 defined in it. Warning! If you had other variables named y and Hejw2 they will be overwritten by the load command.

Figure 2: Plots of the data in pj2data.mat

As a first exercise you may want to plot the data y and Hejw2. They are plotted above for your convenience.  Note that the k-indices on Hejw2 represent sampled values in the interval  (i.e. ).

In this project you are using the data  to construct estimates  of , where .

Error Criterion: The error criterion is defined as the discrete-squared-error in the frequency domain. An expression for the estimation error is given in (1):

 

Autocorrelation Estimate: The autocorrelation function of a wide-sense stationary signal $y[n]$ is defined as:

 

In this project, the autocorrelation estimate provided by the observed data of length  is:

 

Project Write-Up

This project contains several problems, each of which address some aspect of spectral analysis. The write-up should include a summary of how you went about the project, any relevant plots and, at minimum, address all the questions asked in each of the problems. The total, including your description and plots, should not exceed 14 pages (i.e., 14 sides).

In addition, this project will be graded not only on content and correctness, but also on neatness and style of presentation.  Thus, although not a requirement, it is recommended that you typeset this report using Latex, Microsoft Word, or some other typesetting program.

You can save your figures off the figure window in MATLAB by clicking export and saving as .eps, .bmp, .tiff, etc.  If embedding images into MSWord, be warned that JPEG quality is low. Bitmaps or TIFFs might be better than JPEGs.

MATLAB Notes:

Functions you may find useful in this project include fft(), fftshift(), fliplr(), conv(), downsample(), xcorr(), toeplitz(), levinson(), freqz(), sum(), abs(), transpose(), ones(), zeros(), triang(). Look at the appropriate help files in MATLAB to see how to use them.

1. Autocorrelation and MATLAB

The function xcorr() in MATLAB calculates the finite autocorrelation between two sequences. Calculate the following sequence

= xcorr(y, y, ‘biased’)

This is similar to the convolution of two finite sequences  and .  For the first part of this exercise, show this fact using the following example:

1. Take the first 32 points of (we define this signal as ).
2. Plot the autocorrelation of using xcorr().
3. Then plot the autocorrelation using the conv()

The two plots you get should be the same, except for a scaling factor.

Problem A.1

What does xcorr calculate if you replace ‘biased’ with ‘unbiased’?

Problem A.2

• The deterministic autocorrelation of is

 

where is the length of . Explain why the Fourier transform of this autocorrelation function is a positive, real function.

• Take the 64-point DFT of using fft command in MATLAB. Plot the absolute value and the phase of the DFT of . What is the relationship between the 64 point DFT of  and the 64 point DFT of ?
• Use and generate a signal whose 64 point DFT, which is calculated by fft, is the same as the 64 point DFT of .

Problem A.3

The PDS is the Fourier transform of the autocorrelation function defined in . The Fourier transform of  is an estimate of the PDS, which is the 64-point periodogram of  using a rectangular window. The periodogram can be calculated directly from  without first convolving it with itself: Take a 64-point DFT of  and then find its magnitude squared. This will give you a 64-point DFT representing .

For the last part of this exercise, plot three figures:

• The absolute value of the 64-point DFT of , the deterministic autocorrelation sequence.
• The magnitude squared of the 64-point DFT of .
• The magnitude squared of the 64-point DFT of the first 64 points of .

What is the relationship between the three signals in A.3.a, A.3.b and A.3.c ?

B. Nonparametric PDS Estimation

Problem B.1

Estimate the PDS of  using only the first 32-points of the data. Do this by taking the 64-point periodogram. Plot the 64 DFT points of the desired frequency response, , and your estimate of the PDS at  on the same plot. Calculate the estimation error as given in .

Problem B.2

Estimate the PDS using all 512 points. Do this by taking the 1024-point periodogram of . Similar to Problem B.1, plot the 64 DFT points of the desired frequency response, , and your estimate of the PDS at those frequencies on the same plot. Calculate the estimation error as given in .

Problem B.3

Estimate the PDS using periodogram averaging. Write a script to find the 64-point periodogram of every 32 non-overlapping samples and average the results (i.e. there will be 512/32 = 16 periodograms to average). The  expression for the periodogram averaging is:

Similar to Problem B.1, plot the 64 DFT points of the desired frequency response and your estimate on the same plot. Comment on any differences between this estimate and the previous two. Calculate the estimation error as given in .

Problem B.4

In this problem you will use the indirect Blackman-Tukey method to estimate the PDS of .  The Blackman-Tukey method requires you to first estimate the autocorrelation of . Use the following three steps to generate your Blackman-Tukey estimate:

• Estimate the autocorrelation of using and .
• Truncate your autocorrelation estimate by using only for  (i.e. Do not use the estimate of autocorrelation for .)
• Take the 64-point DFT of the truncated autocorrelation function.

Similar to Problem B.1, plot 64 DFT points of the desired frequency response and your estimate on the same plot. Calculate the estimation error as given in .

Problem B.5

Format the errors you found in Problems B.1-B.4 in a table and discuss your results.

• Which method performed the best?
• Explain why the Blackman-Tukey method does so well without any averaging.
• Multiply your estimated autocorrelation obtained in part B.4.(2) with , a triangular window of size (nonzero for ) with unit center tap. Take a -point DFT of this new estimate. Calculate the estimation error as given in . Discuss why the performance is different than your result for problem (B.4).

C. Parametric PDS Estimation

Estimate the PDS of  using all-pole modeling. Here the stable impulse response  in Figure 1 is modeled with a class of stable filters with the following structure

 

.

The goal is to fit an all-pole model to the data and use that model to estimate .

Problem C.1

Formulate the Yule-Walker equations for 2nd through 7th order all-pole models, and solve for the coefficients to estimate . Use the values of  as estimated in Problem B.4. In a table, list the coefficients and the prediction error you found for . For which order is the prediction error minimum? (hint: check out the levinson command in MATLAB)

Draw the lattice implementation for .

Problem C.2

In addition to finding the denominator coefficients, you must also find the gain factor . In the Yule-Walker estimation method, the minimum mean-squared value for  in the ^th-order predictor is given by

Plot the estimation error as defined in as a function of . You may find the command freqz() helpful. Comment on any interesting features of this error curve. What does this imply about the system function ? What is your best estimate of the number of poles in ? Compare the prediction error to the estimation error.

Problem C.3

• Another method to find the gain factor is to use the Yule-Walker equations. In this case, the estimate of  is given by:

 

Show how this equation is obtained. Use  and find the estimation error in . Compare s with the estimates s in Problem C.2.

• Another estimator is defined based on the prediction error signal, . The estimator is

 

Explain how this estimator is obtained. Repeat C.3.a for this estimator.

Problem C.4

MATLAB contains functions called

• periodogram(),
• pwelch(),
• pyulear().

View the help files for these functions and explore their use in the context of Problems B.1, B.2, B.3 and C. Show how these commands can provide the same results which you obtained for these relevant problems . Do not use these functions to find answers in the previous problems in the project.