Frequency Response Functions (FRF) in Julia

Yes, the FFT is sometimes presented in this way, but an FRF never. Where have you seen an FRF with both negative and positive frequencies, if I may ask?

Edit: Sorry, I just realized I misread the post to which I was responding. The graph was correctly labeled as FFT, not FRF. I withdraw my comment.

A couple of comments. The fft algorithm gives as many frequency points as time points. You display all frequency points, both the positive frequencies and the negative frequencies. Plotting the negative frequencies does not add information and takes up plotting space. So take your frequency array, truncate it to half length, and multiply the amplitude by 2 to retain all the energy from the negative frequencies. Most analyzer companies will keep the frequency points length(Gyy)/2.56, to eliminate the points that are potentially affected by aliasing, depending upon the filtering of the time data.

Rafael has the calculation of the frequency bin spacing correct, which in the example is 1 Hz. The simplest relationship is

Δf = 1/t = f0/N    # where t is the duration of the block of data analyzed

His point about physical frequencies is that it has units of [Hz] which means cycles/second, or revolutions/second when dealing with rotating machines, not radians/second.

The amount to pretrigger is usually in the order of 4 - 10 samples. A number of reasons for this; It allows more time for the signal to decay and have the complete transient in the block of data that is being analyzed. It will include the samples that are influenced by the impact, prior to the trigger. Some will use the average amplitude of these pretrigger points for mean removal of the time data prior to the fft.

Compare the values of

mean(Hammer)
mean(Hammer[1:6])  # where the pretrigger is set to 7 samples

Because it is a transient, I expect that there will be a substantial difference and the later will be more accurate.

With what you are doing with your data, not using windows is perfectly acceptable. If your acceleration signal is completely decayed by the end of the block of data going to the FFT, then it is 0 from -inf in time to the beginning of the block and again 0 from the end of the block to +inf in time, thus making it a self-windowed function.

The reason for using the exponential window is in case the data has not completely delayed by the end of the block of data. The second reason is if there is some accelerometer cable noise, the accelerometer signal is a fairly low amplitude by the end of the block, so it becomes a low signal to noise ratio. The exponential window can help alleviate this. Again for where you are in this project, it is not essential.

Blockquote
How do I apply a force window to the input data, and how do I apply an exponential window to the response data?
I am familiar with DSP and applying filters to signals, but not a specific waveform as depicted for force windows.

To apply a window, take the block of time data that you are working on and multiply it, element by element, with the window desired. When generating an FRF, both the force and response data is multiplied by the same window and thus any window correction factor will cancel out.

How does Gxy differ from Gyx?

Gyy = Sy .* conj(Sy)
Gxy = Sx .* conj(Sy)
Gyx = Sy .* conj(Sx)

The pattern is that the first subscript is the Spectrum associated with that subscript and it is multiplied by the complex conjugate of the spectrum associated with the second subscript.

I have coded windows in Julia and am willing to donate them to the DSP package. I currently don’t have time to do that myself, so if there is a volunteer out there, let me know.

FRFs are usually only displayed with positive frequencies due to symmetry, but some systems does not have this symmetry and you have to display the negative frequencies as well. This kind of systems commonly arise in control of accelerators and power systems, a friend of mine, @olof3 , has written many nice articles on the subject, for example

Thank you, that would be interesting to explore. The paper that you linked to
does not show that kind of FRF. Perhaps you would know where to look?

The positivity of the frequencies in FRF graphs has little to do with symmetries, but rather with the physical meaning of the forcing. At least that’s what I am familiar with from mechanical systems. I’d be happy to expand my horizons.

Many thanks,

Found this: Negative Frequency - an overview | ScienceDirect Topics
Not a mechanical system…

1 Like

No, but let me give some background. I am using this exercise to learn Julia, Pluto, DSP, and specifically modal analysis. I forgot to reset both accelerometers to 200 Hz, so one booted up to its default of 256 Hz, and the other to 100 Hz (the one I put on the hammer for the input force or excitation force). I have long data sets that I was using, but thought you needed to window or capture only the event. 256 samples allowed for a lead-in and a fade out of the output response visually, so I clipped the data to that bucket.

I will provide the raw data and give a cheat sheet of what was what in the next 24 hours for you to play with. I purchased Peter Avatabile’s book Modal Testing per @PetrKryslUCSD suggestion, and it is a great book for me. Practical and enough theory to hurt myself! I would hope I could use any simulations you run as a sanity check of my efforts. I am sure I missed the mark, but I sure learned a lot! Special thanks to @Jake for his lengthy and detailed replies!

