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

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).

Answer 1

[Yazan]

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
X = fftshift(fft(inpSigPad));
Y = fftshift(fft(outSigPad));
% 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'})

Answer 2

[Paul]

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[0] in terms of x[n] and h0. Then you can solve for h0. Repeat to define y[1] 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