Applied Signal Processing ETB006

Computer Assignments
Applied Signal Processing
ETB006
c Ingvar Claesson Mattias Dahl Sven Johansson
Translated to English by Jan Mark de Haan
October 2, 2003
Chapter 1
Elementary Signals and
Input-Output-Signal-Relations
Preparatory Assignment 1.1 Start your computer and create a sub-directory
named ”matlab” and start Matlab. Run all demonstrations you can find under
demo. Try the built-in Matlab Editor. Using the editor you can save your Matlab
files in your new sub-directory. You are free to use any other text editor. Generally
you do not write your commands directly in Matlab but you write a file which is
executed in Matlab. A Matlab command file has ‘extension’ .m and Matlab searches
in a path for these files. You executed a file by writing the file name without
‘extension’. The comment sign is % and comments that follow directly after the
first line will automatically be part of the help-system. Get acquainted with the
commands help and lookfor.
Assignment 1.1 Generate a time vector which contains the numbers 0…99. Create the unit impulse function δ(n) and the unit step function u(n) with zeros and
ones for the sample interval 0 ≤ n ≤ 99. Plot the signals in different intervals
with stem and plot respectively. Always use help to learn the command syntax
it’s function, for example for the command zoom.
Preparatory Assignment 1.2 Write down expressions for finite and infinite geometric series and give the requirements for convergence.
Assignment 1.2 Generate an exponentially decaying series x(n) = an · u(n) with
a = 0.5, 0.9, 1.0, 1.1, 2.0 for the previous time interval and create a plot. Elementwise operations are done with the point sign, for example should you use ˆ. in
this assignment. Also, calculate the sum of the samples without using loops and
compare the result with the expression for the geometric series. Use the command
input in this assignment.
Preparatory Assignment 1.3 What is a periodic signal? Give the conditions.
Assignment 1.3 Create the signal x(n) = cos(2πf0n) using a time vector of
length 8192. Plot the first 20-100 samples for the frequencies 0, 1/16, 1/8, 1/4, 1/2
and 1/√5. Did you expected the results? Are they periodic? Listen to the signals
using the command sound. Which frequency do you hear. Repeat the assignment
for a sinusoid with frequency 1/2 and explain the resulting plot.
1
Preparatory Assignment 1.4 Write the sum for two adjacent cosines with frequencies f0 and f1 as a product of terms with sum- and difference frequencies
Assignment 1.4 Continue by adding the signals x1(n) = cos(2πf1n) and x2(n) =
cos(2πf2n) in a vector of length 8192 with the adjacent frequencies f1 = 1/8 –
1/8192 and f2 = 1/8 + 1/8192. Listen and plot and rewrite your program using
the results from assignment 1.4. Then, listen only to x1(n). Change the frequency
f1 = 1/2 – 1/8192 and listen again. Can you explain the beating effect although
there is only one cosine???
Assignment 1.5 Plot x(n) = cos(2πf0n) for frequencies 1/8, 7/8 and 9/8. In
which way differ these functions?
Assignment 1.6 Create the complex signal x(n) = zn = anej2πf0n = elog(a)n+j2πf0n
where z = |z|ej arg z = |a|ej2πf0 for any a and f0 using a vector of length 8192. Plot
x(n) in a diagram. Also plot |x(n)| and the real- and imaginary part of x(n) in
an appropriate interval in the same figure, use the command hold. Listen to the
signals using sound. See which variables you have in Matlab using whos and
remove one or more variables with clear.
Assignment 1.7 Generate a chirp signal, i.e. a signal where the pitch (frequency)
increases gradually, and listen to it. Also, create sinusoids with beating-, vibratoand echo effects. Read ”Hallelujah” load handel, play it and add echo to it. Can
you change the pitch of the music?
Assignment 1.8 Generate a siren sound, i.e. a signal that alternates between
two tones. Also, try to generate a siren that glides back and forth between tones
(American siren).
Preparatory Assignment 1.5 Explain the terms impulse response and step response. Explain the system properties linear, causal and stable. Write down the
definitions and interpret these using your own words. Are the names of the definitions plausible?
Assignment 1.9 Generate a sequence x(n) = 1+cos(2πf0n) of 8192 samples with
f0 = 0.0625. Let a filter calculate the 1st difference, i.e. y(n) = x(n) – x(n – 1).
Use loops for..end or use subtraction with ‘two ‘vectors. Study the signal and
explain the result. What is the impulse response of the system calculating the 1st
difference? Is it causal?? Stable??
Assignment 1.10 Let a filter calculate the cumulative sum
y(n) =
n k
=0
x(k)
. Is this system stable?? Causal?? What is the impulse response h(n)??
2
Assignment 1.11 Finally, let a filter calculate finite cumulative sums,
y(n) =
n
k=n-4
x(k)
i.e. a gliding mean. Is this system stable?? Causal?? What is the impulse response
h(n)??
Assignment 1.12 Add 10δ(n – 10) – 10δ(n – 20) to x(n) = 1 + cos(2πf0n) for a
time vector with length 8192 where f0 = 0.025. Perform computation of a gliding
median,
n
y(n) = Median x(k)
k=n-4
i.e. a different gliding ‘mean ‘. Is this system stable?? Causal?? Linear?? What is
it’s impulse response??
Preparatory Assignment 1.6 Convolve the signal [1 2 3 4 5] with [1 3 2] with
a method of manual computation you have learned.
Assignment 1.13 Convolve signals [1 2 3 4 5] and [1 3 2] with the Matlab functions conv and xcorr. Also, calculate the correlation between the signals, both
with conv, and withxcorr.
Assignment 1.14 Redo assignments 1.9-1.11 by using conv. compare the and
convince yourself of good results by studying the samples and the vector dimensions
with whos or by plotting in one diagram. Use hold on and hold off, and plot
with different colors.
Assignment 1.15 Redo assignment 1.9-1.11 by using the Matlab function filter.
Compare the results and convince yourself that they are the same.
Assignment 1.16 Load ”Hallelujah” and filter it using the filter in assignment
1.9, and with the filter in assignment 1.11 using impulse responses of different
length. Can you hear the effect of the filter?
Preparatory Assignment 1.7 A system has step response s(n) = 0.5n cos(πn/6)·
u(n). Calculate the impulse response analytically.
Assignment 1.17 A system has step response s(n) = 0.5n cos(πn/6) · u(n). Determine the impulse response using Matlab and plot both the step response and
the impulse response in one and the same diagram. Compare with your analytic
expression.
Preparatory Assignment 1.8 A circuit has impulse response h(n) = 0.5n cos(πn/6)·
u(n). Calculate the step response analytically.
Assignment 1.18 A system has impulse response h(n) = 0.5n cos(πn/6) · u(n).
Determine the step response with Matlab and plot both the step response and
the impulse response in one and the same diagram. Compare with your analytic
expression.
3
Assignment 1.19 A helicopter has a rotor which rotates with 780 rotations per
minute. The rotor with accompanying propeller generates harmonics and together
they constitute helicopter sound (noise). An analysis shows that harmonics 6…15
dominate the sound. Generate helicopter sound and listen! What happens when
you remove one of the harmonics?
4
Chapter 2
Transforms and Frequency
Functions
Preparatory Assignment 2.1 A system has z-transform
1 – 0.5z-1
1 – z-1 + 0.5z-2
Calculate manually it’s poles, zeros and impulse response. What is the region of
convergence of the z-transform if we want it to be a causal system? Is the system
stable?
Assignment 2.1 A system has z-transform
1 – 0.5z-1
1 – z-1 + 0.5z-2
Determine and plot its poles and zeros using the command roots and plot its
impulse response using impz or filter. Compare with your previous results.
Preparatory Assignment 2.2 A system has z-transform
1
1 – z-1 – z-2 + z-3
Calculate its poles, zeros and impulse response. Use command residuez to verify
the partial fraction expansion. What is the region of convergence if the system is
causal? Is it stable?
Assignment 2.2 A signal has z-transform
1
1 – z-1 – z-2 + z-3
Determine and plot its poles and zeros and its impulse response. Compare with
your previous result.
5
Assignment 2.3 A signal has z-transform
1 – z-15
1 – z-1
Determine its poles and zeros and its impulse response. Can you simplify the
z-transform? Draw two realizations of the system, one based on the original expression and one based on the simplified expression. Which realization do you
prefer?
Assignment 2.4 Causality is an obvious property for systems working in realtime. However, there are many occasions where the input signal is available as
a whole, for example stored on a hard disk. Here the requirement of causality
can be ignored, for example tape recorder ECG. Here, zero phase filtering can be
employed where no delay will be introduced in the output signal compared to the
input signal. Load an ECG signal load ecg.mat and plot it. Filter this signal
partly once, partly twice with the filter
H(z) = 0.2
1 – 0.8z-1
first as it is (forward) and twice mirrored in origin (backward). Here you can use
fliplr or flipud. Also, check what the function filtfilt does and express the ‘equivalent ‘impulse response when you use this function. Recalculate the denominators
so that the output signals get the same amplitude as the input signal. Plot the
input signal and both output signals in one plot and try the command axis. Comment your results. If you prefer, you can apply the filter on a noisy ECG signal,
add some noise randn and try a different decay in your LP-filter H(z).
Assignment 2.5 What is the equivalent H˜(z) and its amplitude- and phase function when you use forward and backward filtering as in the assignment above??
Verify!
Preparatory Assignment 2.3 Determine the frequency function of H(z) above
analytically. Also, give expressions for the amplitude- and phase functions.
Assignment 2.6 Calculate the frequency function of H(z) with freqz, and verify
by programming your manually calculated expression. Also, calculate the amplitude and phase function and compare with your manually calculated expressions.
Assignment 2.7 Determine the impulse response of H(z) above in 4, 8, 16, 32
points and calculate the frequency function for this truncated filter. Compare with
the untruncated filter of the previous assignment.
Assignment 2.8 The Discrete Fourier Transform (DFT) is given by the expression
Hk =
N-1
n
=0
h(n) · ej2πkn/N
6
Write this relation as a matrix multiplication with a matrix of exponentials E (N ×
N) and vector h (1 × N) program your expression for arbitrary N, for the signal
in the assignment above.
Observe that the matrix E is orthogonal since EHE = N · I. Verify this and
utilize this fact to determine the matrix for the IDFT.
Assignment 2.9 When a signal h(n) has a short duration is it usual to add zeros
at the end of the signal (zero-padding) when the frequency function H(f) is calculated. Use the FFT function in Matlab fft to calculate and plot the frequency
function in 512 points for the truncated h(n) in the assignment above for the truncation lengths 4, 8, 16 and 32.
Assignment 2.10 You can measure the execution time for Matlab commands using tic . Repeat the previous assignment for the truncation lengths 7, 15, 31, . . . (214-
1). Compare the elapsed time for these lengths with the elapsed time for the truncation lengths 8, 16, 32, . . . (214). Study the number of floating point operations
with command flops.
Preparatory Assignment 2.4 Convolve thee signals [1 -1] and [1 0.5 0.25 0.125]
by manual calculation.
Assignment 2.11 Convolve signals [1 – 1] and [1 0.5 0.25 0.125] using DFT
and IDFT, partly with length N = 4, partly with N = 8. Study the result and
investigate which samples will have errors when you use N = 4.
Preparatory Assignment 2.5 Write down expressions for overlap-add techniques
for convolution and explain what the formulas mean.
Assignment 2.12 Load ”Hallelujah ” and filter the signal with truncated H(z) as
above with 64, 32… points. Perform the filter operation in the frequency domain,
partly with full length FFT and IFFT, partly with overlap add technique with
FFT length 64, 32 … without cyclic convolution. Choose the segments so that a
correct result is obtained. When do you hear the truncation effect? Verify that the
overlap-add method gives exact the same result as a full-length FFT.
7
Chapter 3
Measuring Frequency Content in
Signals
I this section we will study some non parametric methods for spectrum estimation
of a stochastic process. These methods are described in the literature.
All methods are based on the Periodogram which is defined for a sequence x(n)
with length N according to
P
xx( k
N
) = 1
N

