95 views (last 30 days)

I have measured time signals of input (pink noise) and output of a system. I need to calculate TF of a system and it's impulse response.

In order to calculate TF I can do in Matlab directly

TF = tfestimate(in, out)

or compute `fft` of both then do

TF = fft(out)./fft(in).

At this point everything is clear and expected for me, but then when I calculate Impulse response of the system by doing

ifft(fft(out)./fft(in))

I get impulse response which has a tail at the end. I was thinking maybe something wrong with my measurements and decided to repeat experiment in matlab synthesizing exponential chirp and passing it trough a filter with following impulse response calculation.

Then I get the "same" tale at the end:

Please see the code below

t = 0:1/48000:10;

fo = 1;

f1 = 24000;

x = chirp(t,fo,10,f1,'logarithmic');

x = x(:);

y = lowpass(x,1000,48000);

y = y(:);

% playing with FFT length here

% nfft=length(x)*2-1;

% nfft = 4*1024;

nfft=length(x);

tf = tfestimate(x, y, hann(nfft), nfft/2, nfft, 48000);

if mod(length(tf),2)==0 % iseven

tf_sym = conj(tf(end:-1:2));

else

tf_sym = conj(tf(end-1:-1:2));

end

ir_from_tfestimate = ifft([tf; tf_sym]);

in = fft(x, nfft);

out = fft(y, nfft);

tf1 = out./in;

ir_from_fft = ifft(tf1);

figure

stem(ir_from_tfestimate)

hold on

stem(ir_from_fft)

set(gca, 'XScale', 'log')

- So I am wondering what is the nature of this tale? Maybe there is a mistake in my code? I was google and looking for this problem/explanation, but did not find anything
- Can it be related how I define nfft and what is the right choice for number of points? Can we say if nfft is defined not correctly, impulse response will not be correct?
- When I am trying to recover symmetric part of a spectrum obtained from `tfestimate`, am I doing it correctly: ?

tf = tfestimate(x, y, hann(nfft), nfft/2, nfft, 48000);

if mod(length(tf),2)==0 % iseven

tf_sym = conj(tf(end:-1:2));

else

tf_sym = conj(tf(end-1:-1:2));

end

Thank you!

David Goodmanson
on 13 Jan 2021

Edited: David Goodmanson
on 13 Jan 2021

Hi Yurii,

It looks like the 'lowpass' function is doing something unexpected, which causes this. In order to get the transfer function (tf) you are using a chirp input into lowpass and dividing ffts, and then of course doing the ifft to get the impulse response (ir) in the time domain.

If you zoom in on the resulting ir you will find that the 'tail' is actually a perfect time-reversed copy of the initial pulse. Since the fft assumes the functions are periodic, one can swap halves of the array with fftshift and go with a shifted time array that puts t=0 at the center. That is similar to what is often done in the frequency domain, and that is done here as well.

The first part of the code below takes the transfer function you got from chirp (and ignores the tf_estimate method). Plot 1 shows the tf and the resulting ir, expanded a bit with the use of xlim. Note that the tf is a symmetric real function, plus a very small imaginary part. I did not expect this at all. Since [symmetric and real] fourier transforms to [symmetric and real], the resulting ir is symmetric about t=0 and has a lot of extent into negative time. This means that the filter is not causal.

In the fine print for the documentation of lowpass, it says that it "compensates for the delay introduced by the filter" so I am guessing that they multiply the filter response by exp(i*w*tau) (or perhaps someting oscillatory but more complicated) in order to perform a time shift. I suppose they have their reasons, maybe good ones. But the bigger the company gets, and the more convenience functions and black boxes they introduce, the more opaque things become. That's just the way it is.

The second part of the code uses a simple lowpass butterworth filter as the tf. The frerquency response is similar to the lowpass filter you used, but I did not attempt a strict AB match. In this case the tf has as much imaginary part as it does real part, which allows the filter to be causal, and you only have an ir for t>=0.

You have an odd number of points, which is very good for accuracy in this calculation, but it does mean that fftshift and ifftshift are different functions so you have to be careful which one to use.

chirp

butterworth

t = 0:1/48000:10;

fo = 1;