2 Likes

@Jake thanks again!
It all makes sense, and between you and Peter’s book I am starting to grok it all.

I am going to spend some time tonight with different parts of the data I captured and see how the graphs change in the interactive Pluto notebook. I am also going to play with windows and cutting out the right portion of the FRF’s symmetry and multiplying by 2. I had to comment out:

newfft = 2 .* [1:floor(Int, length(newfft)/2)]

because I had type errors or local vs. global scope (why I started naming things differently in some of my begin/end blocks). I will figure these out on the Julia side. I also didn’t utilize the ‘bs’ (bin size) variable and portion you had, since this particular exercise had a 1-second sample at 256 Hz, so it was 1 Hz per sample on the Freq. axis. I may carry out this exercise on a local steel structure in the playground across the street to get a few more things right and get some practice with my methodology. I would love to play with a real instrumented hammer/DAQ deck, but the software and all-in kits are well over $15k USD. I have one expensive accelerometer from Lord Microstrain that is getting old, and several cheap ones from Witmotion. I am concerned that the Witmotion sensors have too much built-in filtering for true, raw data. I believe they are geared towards app/drone markets where people do not need high frequency sample rates and want smooth pretty graphs and data to work with when integrating them into their systems.

On a personal note, I am championing a DAQ department at my engineering company. We have strain gauges that we purchased with the proper DAQ box and wireless nodes, and it proved to be very helpful. Our engineers had done a very complicated analysis of a steel structure in dynamic loading, and they wanted confidence in their calculations. I gathered 16 strain gauge readings in real time that showed how the structure played out during the variational dynamic loading. It was a real appetizer to do more of this work. I love math, programming, engineering, and actually working in the field, so a win-win for me. Thanks for all of your help!

1 Like

For linear systems with real coefficients (certainy mechanical systems, but most other systems as well), the frequency-response functions is symmetric w.r.t. the zero frequency so the negative frequencies, while there, are usually not shown. If the system has complex coefficients the FRF no longer exhibit this symmetry. The article I linked does have those systems, but displayed the system as a Nyquist plot instead of a Bode plot. The Nyquist plot usually exhibits symmetry w.r.t. the real axis, but not for comple-coefficient systems. Page 94 of the thesis

will show you the Bode curves of such systems, representing parts of the accelerator at the European spallation source.

1 Like

@baggepinnen, the first article you linked shows a Bode plot in Figure 2, where the asymmetry can be seen. Please note the subtle display, which uses dashed lines for ω < 0.

2 Likes

I don’t know of that many mechanical examples, but here are a few: Active control of magnetic bearings. MEMS gyroscopes. The Foucault pendulum (probably the most classical example, mentioned in, e.g., Arnold’s book).

Complex-signal models and transfer functions are essential for telecommunications, but there are not many applications outside that field. Although when a complex-signal modeling is applicable, it usually gives a lot more insight.

3 Likes

In addition to what @olof3 said, it is quite common to use two-sided spectra when studying rotating systems (e.g. turbines, compressors …) because due to the symmetry properties of the structure, one generally (at least numerically, and experimentally in the right conditions) obtains propagative modeshapes.
Due to the rotation, the frequency of the forward and backward propagating modeshapes are different when you observe your structure from another reference frame (for instance, observe a rotor from a stationary frame) or for other reasons as Gyroscopic effects.
These can be shown with a single sided spectrum, but with a two sided spectrum (either FRF, spectrogram or else), you can actually make a difference between the propagation directions.

3 Likes

I confess to becoming a bit confused here. Bentley Neveda pioneered full spectrum analysis, where the x & y components of vibration at a bearing are shown on the same spectrum. Is this what you are referring to @BambOoxX? That is different from the positive and negative parts being different.

Now back to @Eggy. Proper data acquisition for modal analysis would be simultaneous sampling of all transducers. This implies sampling at exactly the same rate for all signals. And all signals are treated identically. With your data sampled at different rates, your data is not good. For learning or for data that cannot be collected again, you could try and digitally filter and resample the high sample rate data to get it at the same sample rate as the low sample rate data. This will at least allow you to create cross power spectra. However since it is not sampled at the same time your phase is likely to be bogus, but at least a consistent delay time between force and response. If your DUT (device under test) is convenient it is much better to collect the data a second time. For your purposes, using the hammer with accelerometer is fine. Using F=ma, your acceleration is proportional to Force. There is a force hammer calibration technique that uses this equation. See slide 13 of https://www.modalshop.com/filelibrary/Modal%20Excitation%20Tutorial.pdf.