N-1
n
=0
x(n)e-j2πnk/N

2
k = 0, 1, …, N – 1 (3.1)
This is an asymptotically unbiased estimate of the true spectral density. On the
other hand, it is not consistent, i.e. the estimate does not go towards zero when
the number of samples N → ∞. To avoid this, more advanced methods need to be
used. In this laboratory we will study three different methods, namely Bartlett,
Welch and Blackman/Tukey.
Bartlett: To reduce the variance of the periodogram proposed Bartlett to divide
the N samples into K segments with length M = N/K. For every segment, a
length M periodogram is calculated, whereafter the K periodograms are averaged,
i.e. the estimate becomes
P B
xx( k
M
) = 1
K
K-1
i
=0
1 M

M-1
n
=0
x(n + iM)e-j2πnk/M

2
k = 0, 1, …, M – 1 (3.2)
With that, the variance is reduced by a factor 1/K at the expense of reduced
frequency resolution, ∼ 1/M instead of ∼ 1/N.
Welch: To reduce the variance as much as possible, as many periodograms as
possible should be averaged. If the total number of samples N is given and a certain
frequency resolution is specified, the number of averages van be increased by letting
the segments overlap. To reduce the correlation between adjacent segments, each
segment is multiplied with a window w(n) before calculating the DFT. The result
are modified periodograms according to
P W
xx( k
M
) = 1
L
L-1
i
=0
1
MU

