George Mason University
Signals and Systems I
Spring 2018
Laboratory Project #4
Assigned: March 5, 2018
Due Date: Laboratory Section on Week of March 26, 2018
Lab Report
Your report for this lab will consist of all the analytical (i.e, pencil/paper) work,
MATLAB plots and code, and relevant explanations. Each student must do his or her
own work on this lab. However, you may ask other students or any of the teaching staff for
advice.
1
Prelab
(a) We have seen that an important class of linear time-invariant systems is one in which
the input and output are related by a linear constant coefficient differential equation,
such as
d2 y(t)
dy(t)
+2
+ 10y(t) = x(t)
2
dt
dt
or, more generally,
N
M
X
dk y(t) X dk x(t)
ak
=
bk
(1)
dtk
dtk
k=0
k=0
We have also seen that if the input to a linear time-invariant system is an exponential,
x(t) = est
where s = σ + jω is a complex number, then the output will be an exponential of
exactly the same form, but scaled by a complex number whose value depends on s
and the impulse response of the system, h(t). More specifically, the output is
y(t) = H(s)est
where
Z∞
H(s) =
h(t)e−st dt
−∞
is the system function. Thus, exponentials are eigenfunctions of linear time-invariant
systems. In the special case when s = jω, x(t) is a complex exponential of frequency
ω,
x(t) = ejωt
and
y(t) = H(jω)ejωt
In this case,
Z∞
H(jω) =
−∞
h(t)e−jωt dt
and H(jω) is called the frequency response of the system (filter). Clearly, H(jω) may
be found from the system function by setting s = jω,
H(jω) = H(s)
s=jω
For an LTI system described by a LCCDE of the form given in Eq.(1), we have seen
that H(s) is a ratio of polynomials
M
X
H(s) =
sM
sM −1
B(s)
bM
+ bM −1
+ · · · + b1 s + b0
=
= k=0
N
N
−1
N
A(s)
aN s + aN −1 s
+ · · · + a1 s + a0
X
bk sk
=
ak sk
P (s)
Q(s)
k=0
where the coefficients ak and bk are the coefficients of the differential equation. From
this, the frequency response is found by setting s = jω. The magnitude of H(jω) at
a given frequency ω tells us how much a complex exponential is attenuated or amplified in amplitude, and the phase (angle) of H(jω) indicates how much the complex
exponential is shifted in phase (delayed or advanced in time).
MATLAB has a number of useful m-files to find and plot the frequency response of
a system and to filter signals. One of these that you will be using in this lab is
freqs, which returns the frequency response H(jω) of a continuous-time system that
is specified by the coefficients ak and bk in the numerator and denominator polynomials
of H(jω). For example, if we have a linear time-invariant system defined by the
LCCDE
d3 y(t)
d2 y(t)
dy(t)
+
2
+ 10
+ y(t) = x(t)
3
2
dt
dt
dt
and if we enter the following commands in MATLAB ,
>> a = [1 2 10 1]; b = [1];
>> w=linspace(0,5,100);
>> H=freqs(b,a,w);
then H will contain the values of the frequency response in rad/s at the frequencies
specified by the vector w. In this case, H contains 100 samples of H(jω) within the
range 0 ≤ ω ≤ 5 rad/s. The magnitude of the frequency response may then be plotted
as follows,
>> plot(w,abs(H))
>> xlabel(’Frequency (rad/s)’)
>> ylabel(’|H(jw)|’);
and the plot will be as shown in Figure 1. There are other ways to use MATLAB and
these may be found by typing help freqs.
Other useful MATLAB functions that you will be using in this lab that you have used
before include tf(b,a), impulse(sys,tfinal,dt), and lsim(sys,u,t,x0). Again,
for documentation on any of these use the help function.
Figure 1: Plot of the magnitude of the frequency response of a filter with coefficients b=[1]
and a=[1 2 10 1]
.
(b) Knowing the systems’ frequency response H(ω) is particularly helpful when your
input is a continuous periodic signal. As you have studied, a periodic signal can be
represented as a Fourier Series
∞
X
x(t) =
Xn e nωo t
(2)
n=−∞
where
Xn =
1
To
Z
x(t)e−nωo t dx.
(3)
To
Notice that x(t) is now made up of a constant, Xn , times an ”everlasting” exponential,
enωo t . It is an easy process then to calculate the output to a system knowing its
frequency response:
∞
X
y(t) =
H( nωo ) Xn e nωo t
(4)
n=−∞
2
Periodic Signals
In this section we will focus on the output to an LTI system defined by its system differential
equation when the input is periodic. This will involve determining the Xn coefficients of
a Fourier Series expansion for use as an approximation to the input signal. The next step
will be to calculate the frequency response H( ω) for use in determining the output y(t).
2.1
Pulse Wave Signal
You are given a system whose linear constant coefficient differential equation is
d3 y(t)
d2 y(t)
dy(t)
+
59
+ 1951
+ 9179 y(t) = 9179 x(t)
3
2
dt
dt
dt
(5)
The input to this system is the periodic signal shown in Figure 2.
Pulse Wave
2
1
t
-1
1
2
3
4
-1
Figure 2: Plot of Periodic Input
(a) Calculate via hand analysis the Fourier Series coefficients Xn . Now write a MATLAB
m-script file to plot the Xn coefficients using the stem function
>> nmax=20;
>> n=[-nmax:nmax];
>> L_n=length(n);
>>
>>
>>
>>
>>
>>
%
%
%
%
omega_0=?;
%
%
X_n = zeros(1,L_n); %
for i=1:1:L_n
X_n(i) = provide
end
X_n((L_n+1)/2)=?;
%
%
You determine this number
Identifies the coefficients to calculate;
Coefficients to calculate.
Note: zero is at (L_n+1)/2
The radial frequency of the square wave
shown in Figure 2.
Initialize the coefficients
you calculation;
Provide a specific answer for X_o if
the calculation is faulty (e.g., 1/0).
>> figure(1)
>> stem(n,abs(X_n));
%Plot the magnitude of X_n
From the stem plot, determine nmax where Xnmax is greater than 1% of the largest
Xn . How fast do the Xn0 s fall off compared to n? (Hint: you may want to use the
”data cursor” under the ”tools” menu)
Instructor Verification (separate page)
(b) Approximate the input x(t) using the Fourier System coefficients X−nmax to Xnmax
that you calculated above.
>>
>>
>>
>>
>>
>>
t=0:0.005:2*To;
% Determine To from Figure 2.
x=X_n*exp(j*omega_0*n’*t);
figure(2)
plot(t,real(x)),xlabel(’t’), ylabel(’Partial Sum’)
% Plot real(x) to remove any small imaginary part due to
% computational error
axis([0 2*To -1.5 2.5]) %You provide the period To
text(0.05,-0.25,[’max. har. = ’ num2str(nmax)])
How does the Fourier Series approximation to the square wave compare with the
actual square wave? What discrepancies do you notice?
Instructor Verification (separate page)
(c) Determine H( ω) from Eq. (5). For the Fourier Series, ω will be replaced with n ωo
so that
H( ω) = H( n ωo ).
Now calculate the output of your system y(t) given the input x(t)
>>
>>
>>
>>
>>
>>
>>
for i=1:1:L_n
H_n(i)= ? ; Use your derived solution from Equation (5).
end
y=(H_n.*X_n)*exp(j*omega_0*n’*t);
figure(3)
plot(t,real(y))
% Be sure to label your plot!
axis([0 6 -1.5 2.5])
How does the output compare to the input? What type of filter is h(t)?
Instructor Verification (separate page)
Triangular Wave
2
1
-1
1
2
3
4
5
t
Figure 3: Plot of Second Periodic Input.
2.2
Triangular Wave Signal
In this section consider the triangular wave signal shown in Figure 3.
(a) Using hand analysis calculate the Fourier Series coefficients Xn .
(b) Plot the Xn coefficients using the stem function.
(c) What value of nmax did you find so that Dnmax is greater than 1% of the highest
Xn ? How does this compare with Part 2.1(a) above? How fast do the Xn0 s of the
triangle wave fall off compared to those of the pulse wave?
Instructor Verification (separate page)
3
Power Spectrum
A power spectrum shows where the power of a signal is distributed across the frequency
spectrum. A single sinusoid would be expected to have a majority of its power at the
frequency of the sinusoid while other signals may be more uniformly distributed. In this
section we’ll look at several wave forms to see how their power is distributed.
3.1
Sinusoid
Consider a vector x consisting of a unit-amplitude sinusoid with frequency fo = 330 Hz.
Given the sampling frequency fs = 33, 100 Hz, the spacing between samples would then
be dt = 1/fs = 3.02115x10−5 seconds. Set the number of samples of your waveform to
L = 216 = 65, 536.
(a) Write a MATLAB program to create the sinusoid x and plot the first 500 points. Next,
listen to the signal x using the soundsc command. The format for this function is
soundsc(x,fs). Failure to add the fs will mean you will be listening for a long time!
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
fs = 33100;
% Sampling frequency
fo = 330;
% Sinusoid frequency
dt = 1/fs;
% Sampling interval
L = 2^(16);
% Number of points in sinusoid
t = [0:dt:(L-1)*dt];
x = sin(fo*2*pi*t);
figure(5)
plot(t(1:500),x(1:500));
% Plot the first 500 points
xlabel(’Time (sec)’);
title(’Sine’);
ylim([-1.25 1.25])
% Expand y-axis
grid
soundsc(x,fs)
(b) The MATLAB code below computes the spectrum of the sinusoid x(t) and graphs the
spectrum from 0 to 3000 Hz.
>>
>>
>>
>>
>>
>>
>>
>>
>>
Fx = 20*log10(abs(fft(x,L)));
% Spectrum of signal x
ff = fs*linspace(0,1,L);
% Frequency values
figure(6)
plot(ff,Fx)
xlim([0 3000]) % Limits frequencies to 0 -> 3,000 Hz.
xlabel(’Frequency (Hz)’);
ylabel(’Spectrum Magnitude (dB)’)
title(’Sine’)
grid
What can you say about the distribution? Why isn’t all the power located at fo ?
3.2
Square Wave
(a) Next consider a square wave. MATLAB provides a function titled square(x) which
creates a unit square wave with frequency 2π Hz. Using the same parameters as given
in Section 3.1 above, use the code below to plot the first 500 points of the square
wave. Also, listen to the square wave signal using the MATLAB command soundsc.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
fs = 33100;
% Sampling frequency
fo = 330;
% Sinusoid frequency
dt = 1/fs;
% Sampling interval
L = 2^(16);
% Number of points in signal
t = [0:dt:(L-1)*dt];
x = square(fo*2*pi*t);
figure(7)
plot(t(1:500),x(1:500));
% Plot the first 500 points
xlabel(’Time (sec)’);
title(’Sine’);
>> ylim([-1.25 1.25])
>> grid
>> soundsc(x,fs)
% Expand y-axis
How does the square wave signal sound when compared to the sinusoidal signal in
Section 3.1(a)?
(b) Use the MATLAB code below compute and graph the spectrum of the square wave.
>>
>>
>>
>>
>>
>>
>>
>>
Fx = 20*log10(abs(fft(x,L)));
% Spectrum of signal x
ff = fs*linspace(0,1,L);
% Frequency values
figure(8), plot(ff,Fx)
xlim([0 3000]) % Limits frequencies to 0 -> 3,000 Hz.
xlabel(’Frequency (Hz)’);
ylabel(’Spectrum Magnitude (dB)’)
title(’Square Wave’)
grid
How does this spectrum compare to that found for the sinusoid in Section 3.1(b)?
3.3
Speech
In this section we’ll look at the spectrum of a noisy signal. Noise may be defined as an
unknown and undesirable signal, n(t), that is added to a desired signal x(t), i.e.,
y(t) = x(t) + n(t)
An important application of signal processing is to reduce the noise in y(t) by filtering out
as much of the noise as possible. Of course, in able to do so, something must be known
about the properties of the noise and/or signal. In this lab, you will investigate ho one
might remove additive noise from a speech.
(a) Load lab04_signal.mat into MATLAB as shown below. The following two variables
will be loaded, signal and fs .
load lab04_signal.mat
To listen to the signal, use soundsc(signal,fs) being sure to include fs . Describe
what you hear.
(b) Plot the magnitude spectrum of the signal using the same format as before.
>>
>>
>>
>>
>>
>>
L=length(signal);
Fx_filtered=20*log10(abs(fft(signal,L)));
ff=fs*linspace(0,1,L);
figure(9)
plot(ff,Fx_filtered);
% Be sure to label your plot!
xlim([0 fs/2]);
Compare this spectrum to both the sinusoid and square wave. What can you say
about the frequency distribution of this signal?
(c) There are a number of standard filters that have been created over the years. The
Butterworth filter is one that has also been made into a function in MATLAB called,
appropriately, butter. It consists of three parameters butter(n,wc,’ftype’) where
n is the order of the differential equation, wc is the cutoff frequency, and 0 f type0 is
the type of filter: ’low’,’high’, or ’stop’. We will be using the ’low’ pass filter. Use
the code below to plot H(jω) for the filter being used. Next apply the filter as shown
in the code and listen to the results.
>>
>>
>>
>>
>>
>>
wc = 0.5;
% Cutoff frequency is the range 0 < wc < 1.
[b,a] = butter(5,wc,’low’) ;
figure(10)
freqs(b,a,100);
% Plots to characteristics of H(jw)
signal_filtered = filter(b,a,signal);
soundsc(signal_filtered,fs);
How does this compare with the original signal? Try several cutoff frequencies (wc)
and note any differences.
(d) Finally, calculate and plot the spectrum of the filtered signal using the code below.
>>
>>
>>
>>
>>
>>
L=length(signal_filtered);
Fx_filtered=20*log10(abs(fft(signal_filtered,L)));
ff=fs*linspace(0,1,L);
figure(11)
plot(ff,Fx_filtered);
xlim([0 fs/2]);
How does this differ from the original signal’s spectrum?
Lab #4
ECE 220: Spring 2018
Instructor Verification Sheet
For each verification, be prepared to explain your answer and respond to other related questions that the lab TA’s or professors might ask. Turn this page in along with your lab report.
Name:
Date of Lab:
Part 2.1(a): What value of nmax did you find so that Xnmax is greater than 1% of the
highest Xn ?
Part 2.1(b): How does the Fourier Series approximation to the pulse wave compare with
the actual square wave? What discrepancies do you notice?
Part 2.1(c): What was the impact of passing the pulse wave through the system defined by
the differential equation in (5)?
Part 2.2: What value of nmax did you find so that Xnmax is greater than 1% of the highest
Xn ? How does this compare with Part 2.1(a) above? How fast do the Xn0 s of the triangle
wave fall off compared to those of the square wave?
Lab Report Format
Signals and Systems 1
ECE 220 Spring 2018
Each lab report that you submit for this class is a formal report that is well-prepared, carefully
written, and follows the format given below. The report must include descriptions and analysis
of all assigned portions of the lab. Lab reports are to be typed and submitted through the class
Blackboard website. Once you have completed your report, please submit it via Blackboard prior
to the start of the next lab. Thus, Lab #1’s report is due prior to the start of Lab #2.
General Guidelines
• Your report should be typewritten, neat and well-organized.
• The report must follow the format given below. All analytical work and calculations made
should be clearly explained.
• MATLAB code is to be placed in report appendices.
• It is expected that the report will be grammatically correct with no spelling errors. Points
will be deducted for poor grammar and/or spelling errors. You are encouraged to use a spell
checker.
• Explanations that describe your work must be included in the report. Plots that are submitted
within the report must include properly labeled axes and a title. Each plot should be given a
figure number and a caption. Points will be deducted if you include unnecessary graphs and
plots.
• Within the text of your report, when referring to a particular plot, refer to it by number.
Detailed Lab Report Format
Your report must contain the following sections. Subsections may be added for purposes of
clarification.
1. Title Page: This page contains an identification of the Laboratory by number and title. It
also must include the name and G-Number of the author and the date of submission.
2. Introduction: This section contains a description of the purpose and objectives of the
laboratory project. Do not simply copy text from the assignment. The Introduction should
summarize the topics covered in the laboratory and a brief summary of the key results
obtained.
3. Main Body: This section is the most detailed part in your report and it will generally contain
subsections. Note that it should not contain any MATLAB code. What it does contain is a
description of all results obtained during the lab as well as any figures that were generated
with MATLAB or other sources. If theoretical calculations are required, these should be detailed
and, where appropriate, compared with the experimental portion of the lab.
4. Conclusions: In this Section, summarize the conclusions that you made once the results
were obtained. The conclusions should be tied back to the objectives of the Lab Project.
This should be a concise section that focuses on the important results obtained and lists the
conclusions that resulted from these results.
5. References/Appendices: Include a description of any external references that you used
(e.g. web pages, technical papers, etc.) during the lab. The appendices should also contain
listings of your MATLAB code. You should document your code with comment lines, where
appropriate. This holds for both the primary code as well as any scripts and subroutines that
you have written. Remember, if a grader has trouble understanding your MATLAB code,
you will not receive full credit for that portion of your report. Adding carefully chosen
comments within the code will make it much more likely that your code is fully understood.
1
Purchase answer to see full
attachment