Finding Impulse response of LTI system when input signal and output are given.
145 views (last 30 days)
Show older comments
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).
0 Comments
Answers (2)
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
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'})
0 Comments
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[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)
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)
But we also know that x[n] = 0 for n > 2. Therefore
y(n) = symsum(x(k)*h(n-k),k,0,2)
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
0 Comments
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!