For a nice cost effective DAQ for DSP, take a look at IEPE Measurement DAQ HAT for Raspberry Pi from Measurement Computing. This is a $200 US per channel solution, which is the lowest cost solution that I am aware of. I have played a bit with it on an 8 Gbyte Raspberry PI, and put most of the c function calls into Julia. If someone wants to try it perhaps we can make a Julia package to interface to it directly.

1 Like

Good catch!

1 Like

I really like your thesis! Especially that it has lots of pictures. (As the saw goes, one picture is worth a thousand words, but it could be extended to equations as well :slight_smile: .)
I did not follow the Foucault pendulum example reference though. The system is unforced. What FRF could be constructed?

2 Likes

I am sorry if I have hijacked this thread. I do enjoy the discussion though. It is so good to have one’s horizons expanded by interaction with researchers in other subject areas! Thanks for the resources you all compiled!

3 Likes

No, I am not referring to direct bearing measurements, though we also talk about full spectrum analyses. I was more considering bladed disks.

No you’re right, I forgot that was the question. But it is a nice and classic example of complex-valued modeling though :slight_smile:

1 Like

@Jake It’s still not right, but the original data for the hammer accelerometer was actually 200 Hz, not 100 Hz, so It is 256 Hz on the structure, and 200 Hz on the hammer now. Somehow saving it as CVS from the accelerometer’s software downsampled it to 100 Hz. Fortunately I have the raw data at 200 Hz. At this point, I really need to know the natural frequency of the structure and nothing else. Perhaps I can use an OMA technique with the knowledge of the hammer’s force? I realize phase will still be invalid with resamping 200 Hz to 256 Hz, but maybe another technique will help here? I will have a second bout, but I need to work with this for now. It’s been a great learning experience.

Are you recommending downsampling the 256 Hz data to 200 Hz instead of the hammer single impulse to 256 Hz? I feel like the response data at 256 Hz is more complex than the hammer impulse and may affect the calculation more negatively than the change to the hammer impulse. I realize it is best to downsample than upsample here, but asking anyway.

Also one important point I forgot: If I have 2 seconds of data at 200 Hz does it matter if my first vector is the sample number vs. the accumulated time? In other words, I have been using the Z-axis accelerometer, single vector to perform my FFT and other calculations and it defaults to the sample number on the X-axis vs. the Z-axis (g) on the Y-axis of my plots. Does the process require two vectors (or 2 column array), or my single vector approach is OK?

I plan on restaging the test on the same structure, but not for a while (~1 month from now). I will try and put together a system and practice with it before then.

The DAQ HAT for Raspberry Pi looks promising. Can the Raspberry Pi handle two or more of these things sampling at 51.2 kS/s - a 4 channel system? The Lord wireless system allows for the connection of many accelerometers, and it is relatively inexpensive, but the accelerometers are $500 plus. They are wireless, have a lossless protocol, and their SensorConnect software is opensource .NET/C#. Again, I was pricing out a DAQ system, accelerometers, and impact hammer. Over $12k! Renting is pretty pricey too. @Jake Does your research or job involve DSP and vibration analysis?

If you take a look at slide 15 of https://www.modalshop.com/filelibrary/Excitation_Techniques_IMAC_30_by_Pete_Avitable.pdf, you will see that a good impact excitation force is approximately equal across frequency bins. Your plot of Gxx has a dip at sample ~10 which we concluded is 10 Hz. Your response has a peak at about the same frequency. From this I would conclude that the excitation force was not good, and that there is a resonance at ~10 Hz. If your frequency spacing of the response and force signals are different, which they will be due to the different sampling rates, you cannot obtain valid cross spectra, so just view the response spectra.

The Raspberry PI can handle up to 8 of these HATS. Its bus limits throughput, so by recollection it can handle 3 channels at 51.2 kHz, and the maximum sampling rate decreases with the number of channels. For modal you are not likely to need these high throughputs, perhaps limiting for envelop analysis of bearings at higher channel counts.

I do a lot of vibration analysis, including modal from time to time.

1 Like

@Jake Thanks again!
Just to be clear for me: the 10 Hz is not the natural frequency but one of a few. When you say “just view the response spectra”, do you mean just the FFT of the response accelerometer data, or Gyy?
In that case it is around 10 to 12 Hz. However, what is the line to infinity on both sides of the “Output Power Spectra” of Gyy?

I am going to definitely check out the Raspberry Pi HATS. I have gotten spoiled by the wireless unit from Lord Microstrain, however, adding more sensors quickly adds up.

Thanks!