MATLAB Answers

Why when calculating impulse response of a filter (system) I get a “tail” at the end

95 views (last 30 days)
Yurii Iotov
Yurii Iotov on 11 Jan 2021
Edited: David Goodmanson on 17 Jan 2021
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')
  1. 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
  2. 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?
  3. 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!

  1 Comment

Mathieu NOE
Mathieu NOE on 11 Jan 2021
hi
could not really distinguish the 2 plots with stem , so I use the regular plot - and both traces show the tail
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 = real(ifft([tf; tf_sym]));
in = fft(x, nfft);
out = fft(y, nfft);
tf1 = out./in;
ir_from_fft = ifft(tf1);
figure
plot(ir_from_tfestimate,'*b')
hold on
plot(ir_from_fft,'or')
set(gca, 'XScale', 'log')
hold off

Sign in to comment.

Answers (1)

David Goodmanson
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])

  2 Comments

Yurii Iotov
Yurii Iotov on 14 Jan 2021
Hi David
First off all thank you very much for your comment and explanation.
I was recalling and analyzing what you explained for the past day and have some comments and questions, let me please ask you.
1. Regarding your 2nd paragraf, the tail indeed is a perfect time-reversed copy of the initial pulse, and you wrote: 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.
What gives us the right to do this operation in time domain, is there any property? Because I do understand why from a real function whe get from its fft also conj symmetric copy, but then going back to time domain with ifft, this symetry vanishes (should)? So from physical or math meaning I do not understand this point, please explain.
2. ...Note that the tf is a symmetric real function, plus a very small imaginary part. ... This means that the filter is not causal.
...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.
Is there any property or it's from personal experience regarding the distribution of real and imaginary part and how it corresponds to the definition that filter is causal or not? For example if we have a black box and we use pink noise as a input, calculating a TF and then looking at real and imaginary part can give us an indication that filter might not be causal?
3. 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.
I do understand that these are different functions, but what do you mean in relation to number of points and these functions use?
And coming back to my real measurement and signals, I would like to get you opinion: I have measured time signals (fs = 48000), as an input of a system I used pink noise. The system 1 is sound propagation trought a small duct (actually it has a bit more complex shape).
The system 2 is the same but involves also a small loudspeaker (I imagine it as loudspeaker(LS) - 1st filter, propagation from LS to another mic is a 2nd filter). I attached mat files as well.
Then I did TF and impulse response calculation in the same way, and below what I have got from the following code:
load('Measurements.mat') % get mat file from here,
% since its more than 5 mb and I couldnot attach here https://dropmefiles.com/809zm
nfft = 4*1024;
% nfft = 1200000;
win = hann(nfft);
% win = [];
fs = 48000;
[TF1, freq_vector] = tfestimate(System1_IN, System1_OUT, win, nfft/2, nfft, fs, 'twosided');
[TF2, freq_vector] = tfestimate(System2_IN, System2_OUT, win, nfft/2, nfft, fs, 'twosided');
ImpResp_tfestimate1 = (ifft([TF1]));
ImpResp_tfestimate2 = (ifft([TF2]));
t = 0:1/48000: length(ImpResp_tfestimate1)/48000 - 1/48000;
fcen = 48000*(ceil(-length(t)/2):ceil(length(t)/2)-1)/length(t); % f=0 at the center
tcen = t - (t(1)+t(end))/2;
figure('Name','tfestimate','NumberTitle','on');
tiledlayout('flow')
nexttile;
plot(fcen,real(TF1),fcen,imag(TF1))
title('System 1')
nexttile;
plot(fcen,real(TF2),fcen,imag(TF2))
title('System 2')
nexttile;
plot(tcen, fftshift(ImpResp_tfestimate1))
title('System 1')
xlim([-.001 .002])
grid
nexttile;
plot(tcen, fftshift(ImpResp_tfestimate2))
title('System 2')
xlim([-.002 .002])
grid
%%
IN1 = ( fft(System1_IN) );
OUT1 = ( fft(System1_OUT) );
TF_fft1 = OUT1./IN1;
ImpResp_fft1 = (ifft((TF_fft1)));
IN2 = ( fft(System2_IN) );
OUT2 = ( fft(System2_OUT) );
TF_fft2 = OUT2./IN2;
ImpResp_fft2 = (ifft((TF_fft2)));
t = 0:1/fs: length(ImpResp_fft1)/fs - 1/fs;
fcen_fft = fs*(ceil(-length(t)/2):ceil(length(t)/2)-1)/length(t); % f=0 at the center
tcen_fft = t - (t(1)+t(end))/2;
figure('Name','fft','NumberTitle','on');
tiledlayout('flow')
nexttile;
plot(fcen_fft,real(TF_fft1),fcen_fft,imag(TF_fft1))
title('System 1')
nexttile;
plot(fcen_fft,real(TF_fft2),fcen_fft,imag(TF_fft2))
title('System 2')
nexttile;
plot(tcen_fft, fftshift(ImpResp_fft1))
title('System 1')
% xlim([-.001 .002])
grid
nexttile;
plot(tcen_fft, fftshift(ImpResp_fft2))
title('System 2')
% xlim([-.002 .002])
grid
% sound(filter(ImpResp_fft1,1, yoursig), 8000)
% sound(filter(ImpResp_tfestimate1,1, yoursig), 8000)
Looking at these figures, and taking into account real and imaginary part of TF, could you please tell me what might be the meaning of the impulses before 0? Since these systems is a propagation and loudspeaker involved I do not expect any non-causal filter (from physical point of view), but there might be propagation time (delay) involved, but I am not sure how to see it/measure correctly from the TF and impulse response.
Thank you! I would be gratefull and appreciate your help.
David Goodmanson
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

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by