f1 = 24000;

x = chirp(t,fo,10,f1,'logarithmic');

x = x(:);

y = lowpass(x,1000,48000);

y = y(:);

in = fftshift(fft(x));

out = fftshift(fft(y));

tf = out./in;

ir1 = fftshift(ifft(ifftshift(tf)));

% create arrays with f=0 at the center, t=0 at the center

n = length(t);

fcen = 48000*(ceil(-n/2):ceil(n/2)-1)/n; % f=0 at the center

tcen = t - (t(1)+t(end))/2; % t=0 at the center

figure(1)

subplot(2,1,1)

plot(fcen,real(tf),fcen,imag(tf));grid on

xlim([-20000 20000])

subplot(2,1,2)

plot(tcen,ir1); grid on

xlim([-.004 .004])

% butterworth filter

n = 5;

f = 1000;

[zb,pb,kb] = butter(n,2*pi*f,'s');

[bb,ab] = zp2tf(zb,pb,kb);

[hb,wb] = freqs(bb,ab,fcen*2*pi);

ir2 = fftshift(ifft(ifftshift(hb)));

figure(2)

subplot(2,1,1); grid on

plot(fcen,real(hb),fcen,imag(hb))

xlim([-1e4 1e4])

subplot(2,1,2)

plot(tcen,ir2); grid on

xlim([-.006 .006])

David Goodmanson
on 16 Jan 2021

Hi Yurii

For function names It's not that easy to find letters that do not have other connotations. Here I will go with p(t) and q(f),

[1] The fft is definitely not the same as the ordinary fourier transform. In that transform, both t and f run from -inf to inf, and p(t) and its transform q(f) are both continuous functions. Next in line is the transform from freshman physics and engineering (I imagine you've been through that before), where time goes from 0 to T and p(t) is continuous and periodic with period T. That is all that's happening here.

Because of the periodicity, the transform q(f) is an infinite sum of exponentials at evenly spaced discrete frequencies. Since p(-t) = p(T-t), for small negative t the signal ends up at the top end of the array. Also since p(t) is periodic, you can slide the signal to the left or right in the window and have the same information, as is done with fftshift and ifftshift.

Next in line is the fft, with both t and f being discrete and evenly spaced, each one periodic from 0 to T or 0 to F, each with N points.

[2] In the frequency domain, the transfer function A(w) of a causal filter obeys the Kramers-Kronig relations (see e.g. WIkipedia), by which Im(A(w)) equals an integral involving Re(A(w)), and vice versa. So if A is real only, it cannot be causal. And if A is causal, then the K-K relations should work. But one of the easier ways to show causality is to just fft the transfer function and see how it comes out for negative t, as with the plots in the answer.

In the causal example I used, I had the advantage of an analytic function (butterworth). IF you had a causal black box, then by getting A(w) from dividing output by input, things are not so perfect and a few points with small amplitude might straggle into negative time. That would not be a big deal.

[3] Consider the transform from p(t) to q(f).

'f=0' refers to the point in the q array that corresponds to f=0, not the f array itself.

For the transform q = fft(p), f=0 is at point 1. fftshift puts f=0 up at the center of the array, and ifftshift does the reverse, put the center point down to point 1. But for even n, the array does not have a true center point. In that case both fftshift and ifftshift swap the two halves of the array. fftshift puts point 1 up to point n/2+1, ifftshift puts point n/2+1 back down to point 1. For even n, fftshift and ifftshift are their own inverses.

For odd n, swapping halves obviously not work, but odd n has the advantage of having a true center point, (n+1)/2. fttshift puts point 1 up to (n+1)/2 and ifftshift puts it back down to point 1. For odd n. fftshift and ifftshift are inverses of each other, but they are not their own inverses.

Since the fft is basically symmetric in f and t, everthing said above is the same for t <--> f, p <--> q

[4] on your most recent plots, TBD

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 1 Comment

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/714033-why-when-calculating-impulse-response-of-a-filter-system-i-get-a-tail-at-the-end#comment_1257408

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/714033-why-when-calculating-impulse-response-of-a-filter-system-i-get-a-tail-at-the-end#comment_1257408

Sign in to comment.