×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Log In

Come Join Us!

Are you an
Engineering professional?
Join Eng-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!
  • Students Click Here

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Jobs

Computing Impulse Response Function from FRF

Computing Impulse Response Function from FRF

Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

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.

RE: Computing Impulse Response Function from FRF

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?

RE: Computing Impulse Response Function from FRF

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.

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

(OP)
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?

RE: Computing Impulse Response Function from FRF

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.

RE: Computing Impulse Response Function from FRF

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().

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

(OP)
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?

RE: Computing Impulse Response Function from FRF

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.

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

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?

RE: Computing Impulse Response Function from FRF

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?

RE: Computing Impulse Response Function from FRF

(OP)
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

RE: Computing Impulse Response Function from FRF

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.

RE: Computing Impulse Response Function from FRF

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

RE: Computing Impulse Response Function from FRF

(OP)
Dr. Platten: Your code works except that the IRF is flipped!! Sampling rate for the signal is 1280Hz and DeltaF is 0.078125 Hz..

This IRF business is more complicated that I thought..

RE: Computing Impulse Response Function from FRF

Hmm. It sounds like there may be something wrong with your original FRF.

Does the IRF decay initially and then grow towards the end or does it start from zero and grow continuously?

I would expect some non-causal response at the far end of the IRF. This is normal, but I would expect the signal to decay for at least the first half of the IRF before it started to grow again.

How was the FRF measured? Random excitation? Hanning window? Number of averages? Averaging Overlap?

Is the structure very lightly damped? What is the typical critical damping ratio for the modes?

M

--
Dr Michael F Platten

RE: Computing Impulse Response Function from FRF

(OP)
I really appreciate all the responses. I put up the data at the following link http://www.cnuga.com/119-1.html  it is in real imaginary format, first column being real.

The image of IRF I got from using Dr. Platten's code can be found at http://www.cnuga.com/irf.html

The FRF was measured using impact hammer - rectangular window was applied on impact, no window on time response. Number of averages is 4. The data was zero-padded before evaluating the FRF.

RE: Computing Impulse Response Function from FRF

Your frf is showing wild multiple resonances.  So would not get anything that looks like a modified impulse (dirac).  However, it looks like you have the result of an analysed system and the model went off to lala land.

If you have strange resonant system you are looking at, then you would want to carefully excite it with a stepped sine system, changing the input drive as a funcion of the resonant gain, so as to stay away from non-linear problems.

Now, if your first column was dB and your second degrees, then you might have something.

RE: Computing Impulse Response Function from FRF

VisiGoth,
The data looks fine to me. The columns are the real and imaginary parts of the FRF as the original poster stated.

SrinivasAluri,
I think you must still have a mistake in your code somewhere. The results I got look OK. There is a very small non-causal response at the end of the IRF, but I have seen much worse on my own data.

Here is the code that I used

CODE

load frf_et.txt % load the text file containing the data

frf = frf_et(:,1)+i*frf_et(:,2); % form the complex FRF
fn = 640; % Nyquist frequency
df = fn/8192; % frequency spacing
freq = [0:df:fn-df]; % frequency vector
frf2 = [0 frf(2:end).' 0 frf(end:-1:2)']; % form 2-sided FRF
irf = ifft(frf2); % calculate IRF
dt = 1./2./fn; % sampling interval
time = 0:dt:16383*dt; % vector of sampling times
% plot the results
figure;
subplot(3,1,1);
plot(freq,abs(frf));
xlabel('Frequency (Hz)');
ylabel('FRF Magnitude');
subplot(3,1,2);
plot(freq,angle(frf));
xlabel('Frequency (Hz)');
ylabel('FRF Phase (rad)');
subplot(3,1,3);
plot(time,real(irf));
xlabel('Time (s)');
ylabel('IRF');

M

--
Dr Michael F Platten

RE: Computing Impulse Response Function from FRF

The only difference between this code and the code I posted previously is that I assumed in the earlier post that that the FRF data was in a row vector. The FRF data in my last post is in a column.

M

--
Dr Michael F Platten

RE: Computing Impulse Response Function from FRF

(OP)
Thank you Dr. Platten for the code. I wanted to generate this IRF so that I could write a program fof least square complex exponential method and get a feel for it.

RE: Computing Impulse Response Function from FRF

In that case, the non-causal blip at the end of the IRF shouldn't cause you any problems as you will only need to use the start of the decay in the LSCE fit anyway. If you plot the magnitude of the IRF on a log scale, it is easy to see where the decay disappears into the noise floor of the measured data. You should only fit the part of the IRF data up to that point.

Also, if you are just using LSCE to get frequency and damping values, I wouldn't worry about the scaling factors discussed in earlier posts. The scaling will make no difference to the estimated frequency and damping values.

Finally, I would add that the FRF data you posted had rather a lot of modes in your frequency range. This means you will probably have to do the LSCE fit in smaller frequency ranges. In that case, you should simply band-pass filter the FRF prior to performing the IFFT. Simply set the unwanted spectral lines to zero.

Good luck

M

--
Dr Michael F Platten

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Eng-Tips forums is a member-only feature.

Click Here to join Eng-Tips and talk with other members!


Resources