M-1
n
=0
x(n + iD)w(n)e-j2πnk/M

2
k = 0, 1, …, M – 1 (3.3)
8
where iD is the starting point for the i-th segment and N = LD. The factor 1/U
is a normalization factor with respect to the energy of the window
U = 1
M
M-1
n
=0
w2(n). (3.4)
With that, a lower variance is achieved, again at the expense of decreased frequency
resolution, which depends on M and the window. A positive property of the
window is that is reduces the Gibb’s phenomenon and gives increased dynamics in
the estimate.
Preparatory Assignment 3.1
a) An analog signal xa(t) = Aej2πF0t+φ is sampled using sampling frequency Fs =
1/T. Determine a N-point periodogram, see Eq. 3.1, of the signal under assumption that the signal has frequency F0 = l/NT, where l is an integer. Hint: use
geometric series.
b) An analog signal xa(t) = A sin(2πF0t + φ) is sampled using sampling frequency
Fs
= 1/T. Determine a N-point periodogram, see Eq. 3.1, of the signal under
assumption that the signal has frequency F0 = l/NT, where l is an integer. Hint:
use geometric series.
Assignment 3.1 Generate a sinusoid x(n) = A sin(ω0n) for N = 256 points, with
amplitude A = 1 and angular frequency ω0 = 2πk/N where k is integer. Plot the
signal, in order to verify it.
Assignment 3.2 Calculate and plot an N-point periodogram using fft and compare with preparatory assignment 3. Choose for example K = 63.
Assignment 3.3 Repeat the assignment with a sinusoid with frequency ω0 =
2π(k + 1 2)/N where k is integer. Choose for example k = 5 and k = 60.
Assignment 3.4 How large is the amplitude difference in the periodogram between a ‘at the spot ‘and a ‘in between ‘frequency, in linear and in dB scale. Use
max or look at plots to read the maximum values of the two sinusoid spectrum
estimates. How large can the error get if you do not window a sinusoid?
Assignment 3.5 Study the effect of multiplying signal x(n) with a window before
calculating the periodogram. Available windows are bartlett, hamming, hanning, blackman and flattop (see separate information about window functions
and additional m-file for flattop P301 specification on the course homepage). How
much can you reduce the error with these windows?
Assignment 3.6 Investigate the windows both in the time domain and the frequency domain. Plot the windows and change the length.
Assignment 3.7 Determine the signal amplitude A from the periodogram and
your window spectrum estimates.
9
Chapter 4
Measuring System, Time- and
Frequency Functions
The problem of system identification is illustrated in Fig. 4.1.
✚✙
✲ ✲✛✘ ✲

