# Finding Impulse response of LTI system when input signal and output are given.

315 views (last 30 days)
HRITIK YADAV on 12 Aug 2021
Edited: Paul on 13 Aug 2021
Hi, I'm new to Matlab and have no idea how to solve the following problem, Thanks in advance for your kind explanation.
A signal sequence x(n) = {1,1,1} is applied to a LTI system with an unknown impulse response h(n), the observed output was y(n) = {1,4,8,10,8,4,1} write a matlab program to find h(n).

Yazan on 12 Aug 2021
Edited: Yazan on 12 Aug 2021
Below is a code that estimates the impulse response with 3 equivalent approaches:
clc, clear
inpSig = [1 1 1];
outSig = [1, 4, 8, 10, 8, 4, 1];
% pad the input and ouput signals with zeros to a common length (1024 samples)
comLength = 1024;
inpSigPad = [inpSig, zeros(1, comLength - length(inpSig))];
outSigPad = [outSig, zeros(1, comLength - length(outSig))];
% Approach 1
% estimate the impulse response directly using Matlab
% this function uses Welch's averaged periodogram to estimate the signal's psd
% I am using here rectangular windows with 0% overlap (equivalent to regular fft)
[H1, f] = tfestimate(inpSigPad, outSigPad, rectwin(comLength), 0, [], 1, 'centered');
% Approach 2
% You can estimate the cross and auto psd using Matlab
% assume a unitary sampling frequency
pxx = cpsd(inpSigPad, inpSigPad, rectwin(comLength), 0, [], 1, 'centered');
pxy = cpsd(inpSigPad, outSigPad, rectwin(comLength), 0, [], 1, 'centered');
% definition of transfer function
H2 = pxy./pxx;
figure,
plot(f, mag2db(abs(H1)));
hold on
plot(f, mag2db(abs(H2)), '--')
% Approach 3
% estimate the psd of input and output using fft
% reorder samples to have frequency response between (-fs/2, fs/2]
X = [X(2:end), X(1)];
Y = [Y(2:end), Y(1)];
H3 = Y(:)./X(:);
plot(f(1:30:end), mag2db(abs(H3(1:30:end))), 'o'); grid minor
xlabel('Frequency - Hz'), ylabel('Magnitude - dB')
legend({'Using tfestimate', 'using cpsd', 'using FFT'}) Paul on 13 Aug 2021
Edited: Paul on 13 Aug 2021
We know that y[n] = h[n] * x[n] where * denotes convolution. You have x[n] and y[n]. Let h[n] be denoted as [h0 h1 h2 etc.]. Use the defnition of the convolution sum to define y in terms of x[n] and h0. Then you can solve for h0. Repeat to define y in terms of x[n], h0, and h1. Solve for h1. Repeat until done. Recall that the convolution sum is defined as:
syms n k integer
syms y(n) h(n) x(n)
y(n) = symsum(x(k)*h(n-k),k,-inf,inf)
y(n) = We can make a standard assumption that x[n] = y[n] = 0 for n < 0 absent any other information on the indexing of x and y. Therefore
y(n) = symsum(x(k)*h(n-k),k,0,inf)
y(n) = But we also know that x[n] = 0 for n > 2. Therefore
y(n) = symsum(x(k)*h(n-k),k,0,2)
y(n) = Now you can write a program to solve this equation for the unknown h[n]. In doing so, you'll have to figure out what values to assign to h[-1] and h[-2], and also account for the fact that matlab arrays use 1-based indexing but the development above assumes that the first value in the arrays x and y correspond to n = 0.
Once you get a result, you can check that it's correct using
doc deconv