Computing Impulse Response Function from FRF

[SrinivasAluri]

I am trying to compute an impulse response function from a one-sided Frequency Response Function. How do I do that in Matlab?

 

[GregLocock]

IFFT the magnitude of the spectrum.



Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[VisiGoth]

Make sure that what you asked for is what you want. GregLocock is correct, you will get "an" impulse response, just as you requested. But are you sure you want the minimum phase version or do you ned the "correct" version?

 

[GregLocock]

Well OK, I'm interested.

Ah OK, I've got it. If you zero the phase and IFFT the magnitude of the spectrum, then by definition when you re-FFT the resulting impulse there will be no phase information.

What's the 'correct' version, in signal processing terms?











Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[MikeyP]

It can be proved that the FRF and IRF are a Fourier transform pair. You need to take into account both the magnitude and the phase of the FRF in performing the calculation.

In Matlab the steps to compute the IRF are as follows

1) Express the FRF in complex form: A*cos(theta) + i*A*sin(theta), where A is the magnitude and theta is the phase at each individual spectral line.

2) Convert the complex spectrum to the two-sided form required by Matlab. i.e. arrange the spectral lines in the frequency order:
[0 df 2df ... (fn-2df) (fn-df) fn -(fn-df) -(fn-2df) ... -2df -df]
Where df is the frequency spacing and fn is the Nyquist frequency. The complex values at the negative frequencies are the conjugate of the complex values at the associated positive frequencies.

3) Perform the inverse Fourier transform.

4) Check that you have done it correctly, by making sure that the IRF is purely real (there may be some VERY small - 15 orders of magniude smaller - imaginary part which is due to numerical errors and can safely be removed by taking the real part)

5) You may need to apply a scaling factor based on the length of the IRF divided by 2 and possibly a further factor of 2 depending on exactly what you mean by a "one-sided" spectrum.

M

--
Dr Michael F Platten

 

[MikeyP]

There can be only one impulse response function. That is the response of the system to an ideal unit Dirac impulse. A "minimum phase" impulse function is meaningless.

M

--
Dr Michael F Platten

 

[SrinivasAluri]

Thanks for the response. The FRF was computed using LabVIEW 'Transfer Function' function, so I am assuming its properly scaled.

 

[SrinivasAluri]

If there are 8192 spectral lines in the FRF
the array specified by Dr.Platten would have 8192 (including positive and 0 Hz) plus 8191 (negative spectral lines )? Is this correct?

 

[VisiGoth]

MikeyP has given the detailed instructions on how to get the one and only "correct" response. The minimum phase "might" be correct. But since you seem to have a complex FRF single sided I assume the data you want is real. If so the question you have is how to do the MikeyP step 2. Note the instruction to contruct the two sided. You can take your one sided and make up the two sided by knowing that the values of the spectral lines of the negative frequencies are complex conjugates of positive frequencies. You can look at fliplr and conj and concatenating to create an array of twice the size.

A hint on step three is to check the enegy in the time domain and in the freq domain and find the scale factor to make them match. x * x' should do it for a column vector.

 

[SomptingGuy]

You need to reconstruct the double-sided spectrum in the order that Matlab likes it. Here are the pitfalls:

1) Scale your spectrum properly. The absolute value in each bin should be proportional to the number of points in the spectrum.

2) Matlab stores its spectrum starting from DC, going up to nyquist and then the negative frequencies back up to DC. So to reconstruct the spectrum you need to add the flipped conjugate to the end.

Having got the properly scaled double-sided spectrum, take the real part of the ifft().

 

[SomptingGuy]

See also, invfreqz(). If you set NA to zero you'll only get B coefficients - the impulse response of the input spectrum.

 

[MikeyP]

LabView has probably thrown away the value at the Nyquist frequency (this is common practice), so you only have values from 0 Hz to (fn-df) Hz. Just set the value at the Nyquist frequency to zero. That will give you a total of 16384 points.

M

--
Dr Michael F Platten

 

[SrinivasAluri]

Here is how I wrote the code:

frf is my frf row vector with 8192 points.
FRF=[frf,fliplr(conj(frf))];
irf=ifft(FRF);
irf=irf+conj(irf);
plot(irf)

But there is a problem with the plot, because it shows a exponentially decayed time signal at the beginning of x-axis and also a reversed exponentially decayed time signal at the end of x-axis. Am I doing something wrong?

 

[VisiGoth]

Per MikeyP's advice try
FRF=[ frf 0 fliplr( conj( frf( 2 : end ) ) ];

then replace
irf=irf+conj(irf);
with
irf = real( irf );

Could you give us the first 4 complex numbers in the frf?
They might give us some insight.

 

[SrinivasAluri]

-7.19579E-16+8.81201E-32i
-8.32218E-3+1.33314E-3i
-3.37827E-2+2.03763E-2i
-4.38152E-2+7.85120E-2i

 

[VisiGoth]

It looks like you have no DC, so you might want to force the first data point to zero, or at least make the imaginary part exactly zero.

It looks like a phase roll of about 10 degrees per picket so I would guess that you should have your impulse spike (assuming you have a basically flat freq response over most of the band) at about the 5th or 6th time data point.

Is the freq response flat to within a couple dB over 80% of the band?

Or is this some sort of high Q resonant problem? in which case my back of the envelope analysis is inadequate.

So, can you give us the first 14 time data values?

 

[VisiGoth]

At what picket does the frf flatten out? I assumed with a nice long MLSSA we would have seen a max in the frf. So, maybe your system frf actually is HP or we are seeing MLSSA (or whatever is used) artifacts.

Can you give a snippet of a few points around where the HP edge has ended?

 

[SrinivasAluri]

Thanks for the efforts. FYI: the data is obtained from a vibration test of a composite deck. The following is the first 14 points of the FRF.

-0.00000000000000 + 0.00000000000000i
-0.00832218000000 + 0.00133314000000i
-0.03378270000000 + 0.02037630000000i
-0.04381520000000 + 0.07851200000000i
0.00500428000000 + 0.13423600000000i
0.06934280000000 + 0.12662800000000i
0.08197070000000 + 0.09926960000000i
0.09618740000000 + 0.11477100000000i
0.16896000000000 + 0.09298500000000i
0.17899000000000 - 0.02864460000000i
0.04210450000000 - 0.08325290000000i
-0.03197860000000 + 0.04724550000000i
0.11644800000000 + 0.11310100000000i
0.19826300000000 - 0.09455950000000i

 

[GregLocock]

Please ignore my first post and all subsequent discussion referring to it. My first post was completely in error, I don't know where that idea came from.

If you plot the magnitude of that lot you'll see that there is probably far more information to come in that spectrum.




Cheers

Greg Locock

Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[MikeyP]

Try this:

FRF_2side = [0 frf(2:end) 0 conj(frf(end:-1:2))];
irf = ifft(FRF_2side);
plot(real(irf));

Does the plot look OK (an oscillating decaying signal)?

Is max(abs(imag(irf))) very small (<1e-10)?

If the answer to both questions is "yes" then all you need to do is apply a scale factor to the irf (either multiply or divide by 8192 - I never can remember which).

If the answer to either of them is "no" then please tell me the sampling frequency in Hz of the original measurement and the frequency in Hz of the 1st, 2nd and last spectral lines in the FRF.

M


--
Dr Michael F Platten