x(n)
h(n)
s(n)
w(n)
y(n)
Figure 4.1: Model for System Identification.
10
We assume that we can measure signals x(n) and y(n) and based on them
determine the unknown impulse response h(n). Signal w(n) is considered as a
disturbing signal. The resulting estimate hˆ(n) should be a FIR-filter of length M.
Deconvolution: This method is based on the convolution sum which, in absence
of the disturbing signal, can be written as
y(n) =
n k
=0
h(k)x(n – k) (4.1)
The last term in the sum h(n)x(0) is removed and after rearranging the other terms
the result is
hˆ(0) = y(0)/x(0) (4.2)
hˆ(n) =
y(n) –
n-1
k
=0
hˆ(k)x(n – k)
x(0) n ≥ 1 (4.3)
When the number of coefficients M for the filter has been chosen, these equations are used to determine the coefficients. Note that Eq. (4.3) is recursive, i. e.
hˆ(n) depends on hˆ(n – 1), …, hˆ(0).
Correlation: This method is based on the cross correlation between signals y(n)
and x(n). Under assumption that the disturbance w(n) is uncorrelated with the
input signal x(n), the following relation can be derived.

When we take the Fourier transform of this relation and rearrange, we obtain
H(ω) = Ryx(ω)
Rxx(ω) (4.5)
The estimate hˆ(n) for the unknown system is obtained by taking the inverse transform of H(ω)
hˆ(n) = F-1H(ω). (4.6)
The method is thus based on estimates of the cross spectral density Ryx(ω) and the
spectral density Rxx(ω) and is closely related to the spectrum estimation problem.
In practice, calculations are done using DFT (FFT) with a certain length M, which
also is the length of the estimate hˆ(n).
Least Squares Approximation: This method minimizes the quadratic distance
between the real output signal y(n) and the output signal ˆ s(n) which is obtained
when the input signal x(n) is filtered with hˆ(n). Under assumption that hˆ(n) is a
FIR-filter with length M we obtain the linear equations
M-1
k
=0
hˆ(k)rxx(l – k) = ryx(l) l = 0, 1, …, M – 1 (4.7)
11
which can be viewed as a special case of Eq. (4.4). These linear equations can be
compactly written as

