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'})
