%% Filter delay demo % Illustratation of delay induced by an LTI system. %% Start clear % close all %% Define system % The system is defined by the difference equation % % y(n) = 0.1 x(n) - 0.1 x(n-1) + 0.1 x(n-2) + 1.5 y(n-1) - 0.8 * y(n-2) b = 0.1*[1 -1 1]; a = [1 -1.5 0.8]; %% Plot frequency response % Use the 'freqz' function to calculate the frequency response. [H, om] = freqz(b,a); f = om / (2*pi); figure(1) clf subplot(2, 1, 1) plot(f, abs(H)) title('Frequency response magnitude') xlabel('Frequency') subplot(2, 1, 2) plot(f, angle(H)) title('Frequency response phase') xlabel('Frequency') %% Create sinusoid % Create a sinusoid with a specified frequency L = 100; n = 0:L; om1 = 0.15*pi; % Frequency of sinusoid radians/sample f1 = om1 / (2*pi); % Frequency in cycles/sample x = cos(om1*n); % Cosine signal figure(1) clf subplot(2,1,1) plot(n, x) title('Signal') ylim([-1.5 1.5]) %% Output signal % Comput output signal. % Plot input and output signals on the same plot. % Notice the output signal is delayed with respect to the input signal. % How many samples is the input signal delayed? y = filter(b, a, x); % Output signal figure(1) clf subplot(2, 1, 1) plot(n, x, n, y, '--') legend('Input signal', 'Output signal') ylim([-1.5 1.5]) %% Frequency response value at om1 % The value H(om1) of the frequency response at om1 is a complex number. I = sqrt(-1); z1 = exp(I*om1); d = length(b) - length(a); H_om1 = 1/z1^d * polyval(b, z1) / polyval(a, z1) % Value of frequency response at om1 %% Magnitude at om1 % |H(om1)| A = abs(H_om1) %% Phase at om1 % The angle of H(om1) will tell us the delay.. theta = angle(H_om1) %% Show magnitude and phase on frequency response plot figure(1) clf subplot(2, 1, 1) plot(f, abs(H), f1, A, 'o') title('Frequency response magnitude') xlabel('Frequency') subplot(2, 1, 2) plot(f, angle(H), f1, theta, 'o') title('Frequency response phase') xlabel('Frequency') %% Delay (in samples) at om1 % Compute the signal delay in samples M = -angle(H_om1)/om1 % Delay in samples fprintf('The delay is %.2f samples\n', M) %% Verify delay property % Plot the input and output signal together, where the input signal is % delayed by M samples and multiplied by A. % Observed that the output signal % y(n) is in phase with the delayed input signal A x(n-M). figure(1) clf subplot(2, 1, 1) plot(n + M, A*x, n, y) xlim([0 L]) legend('A x(n-M)', 'y(n)') title( sprintf('A = %.2f, M = %.2f', A, M) ) %% Plot phase delay function % Show delay M on plot of phase delay function. figure(1) clf plot(f, -angle( H )./om, f1, M, '.', 'markersize', 20) title('Phase delay') ylabel('Samples') xlabel('Frequency')