In order to compare the different methods, a quality measure is needed. A
simple quality measure can be the (normalized) quadratic distance between the
impulse response h(n) and the estimated impulse response hˆ(n), i.e.
Qh =
M-1
n
=0
[h(n) – hˆ(n)]2
M-1
n
=0
[h(n)]2
(4.10)
A problem with this quality measure is that knowledge about h(n) is needed.
A quality measure which can be calculated without knowledge about h(n) exists
implicitly in the method of least squares, namely
Qy =
n
[y(n) – sˆ(n)]2
n
[y(n)]2 (4.11)
This quality measure neither assumes that the unknown system h(n) is a FIR-filter
but works generally.
Preparatory Assignment 4.1 Derive an expression for the signal-to-noise-ratio
(SNR) for the output signal y(n) in Fig. 4.1 assuming that x(n) and w(n) are
uncorrelated. SNR is defined as
SNR = Useful power in y(n)
Not useful Power in y(n)
Preparatory Assignment 4.2 In which cases are the quality measures Qh and
Qy in Eq. 4.10 and 4.11 equal?
4.1 Laboratory Assignments
Assignment 4.1 Matlab program DECON performs a system identification based
on the method of deconvolution. Execute the program for the disturbance free case
(SNR>1000). Vary the length of the impulse response estimate and study the properties of the deconvolution method by studying the quality measure Q. Plot the
true impulse response h and the estimate hhat.
12
Assignment 4.2 Executed program DECON with a disturbance (SNR10–50).
Vary the length M of the impulse response and investigate the properties of the
deconvolution method by studying the quality measure Q. Plot the rue impulse
response h and the estimate hhat.
Assignment 4.3 In Matlab programs CORR64,CORR128,CORR256 (see appendix), system identification is performed based on the correlation method with
different lengths for the FFT. Execute the program for the disturbance free case
(SNR>1000). Vary the length of the impulse response estimate and study the
properties of the deconvolution method by studying the quality measure Q. Plot
the true impulse response h and the estimate hhat.
Assignment 4.4 Execute CORR64,CORR128,CORR256 with a disturbance
(SNR10–50). Vary the length M of the impulse response and investigate the
properties of the deconvolution method by studying the quality measure Q. Plot
the rue impulse response h and the estimate hhat.
Assignment 4.5 In Matlab programs LS (see appendix), system identification
is performed based on least squares approximation. Execute the program for the
disturbance free case (SNR>1000). Vary the length of the impulse response estimate and study the properties of the deconvolution method by studying the quality
measure Q. Plot the true impulse response h and the estimate hhat.
Assignment 4.6 Execute LS with a disturbance (SNR10–50). Vary the length
M of the impulse response and investigate the properties of the deconvolution
method by studying the quality measure Q. Plot the rue impulse response h and
the estimate hhat.
Assignment 4.7 Make a concluding comparison of the three methods of identification.
13
Chapter 5
Adaptive Filters
An adaptive filter is a digital filter which is self-controlling, i.e. it follows automatically the changes occurring in the input signal, surroundings or desired signal.
Adaptive filters are the core for the subject of adaptive signal processing. In this
laboratory we will study a filter which is controlled by an algorithm called LMS.
This algorithm is a stochastic gradient algorithm which has very low computational
complexity. However, it has slow convergence. In many situations, this is not so
critical and methods exist to increase the convergence speed.
5.1 LMS Algorithm
The classical filtering problem is depicted in Fig. 5.1. An input signal is x[n]
is filtered by a filter f[n] which is to be designed so that the output signal y[n]
approximates the desired signal d[n] as much as possible, i.e.
ε[n] = d[n] – y[n] = d[n] – f[n] ∗ x[n] (5.1)
thus the goal is to minimize the error signal ε[n] by choosing f[n]. This can be
done by minimizing the Mean Square Error E(ε2[n]). When f[n] is a FIR filter,
this can be rewritten according to
ε[n] = d[n] – fTx[n] (5.2)
where f is a coefficient vector and x[n] is the input signal buffer at time index n.
The Mean Square Error can be expressed with the following equation
J = E(ε2[n]) = E(d2[n])-E(d[n]xT[n])f-fTE(x[n]d[n])+fTE(x[n]xT[n])f. (5.3)
Since this expression is real and scalar, fTE(x[n]d[n]) can be transposed without a
change. Let P = E(x[n]d[n]) be a cross correlation vector and R = E(x[n]xT[n])
be an autocorrelation matrix for the input signal of the filter. In order to find a
minimum for

