Download Report

George Mason University
Signals and Systems I
Spring 2015
Laboratory Project #4
Assigned: March 18, 2015
Due Date: Laboratory Section on Week of April 06, 2015
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
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
+ 2y(t) = x(t)
dt2
dt
or, more generally,
N
X
k=0
M
dk y(t) X dk x(t)
ak
=
bk
dtk
dtk
(1)
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
A(s)
aN sN + aN −1 sN −1 + · · · + a1 s + a0
X
bk sk
ak sk
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
+2
+ y(t) = x(t)
3
2
dt
dt
dt
and if we enter the following commands in MATLAB ,
>> a = [1 2 2 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(j\omega)|’);
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.
1
|H(jω)|
0.8
0.6
0.4
0.2
0
0
0.5
1
1.5
2
2.5
3
Frequency (rad/s)
3.5
4
4.5
5
Figure 1: Plot of the magnitude of the frequency response of a filter with coefficients b=[1]
and a=[1 2 2 1]
.
2
De-Noising a Trumpet Signal
Noise may be defined as an unknown and undesired 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 the signal. In this lab, you will investigate how
one might remove additive noise from a trumpet signal.
(a) Load lab04.mat into MATLAB . In this file will will find the following two signals:
1. trumpet
2. trumpet_n1:
The file trumpet contains the waveform of a trumpet playing the note middle C, which
is nominally 261.626 Hz. The file trumpet_n1 contains recordings of the trumpet
playing middle C with additive noise corrupting the signal. Use soundsc to listen to
both of these files. The noise should be clearly audible in trumpet_n1. Note that the
sampling rate for both signals is 44.1 kHz, so to listen to these sounds properly, you
need to specify to correct sampling frequency as follows,
>> fs=44100;
>> soundsc(trumpet,fs)
(b) We have seen that a trumpet has ten or more significant harmonics. Shown in Fig. (2)
is the spectrum of the trumpet playing middle C in the absence of any noise. This
plot was generated from trumpet as follows,
>> L=2^(18);
>> V1=abs(fft(trumpet,L));
>> trumpet_spec=V1(1:L/2);
Trumpet C4 (Middle C)
Spectral Magnitude
1000
800
600
400
200
0
500
1000
1500
Frequency (Hz)
2000
2500
3000
Figure 2: Magnitude of the spectrum for a trumpet as a function of frequency in Hz.
The array trumpet_spec contains N = 217 equally spaced samples of the magnitude
of the spectrum between f = 0 and f = 22.05 kHz. Fig. (2) shows a plot of the
spectrum of the trumpet for frequencies between zero and 3000 Hz. This plot was
generated as follows,
>>
>>
>>
>>
>>
>>
>>
ff=(fs/2)*linspace(0,1,L/2);
plot(ff,trumpet_spec)
axis tight;
xlim([0 3000]);
title(’Trumpet C4 (Middle C)’);
xlabel(’Frequency (Hz)’);
ylabel(’Spectrum Magnitude’);
Find the magnitude of the spectrum of the noise-free trumpet trumpet and the noisy
trumpet trumpet_n1, and make a properly labeled plot of each for frequencies in the
range 0 ≤ f ≤ 5000 Hz1 and make sure that your x-axis is labeled in Hz. Compare
the spectrum of the noisy and the noise-free trumpet signals and describe what you
observe.
(c) The MATLAB command
>> [b,a]=butter(N,wc,’s’);
creates an LTI system that is defined by an N th-order differential equation of the
form
N
X
dk y(t)
ak
= x(t)
dtk
k=0
Plot magnitude of the frequency response of the LTI system that is designed using
butter with N=9 and wc=2*pi*3000 as follows:
>> w=2*pi*ff;
>> wc=2*pi*3000;
1
Note that this is over a larger frequency range than the one shown in Fig. (2).
>>
>>
>>
>>
>>
>>
>>
[b,a]=butter(9,wc,’s’);
H=freqs(b,a,w);
plot(ff,abs(H));
grid on; axis tight;
xlim([0 10000]);
title(’Butterworth (LP) Filter’);
xlabel(’Frequency (Hz)’); ylabel(’|H|’);
What happens as you increase and decrease the value of wc used in the call to butter?
(d) As we see from the results of part (c), the filter that is designed using butter is a
lowpass filter because it passes, with very little attenuation, frequencies in x(t) that
are small compared to the cutoff frequency wc, and it attenuates those frequencies in
x(t) that are large compared to wc. With your filter that you designed in part (c)
with wc=2*pi*3000, use the MATLAB commands tf and impulse to find the impulse
response h(t), of your filter. Note that it is important to make sure that the sampling
period is specified when you make the call to impulse, such as
>> th=0:1/fs:0.2;
>> h=impulse(sys,th);
Plot of h(t) for 0 ≤ t ≤ 5 msec and make sure that your plot is properly labeled.
(e) Filter the noisy trumpet signal trumpet_n1 with your filter h designed in part (d)
using conv and listen to the output. Describe what you hear. Compare the signal
trumpet and your filtered noisy trumpet signal in the time domain and frequency
domain and comment on what you observe. Note: In order to make it clear what
the differences are, if any, between the original trumpet, the noisy trumpet, and the
filtered trumpet signals, plot these signals over a range of time no longer than 25 msec
and make sure that you look within an interval where the trumpet is playing.
(f) Experiment with the cutoff frequency wc in your filter and describe the differences in
the filtered trumpet signals. Is there any particular value for wc that seems to be the
best?
(g) Can you think of another type of filter that might do a better job in removing the
noise than the one you designed in part (c) using butter?