Course Work 1
Example 1.1
% Example 1.1
n = 0:40;
x = sin(n * pi /10);
stem(n, x);
title('x(n) = sin(n\pi/10)');
xlabel('n');
ylabel('x(n)');
print -djpeg99 images/eg1p1.jpg
eg1p1.jpg
Exercise 1.1
% Exercise 1.1
% i)
n = -20:20;
delta = inline ('n==0');
x = 0.5 * delta(n - 5);
stem(n, x);
title('x(n) = 0.5\delta(n - 5)');
xlabel('n');
ylabel('x(n)');
print -djpeg99 images/ex1p1.i.jpg
% ii)
n = 0:50;
x = cos(n*pi / sqrt(23));
stem(n, x);
title('x(n) = cos(n\pi / 23^{½})');
xlabel('n');
ylabel('x(n)');
print -djpeg99 images/ex1p1.ii.jpg
% iii)
n = 0:20;
x = (0.8).^n;
stem(n, x);
title('x(n) = 0.8^n');
xlabel('n');
ylabel('x(n)');
print -djpeg99 images/ex1p1.iii.jpg
% iv)
n:-10:10;
x = exp(j * 2 * pi * n / 10);
stem (n, x);
% stem generates Warning: Imaginary parts of complex X and/or Y arguments ignored.
title('x(n) = e^{j2\pi n/10}');
xlabel('n');
ylabel('x(n)');
print -djpeg99 images/ex1p1.iv.jpg
ex1p1.i.jpg
ex1p1.ii.jpg
ex1p1.iii.jpg
ex1p1.iv.jpg
gensin(A, f0, fs, phi, nstart, nend)
function x = gensin(A, f0, fs, phi, nstart, nend)
% This function generates a sampled sinusoid.
%
% gensin(A, f0, fs, phi, nstart, nend);
%
% Input arguments:
% A Amplitude.
% f0 Analogue frequecy in Hertz.
% fs Sampling frequency.
% phi Initial phase.
% nstart First sample to return.
% nend Last sample to return.
n = nstart:nend;
x = A * cos(2 * pi * (f0 / fs) * n + phi);
Exercise 1.2
% Exercise 1.2
x = gensin(1, 1, 2*sqrt(23), 0, 0, 50);
% x should be the same sine wave as in Exercise 1.1.ii.
n = 0:50;
stem (n, x);
% x turns out as expected.
title('gensin: A = 1, f0 = 1, fs = 2 x 23^{½}, \phi = 0, nstart = 0, nend = 50');
xlabel('n');
ylabel('gensin(n)');
print -djpeg99 images/ex1p2.jpg
ex1p2.jpg
Exercise 1.3
% Exercise 1.3
A = 1;
fs = 8 * 1000;
phi = 0;
nstart = 0;
nend = 100;
% i)
t = 0:1/fs:(nend - nstart)/fs; % Set up the continuous time axis.
f0 = [150; 300; 450; 600];
x = [ gensin(A, f0(1), fs, phi, nstart, nend);
gensin(A, f0(2), fs, phi, nstart, nend);
gensin(A, f0(3), fs, phi, nstart, nend);
gensin(A, f0(4), fs, phi, nstart, nend) ];
for i = 1:4
subplot(2, 2, i);
stem(t, x(i, :));
hold on;
plot(t, x(i, :), 'red');
title(['gensin: A = ', int2str(A), ', f0 = ', int2str(f0(i)),', fs = ', int2str(fs) ,',\newline\phi = ', int2str(phi), ',nstart = ', int2str(nstart), ', nend = ', int2str(nend)]);
xlabel('t /seconds');
ylabel('gensin(n)');
end
hold off;
print -djpeg99 images/ex1p3.i.jpg
% ii)
figure; % use clf to clear current figure, use figure to create a new figure. We are comparing 1st and 2nd set, so create 2 figures.
f0 = [7300; 7450; 7600; 7750];
x = [ gensin(A, f0(1), fs, phi, nstart, nend);
gensin(A, f0(2), fs, phi, nstart, nend);
gensin(A, f0(3), fs, phi, nstart, nend);
gensin(A, f0(4), fs, phi, nstart, nend) ];
for i = 1:4
subplot(2, 2, i);
stem(t, x(i, :));
hold on;
plot(t, x(i, :), 'red');
title(['gensin: A = ', int2str(A), ', f0 = ', int2str(f0(i)),', fs = ', int2str(fs) ,',\newline\phi = ', int2str(phi), ',nstart = ', int2str(nstart), ', nend = ', int2str(nend)]);
xlabel('t /seconds');
ylabel('gensin(n)');
end
hold off;
print -djpeg99 images/ex1p3.ii.jpg
% iii)
% Some of the graphs in pt ii look the same as the graphs in pt i, despite
% the fact that the original signals are different. This is because the
% signals in pt ii contain frequencies far above the Nyquist frequency
% which, in this case, is 4KHz. This is known as aliasing.
% iv)
% prediction: Aliasing will occur and the graphs will look the same as, or
% similar to, those in pts i and ii.
figure; % use clf to clear current figure, use figure to create a new figure. We are comparing 1st and 2nd set, so create 2 figures.
f0 = [25100; 25250; 25400; 25550];
x = [ gensin(A, f0(1), fs, phi, nstart, nend);
gensin(A, f0(2), fs, phi, nstart, nend);
gensin(A, f0(3), fs, phi, nstart, nend);
gensin(A, f0(4), fs, phi, nstart, nend) ];
for i = 1:4
subplot(2, 2, i);
stem(t, x(i, :));
hold on;
plot(t, x(i, :), 'red');
title(['gensin: A = ', int2str(A), ', f0 = ', int2str(f0(i)),', fs = ', int2str(fs) ,',\newline\phi = ', int2str(phi), ',nstart = ', int2str(nstart), ', nend = ', int2str(nend)]);
xlabel('t /seconds');
ylabel('gensin(n)');
end
hold off;
print -djpeg99 images/ex1p3.iv.jpg
ex1p3.i.jpg
ex1p3.ii.jpg
ex1p3.iv.jpg
sincinterpolate(x, T, t)
function y = sincinterpolate(x, T, t)
% This function does sinc interpolation on a sampled signal.
%
% sincinterpolate(x, T, t);
%
% Input arguments:
% x Sampled signal. - Must be a row vector.
% T Sampling interval.
% t time values to return. - Must be a row vector.
for tc = 1:size(t,2)
for n = 1:size(x,2)
z(tc, n) = x(n) .* sinc((t(tc) - (n*T)) / T);
%if n == 51
% fprintf(int2str(x(n)));
% fprintf(int2str(sinc((t(tc) - (n*T)) / T)));
%end
end
% y(tc) = sum(z);
% clear z;
end
y = sum(z');
%clear y;
%y = z;
sincinterpolate(x, T, t)
function y = sincinterpolate(x, T, t)
% This function does sinc interpolation on a sampled signal.
%
% sincinterpolate(x, T, t);
%
% Input arguments:
% x Sampled signal. - Must be a row vector.
% T Sampling interval.
% t time values to return. - Must be a row vector.
for tc = 1:size(t,2)
sincf = sinc((t -(t(tc)*T)) / T)
%for n = 1:size(t,2)
% sincf(n) = sinc((t(tc) - (n-1)*T) / T);
%end
y(tc) = x * sincf';
% y(tc) = sum(z);
% clear z;
end
%y = sum(z');
%clear y;
%y = z;
sincinterpolate(x, T, t) sincinterpolate(x, T, t) sincinterpolate(x, T, t)
function y = sincinterpolate(x, T, t)
% This function does sinc interpolation on a sampled signal.
%
% sincinterpolate(x, T, t);
%
% Input arguments:
% x Sampled signal. - Must be a row vector.
% T Sampling interval.
% t time values to return. - Must be a row vector.
for tc = 1:size(t,2)
sincf = sinc((t(tc) - ((0:(size(x,2)-1)) * T)) / T);
y(tc) = x * sincf';
end
Exercise 1.4
% Exercise 1.4
% i)
clf;
T = 1;
X = inline ('n==0'); % This is a delta function.
t = -5:0.1:5;
X = X(t);
x = sincinterpolate(X, T, t);
plot(t, x);
print -djpeg99 images/ex1p4.jpg
%hold on;
%stem(t, X);
% This appears to be broken!
% ii)
ex1p4.jpg
Course Work 2
Example 2.1
% Example 2.1
help fft;
% FFT Discrete Fourier transform.
% FFT(X) is the discrete Fourier transform (DFT) of vector X. For
% matrices, the FFT operation is applied to each column. For N-D
% arrays, the FFT operation operates on the first non-singleton
% dimension.
Exercise 2.1
% Exercise 2.1
% N DFT FFT Reduction
% 16 256 64 192
% 1,024 1,048,576 10,240 1,038,336
% 4.096 16,777,216 49,152 16,728,064
% All values are numbers of complex multiplications.
Example 2.2
% Example 2.2
% subplot(211) - Display a graph in the upper half of a figure.
% subplot(212) - Display a graph in the lower half of a figure.
help abs;
% ABS Absolute value.
% ABS(X) is the absolute value of the elements of X. When
% X is complex, ABS(X) is the complex modulus (magnitude) of
% the elements of X.
% abs - Display magnitude of complex signal.
help phase;
% PHASE Computes the phase of a complex vector
%
% PHI=phase(G)
%
% G is a complex-valued row vector and PHI is returned as its
% phase (in radians), with an effort made to keep it continuous
% over the pi-borders.
% phase - Display phase of a complex signal.
Exercise 2.2
% Exercise 2.2
% i)
n = 0:15;
x = exp((-j*n)/3);
plot(n, x);
% plot generates Warning: Imaginary parts of complex X and/or Y arguments ignored.
title('x = e^{-jn/3}');
xlabel('n');
ylabel('x');
print -djpeg99 images/ex2p2.i.jpg
% ii)
figure;
subplot(211);
plot(n, abs(x));
title('Magnitude spectrum of x = e^{-jn/3}');
xlabel('n');
ylabel('abs(x)');
subplot(212);
plot(n, phase(x));
title('Phase spectrum of x = e^{-jn/3}');
xlabel('n');
ylabel('phase(x)');
print -djpeg99 images/ex2p2.ii.jpg
% iii)
figure;
X = fft(x);
% Plotting a DFT
% In order to plot the DFT, X, that was produced with the fft command, we
% need to generate a frequency vector for the x axis. - Otherwise we will
% be plotting against sample number. The frequency vector that we generate
% will be in discrete-time frequency.
w = (2 * pi / length(n)) * [0 : (length(n) - 1)];
% This generate a frequency vector between 0 and 2pi.
% We can now plot the DFT thusly;
stem(w, abs(X));
title('DFT of x = e^{-jn/3}');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('DFT(x)');
print -djpeg99 images/ex2p2.iii.a.jpg
figure;
% We can also plot the DFT from -pi to +pi in the following manner...
% First, convert the frequency vector from (0 to 2pi) to (-pi to pi) by
% taking the (pi to 2pi) elements, converting them to (-pi to 0) and then
% placing them at the start of the vector;
mid=ceil(length(n)/2)+1; % Find the middle element of w.
w(mid:length(n))=w(mid:length(n)) - 2*pi; % Change (pi to 2pi) to (-pi to 0).
w = [w(mid:length(n)) w(1:mid-1)]; % Put w in numerical order again.
% We can now use the fftshift function to rearrange the samples in the DFT;
X = fftshift(X);
% Then replot the graph;
stem (w, abs(X));
title('DFT of x = e^{-jn/3}');
xlabel('Discrete-time frequency, -\pi to +\pi');
ylabel('DFT(x)');
print -djpeg99 images/ex2p2.iii.b.jpg
ex2p2.i.jpg
ex2p2.ii.jpg
ex2p2.iii.a.jpg
ex2p2.iii.b.jpg
Exercise 2.3
% Exercise 2.3
% i)
delta = inline('x==0');
n = 0:7;
for l = 0:7;
x = delta(n - l);
X = fft(x);
w = (2 * pi / length(n)) * [0 : (length(n) - 1)];
subplot(211);
plot(w, abs(X), 'o');
title(['Magnitude spectrum of DFT of x = \delta(n - ', num2str(l), ')']);
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('abs(DFT(x))');
subplot(212);
plot(w, phase(X), 'o');
title(['Phase spectrum of DFT of x = \delta(n - ', num2str(l), ')']);
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('phase(DFT(x))');
axis([0 6 -20 20]);
eval (['print -djpeg99 images/ex2p3.i.', char(97 + l), '.jpg']);
figure;
end
% ii)
% The DFT of an arbitrary length discrete-time impulse signal is;
% In the general case, the DFT of an arbitrary length discrete-time impulse
% signal is a unit step function.
% iii)
x(n+1) = 1
X = fft(x)
% A constant signal x(n) = 1 is a train of pulses of length n.
% The DFT of this signal is a train of zeros, except the first element.
% This element represents the DC value of the signal and is an impulse
% of height n.
% - i.e. a constant in the time domain turns into an impulse in the
% frequency domain.
% So, from pt ii, the DFT of an impulse is a step and, from pt iii, the DFT
% of a step is an impulse with height equal to the duration of the step.
% This explains the factor of 1/N when doing an inverse DFT.
ex2p3.i.a.jpg
ex2p3.i.b.jpg
ex2p3.i.c.jpg
ex2p3.i.d.jpg
ex2p3.i.e.jpg
ex2p3.i.f.jpg
ex2p3.i.g.jpg
ex2p3.i.h.jpg
Example 2.4
% Example 2.4
% Proakis & Manolakis pg 688 -> 690
% Table 9.1, pg 690.
% If the waveform is even then the DFT will be real.
Exercise 2.4
% Exercise 2.4
% i)
x = [1 1 1 0 0 0 1 1];
X = fft(x)
% X can be seen to be real.
% Plot the DFT in the standard way.
w = (2 * pi / length(x)) * [0 : (length(x) - 1)];
subplot(211);
stem(w, abs(X));
title('DFT of x = [1 1 1 0 0 0 1 1]');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('abs(DFT(x))');
subplot(212);
plot(w, phase(X), 'o');
title('Phase spectrum of DFT of x = [1 1 1 0 0 0 1 1]');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('phase(DFT(x))');
print -djpeg99 images/ex2p4.i.jpg
figure;
% ii)
x = [1 1 1 1 0 0 0 1];
X = fft(x)
% X can be seen to be complex.
% Plot the DFT in the standard way.
w = (2 * pi / length(x)) * [0 : (length(x) - 1)];
subplot(211);
stem(w, abs(X));
title('DFT of x = [1 1 1 1 0 0 0 1]');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('abs(DFT(x))');
subplot(212);
plot(w, phase(X), 'o');
title('Phase spectrum of DFT of x = [1 1 1 1 0 0 0 1]');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('phase(DFT(x))');
print -djpeg99 images/ex2p4.ii.jpg
figure;
% iii)
x = [1 1 1 0 0 0 0 1];
X = fft(x)
% X can be seen to be complex.
% Plot the DFT in the standard way.
w = (2 * pi / length(x)) * [0 : (length(x) - 1)];
subplot(211);
stem(w, abs(X));
title('DFT of x = [1 1 1 0 0 0 0 1]');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('abs(DFT(x))');
subplot(212);
plot(w, phase(X), 'o');
title('Phase spectrum of DFT of x = [1 1 1 0 0 0 0 1]');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('phase(DFT(x))');
print -djpeg99 images/ex2p4.iii.jpg
ex2p4.i.jpg
ex2p4.ii.jpg
ex2p4.iii.jpg
Exercise 2.5
% Exercise 2.5
% i)
n = 0:9;
x = sin((2 * pi * n) / 10);
X = fft(x);
% Plot the DFT from -pi to +pi.
w = (2 * pi / length(n)) * [0 : (length(n) - 1)];
mid=ceil(length(n)/2)+1;
w(mid:length(n))=w(mid:length(n)) - 2*pi;
w = [w(mid:length(n)) w(1:mid-1)];
X = fftshift(X);
subplot(211);
stem(w, abs(X));
title('DFT of x = sin(2\pin / 10), n = 0,1,...,9');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('abs(DFT(x))');
axis([-4 4 0 5.5]);
% ii)
n = 0:10;
x = sin((2 * pi * n) / 10);
X = fft(x);
% Plot the DFT from -pi to +pi.
w = (2 * pi / length(n)) * [0 : (length(n) - 1)];
mid=ceil(length(n)/2)+1;
w(mid:length(n))=w(mid:length(n)) - 2*pi;
w = [w(mid:length(n)) w(1:mid-1)];
X = fftshift(X);
subplot(212);
stem(w, abs(X));
title('DFT of x = sin(2\pin / 10), n = 0,1,...,10');
xlabel('Discrete-time frequency, -\pi to +\pi');
ylabel('abs(DFT(x))');
axis([-4 4 0 5.5]);
print -djpeg99 images/ex2p5.jpg
% The FT for both of these signals should result in an impulse pair at the
% positive and negative frequency of the sine wave. Instead, in the case
% where n = 0:10, spectral leakage has taken place in the DFT.
% The original waveform can be written as
% x = sin((2 * pi * n) / 10) * rect(n).
% The continuous FT of this function is found to be an impulse pair
% onvolved with the sinc function.
% When n = 9 the samples for the DFT coincide with the zeros of the sinc
% function. When n = 10 the DFT samples do not coincide with the zeros of
% the sinc function and spectral leakage results.
% The function has effectively been windowed and the spectral leakage is
% an artefact of the window.
ex2p5.jpg
Exercise 2.6
% Exercise 2.6
clf;
n = 0:31;
w = (2 * pi / length(n)) * [0 : (length(n) - 1)];
x = (0.9).^n;
plot(n,x);
Xhat = fft(x);
X = (1 ./ abs(1 - 0.9 * exp(-j * w)));
% Plot the DFT in the standard way.
plot(w, abs(Xhat), 'oblue');
hold on;
plot(w, abs(X), 'ored');
title('DFT of x = (0.9)^n u(n)');
xlabel('Discrete-time frequency, 0 to 2\pi');
ylabel('Magnitude');
legend('Xhat', 'X');
hold off;
print -djpeg99 images/ex2p6.jpg
% The differences between the two results are negligible. - They probably
% result from rounding errors and the slight variations of the
% implementations of the DFT algorithm in each case.
ex2p6.jpg
Course Work 3
dftr(b, a, x)
function y = dftr(b, a, x)
% This function computes the output of a linear system.
%
% dftr(b, a, x);
%
% Input arguments:
% x n input samples.
% a y coefficients.
% b x coefficients.
if (size(x,2) == length(x))
x = x'; % Ensure that x is a column vector.
end
if (size(a,1) == length(a))
a = a'; % Ensure that a is a row vector.
end
if (size(b,1) == length(b))
b = b'; % Ensure that b is a row vector.
end
x = [x; zeros(length(a),1)]; % Pad x with zeros.
if (a(1) == 1) % Ensure that a_0 = 1.
n = length(x);
x2 = zeros(length(b)-1, 1);
y2 = zeros(length(a)-1, 1);
for i = 1:n
x2 = [x(i); x2(1:length(b)-1)];
y(i, 1) = (b * x2) - (a(2:length(a)) * y2);
y2 = [y(i); y2(1:length(a)-2)];
end
else
fprintf('a(1) must = 1');
y = x;
end
Exercise 3.1
% Exercise 3.1
a = [1 3 5 7];
b = [8 10 12];
x = [1 2 3 4 5];
% First three elements of y should be 8, 2 and 10 respectively.
dftr(b, a, x)
% ans =
%
% 8
% 2
% 10
% -10
% 82
% -168
% 224
% -406
% 1274
% Output as expected.
Example 3.2
% Example 3.2
x = zeros(1024, 1);
x(1) = 1;
Exercise 3.2
% Exercise 3.2
x = zeros(1024, 1);
x(1) = 1;
% i)
a = [1 0];
b = [1, 1, 1, 1, 1];
stem(dftr(b, a, x));
title('Impulse repsonse of a = [1 0] and b = [1, 1, 1, 1, 1]');
print -djpeg99 images/ex3p2.i.jpg
figure;
% ii)
a = [1 0.9];
b = [1];
stem (dftr(b, a, x));
title('Impulse repsonse of a = [1 0.9] and b = [1]');
print -djpeg99 images/ex3p2.ii.jpg
figure;
% iii)
a = [1, (-1.8 * cos(pi/5)), 0.81];
b = [1, 0.5];
stem (dftr(b, a, x));
title('Impulse repsonse of a = [1, (-1.8 * cos(pi/5)), 0.81] and b = [1, 0.5]');
print -djpeg99 images/ex3p2.iii.jpg
ex3p2.i.jpg
ex3p2.ii.jpg
ex3p2.iii.jpg
Exercise 3.3
% Exercise 3.3
help roots
% A(z) = 1 + sum(a_m * z^(-m)) from m = 1 to m = N_a
a = [1, (-1.8 * cos(pi/5)), 0.81];
% A(z) = 1 - (1.8 * cos(pi/5)) * z^(-1) + 0.81 * z^(-2)
% roots requires positive powers of z. Therefore, multiply by z^2 on both
% sides;
% A(z) = z^2 - (1.8 * cos(pi/5)) * z + 0.81 = 0
p = roots ([1, (-1.8 * cos(pi/5)), 0.81])
% p =
%
% 0.7281 + 0.5290i
% 0.7281 - 0.5290i
magnitudep = abs(p)
% magnitudep =
%
% 0.9000
% 0.9000
% Magnitude of both roots is 0.9000
phasep = phase(p)
% phasep =
%
% 0.6283
% -0.6283
% Phase of roots is +/-0.6283 radians.
Exercise 3.4
% Exercise 3.4
% i)
% If the magnitude of at least one of the roots is greater than or equal
% to one then the filter will be unstable and the outut will only increase
% in size over time.
% ii)
a = [1, 2, 0.81];
p = abs(roots(a))
% At least one of the roots is >= 1.
b = [1, 0.5];
x = zeros(1024, 1);
x(1) = 1;
% x and b values from Exercise 3.2, pt iii.
stem (dftr(b, a, x));
title('Impulse repsonse of a = [1, 2, 0.81] and b = [1, 0.5]');
print -djpeg99 images/ex3p4.ii.jpg
% The output becomes enormous as expected and does not decrease.
% iii)
% In order to test the stability of a system the roots can be calculated
% and their magnitudes can be tested. If N_a is large (i.e. greater than
% about 5) it becomes very difficult to calculate the roots of the
% characteristic equation. It is often necessary to resort to
% computationally intensive numerical methods to discover the roots.
ex3p4.ii.jpg
Example 3.5
% Example 3.5
u = ones(1024, 1);
Exercise 3.5
% Exercise 3.5
u = ones(1024, 1);
a = [1, (-1.8 * cos(pi/5)), 0.81];
b = [1, 0.5];
stem (dftr(b, a, u));
title('Step repsonse of a = [1, (-1.8 * cos(pi/5)), 0.81] and b = [1, 0.5]');
xlabel('n');
ylabel('y(n)');
print -djpeg99 images/ex3p5.jpg
% i)
axis([0 100 0 7]);
title('Step repsonse of a = [1, (-1.8 * cos(pi/5)), 0.81] and b = [1, 0.5] for small n');
print -djpeg99 images/ex3p5.i.jpg
% For small n the filter produces an oscillating output that gradually
% decays away
% As n tends to infinity the filter output looks like the input only it has
% been amplified by a factor of about 4.2401.
% ii)
stem (dftr(b, a, u/2));
title('Repsonse of a = [1, (-1.8 * cos(pi/5)), 0.81] and b = [1, 0.5] for x = u/2');
xlabel('n');
ylabel('y(n)');
print -djpeg99 images/ex3p5.ii.a.jpg
% Halving the amplitude of the input halfs the amplitude of the output.
stem (dftr(b, a, 2*u));
title('Repsonse of a = [1, (-1.8 * cos(pi/5)), 0.81] and b = [1, 0.5] for x = 2u');
xlabel('n');
ylabel('y(n)');
print -djpeg99 images/ex3p5.ii.b.jpg
% Doubling the amplitude of the input doubles the amplitude of the output.
ex3p5.i.jpg
ex3p5.ii.a.jpg
ex3p5.ii.b.jpg
ex3p5.jpg
Example 3.6
% Example 3.6
% The output of a Finite Impulse Response (FIR) filter is determined
% solely by previous values of the input and the corresponding
% coefficients b_l. i.e. previous values of the output are not considered.
Exercise 3.6
% Exercise 3.6
x = [1 1 1 1 1];
a = [1 0];
b = [1 2 3];
% i)
% y(0) = 1x1 + 0x2 + 0x3 = 1
% y(1) = 1x1 + 1x2 + 0x3 = 2 + 1 = 3
% y(2) = 1x1 + 1x2 + 1x3 = 3 + 2 + 1 = 6
% y(3) = 1x1 + 1x2 + 1x3 = 3 + 2 + 1 = 6
% y(4) = 1x1 + 1x2 + 1x3 = 3 + 2 + 1 = 6
% y(5) = 0x1 + 1x2 + 1x3 = 2 + 3 = 5
% y(6) = 0x1 + 0x2 + 1x3 = 3
dftr(b, a, x)
% ans =
%
% 1
% 3
% 6
% 6
% 6
% 5
% 3
% ii)
% This operation is called convolution.
Example 3.7
% Example 3.7
help freqz
Exercise 3.7
% Exercise 3.7
N = 256;
% a)
a = [1 0];
b = [1, 1, 1, 1, 1];
[H, w] = freqz(b, a, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1 0] and b = [1, 1, 1, 1, 1] with N = 256');
print -djpeg99 images/ex3p7.a.jpg
figure;
% i) This is a high pass filter.
% b)
a = [1 0.9];
b = [1];
[H, w] = freqz(b, a, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1 0.9] and b = [1] with N = 256');
print -djpeg99 images/ex3p7.b.jpg
figure;
% i) This is a band pass filter.
% c)
a = [1, (-1.8 * cos(pi/5)), 0.81];
b = [1, 0.5];
[H, w] = freqz(b, a, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1, (-1.8 * cos(pi/5)), 0.81] and b = [1, 0.5] with N = 256');
print -djpeg99 images/ex3p7.c.jpg
figure;
% i) This is a high pass filter.
% ii)
a = [1 0];
b = [1, -1, 1, -1, 1];
[H, w] = freqz(b, a, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1 0] and b = [1, -1, 1, -1, 1] with N = 256');
print -djpeg99 images/ex3p7.ii.jpg
% The filter changes from a high pass to a low pass filter.
ex3p7.a.jpg
ex3p7.b.jpg
ex3p7.c.jpg
ex3p7.ii.jpg
Exercise 3.8
% Exercise 3.8
N = 256;
% i)
a = [1];
b = [1 2];
[H, w] = freqz(b, a, N, 'whole');
subplot(211);
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1] and b = [1 2] with N = 256');
ylabel('Magnitude');
subplot(212);
plot(w/(2*pi), angle(H));
title('Phase repsonse of a = [1] and b = [1 2] with N = 256');
ylabel('Phase');
print -djpeg99 images/ex3p8.i.jpg
% This is a notch filter.
figure;
% ii)
a = [1];
b = [2 1];
[H, w] = freqz(b, a, N, 'whole');
subplot(211);
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1] and b = [2 1] with N = 256');
ylabel('Magnitude');
subplot(212);
plot(w/(2*pi), angle(H));
title('Phase repsonse of a = [1] and b = [2 1] with N = 256');
ylabel('Phase');
print -djpeg99 images/ex3p8.ii.jpg
% This is a notch filter.
figure;
% iii)
a = [1 0 0.81];
b = [0.81 0 1];
[H, w] = freqz(b, a, N, 'whole');
subplot(211);
plot(w/(2*pi), abs(H));
title('Magnitude repsonse of a = [1 0 0.81] and b = [0.81 0 1] with N = 256');
ylabel('Magnitude');
subplot(212);
plot(w/(2*pi), angle(H));
title('Phase repsonse of a = [1 0 0.81] and b = [0.81 0 1] with N = 256');
ylabel('Phase');
print -djpeg99 images/ex3p8.iii.jpg
% This is a phase change filter.
% i and ii differ in what they do to the phase of the signal. The filter in
% pt ii has a discontinuous phase.
ex3p8.i.jpg
ex3p8.ii.jpg
ex3p8.iii.jpg
Example 3.9
% Example 2.9
help fir2;
% FIR2 FIR arbitrary shape filter design using the frequency sampling method.
% B = FIR2(N,F,A) designs an N'th order FIR digital filter with the
% frequency response specified by vectors F and A, and returns the
% filter coefficients in length N+1 vector B. Vectors F and A specify
% the frequency and magnitude breakpoints for the filter such that
% PLOT(F,A) would show a plot of the desired frequency response.
% The frequencies in F must be between 0.0 < F < 1.0, with 1.0
% corresponding to half the sample rate. They must be in increasing
% order and start with 0.0 and end with 1.0.
Exercise 3.9
% Exercise 3.9
% Need Attenuation > 55dB fr f < 1KHz. Therefore, set response = 0.
% Need Gain of about 0 for f > 1.4KHz. Therefore, set repsonse = 1.
% Set up the frequencies we want to specify;
F = [0 1/4 1.4/4 4/4]
% F =
%
% 0 0.2500 0.3500 1.0000
% Set up the response at the selected frequencies;
A = [0 0 1 1]
% A =
%
% 0 0 1 1
plot(F, A);
title('Desired frequency response of the filter.');
xlabel('Frequency / 0.5 * fsamp');
ylabel('Gain');
print -djpeg99 images/ex3p9.a.jpg
figure;
N = 256
b = fir2(N, F, A);
[H, w] = freqz(b, 1, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude response of the filter with N = 256.');
print -djpeg99 images/ex3p9.b.jpg
figure;
% OK, this works.
% Now try to reduce the number of taps;
N = 32;
b = fir2(N, F, A);
[H, w] = freqz(b, 1, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude response of the filter with N = 32.');
print -djpeg99 images/ex3p9.c.jpg
figure;
% This no longer works as well.
% Try a different window;
b = fir2(N, F, A, chebwin(N+1, 55));
[H, w] = freqz(b, 1, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude response of the filter with N = 32 and Chebyshev window.');
print -djpeg99 images/ex3p9.d.jpg
% Previous window was better.
% Try another;
b = fir2(N, F, A, bartlett(N+1));
[H, w] = freqz(b, 1, N, 'whole');
plot(w/(2*pi), abs(H));
title('Magnitude response of the filter with N = 32 and Bartlett window.');
print -djpeg99 images/ex3p9.d.jpg
% This is even worse.
% Use standard Hamming window.
ex3p9.a.jpg
ex3p9.b.jpg
ex3p9.c.jpg
ex3p9.d.jpg