dJ
df = -2P + Rf (5.5)
14
d[n]
x[n] y[n] ε
f[n]
+ –

✲ ✲ ✲
✧✦
★✥
Figure 5.1: Classical Filtering Problem
The gradient is set to zero and f is
f
opt = R-1P (5.6)
This solution is called Wiener-Hopf equation. This is the optimal filter, given
that the autocorrelation R and cross correlation P are known. This gives a static
solution to the problem. It is very time consuming to calculate the inverse of R .
We are interested in an adaptive solution, i.e. a solution that follows changes.
The usual solution, in order to avoid matrix inversion is to use a gradient
method, i.e.

The LMS algorithm does not need any knowledge about the autocorrelation R
and the cross correlation P, but uses actual samples directly. Basically, LMS
solves the linear equations 5.6 but it will never reach the optimal solution since the
expectation operator is removed. The LMS solution will approximate the optimal
solution.
15
Preparatory Assignment 5.1 We will estimate the system in the previous assignment using an adaptive filter f(n) = hˆ(n). Draw a block scheme for this
problem. Which is the desired signal d(n) in this case?
Assignment 5.1 Estimate P = E(x[n]d[n]) cross correlation vector and R =
E(x[n]xT[n]) using xcorr. Subsequently, solve the equation system
f
opt = R-1P
Assignment 5.2 Execute the Matlab program lms which solves the problem
adaptively. Vary the step size and plot the decaying error signal, and plot the
converging coefficients.
16

The post Applied Signal Processing ETB006 appeared first on My Assignment Online.

Expert Help With This Assignment — On Your Terms