Tutorial question
A digital filter is described by the expression:
y(n) = 1.5y(n-1) - 0.9y(n-2) + x(n) + 0.7x(n-1) + 0.6x(n-2)
(a) State whether the filter is recursive or non-recursive. Justify your answer.
: The filter is recursive because previous outputs appear in the equation, with the addition of
extra terms involving previous inputs (yn-1, yn-2 etc.).
(b) State the order of the filter.
: The filter is in the Second order.
(c) Derive the filter transfer function.
(d) From the coefficients of this transfer function, run the MATLAB program to compute the
first 100 samples of the impulse response using a sampling frequency of 10kHz, using the filter
function. Display the graph.
MATLAB code
close all;
clear all;
fs = ?; %sampling frequency
L = ?; %input length
y = zeros(L,1); %output signal initialization
a = [?]; %denominator coefficients
w = zeros(length(a),1); %deno internal states
b = [?]; %numerator coefficients
v = zeros(length(b),1);
x = eye(1,100);
output = filter(a,b,x);
figure;
plot(output); axis square;
axis([0 100 min(output) max(output)]);
xlabel('\itn');
ylabel('\ith(n)');
title('Impulse response');
MATLAB Updated code
close all;
clear all;
fs = 10000; %sampling frequency
L = 100; %input length
y = zeros(L,1); %output signal initialization
a = [1 -1.5 0.9]; %denominator coefficients
w = zeros(length(a),1); %deno internal states
b = [1 0.7 0.6]; %numerator coefficients
v = zeros(length(b),1);
x = eye(1,100);
output = filter(b, a, x);
figure;
plot(output);
axis square;
axis([0 100 min(output) max(output)]);
xlabel('\itn');
ylabel('\ith(n)');
title('Impulse response');
(e) From the impulse response, complete your program to compute the magnitude (in dB)
and phase (in °) spectra of the filter, using the fft, abs and angle functions. Then, display the
impulse response, magnitude response and phase spectrum in one graph.
MATLAB code
fft_result = fft(output); %Fourier transform
abs_result = abs(fft_result); %Absolute value
angleInRadians = angle(fft_result); %Angle in radians
dB = 20*log10(abs_result); %Conversion to decibel
angleInDegrees = radtodeg(angleInRadians); %Angle in degrees
(f) Modify your program to set the coefficient for a = [1 -1.9 0.9]; b = [1 0.7 0.6] and
display the corresponding graphs.
close all;
clear all;
fs = 10000; %sampling frequency
L = 100; %input length
y = zeros(L,1); %output signal initialization
a = [1 -1.9 0.9]; %denominator coefficients
w = zeros(length(a),1); %deno internal states
b = [1 0.7 0.6]; %numerator coefficients
v = zeros(length(b),1);
x = eye(1,100);
output = filter(b, a, x);
figure;
plot(output); axis square;
axis([0 100 min(output) max(output)]);
xlabel('\itn');
ylabel('\ith(n)');
title('Impulse response');
fft_result = fft(output);
abs_result = abs(fft_result);
angleInRadians = angle(fft_result);
dB = 20*log10(abs_result);
angleInDegrees = rad2deg(angleInRadians);
figure;
subplot(3,1,1);
plot(output);
title('Impulse Response');
subplot(3,1,2);
plot(dB);
title('Magnitude Response (dB)');
subplot(3,1,3);
plot(angleInDegrees);
title('Phase Spectrum (degrees)');