Frequency Response Functions (FRF) in Julia

I have followed Julia for a while, but I am now just getting into using it more. I have collected data from performing what’s known as an roving hammer test with my accelerometer and an non-instrumented hammer to derive the natural frequency of a structure. Typically you use a pre-configured DAQ, software, hammer system to do this, but this was a learning exercise. I am curious if there exists vibration analysis or other libraries for FRFs (Frequency Response Functions) such as this link shows for MATLAB: Frequency-response functions for modal analysis - MATLAB modalfrf (mathworks.com)

For the last part of my experiment, I placed an accelerometer on the hammer to emulate an instrumented hammer used for such vibration analysis exercises. An FRF can use the input (the accel data from the hammer) and the output (the accel data from the column’s accelerometer) to derive the natural frequency of the structure. I played with the Pluto notebook on FFT from the YouTube video:

What is a Discrete Fourier Transform? | Week 14 | MIT 18.S191 Fall 2020 | Grant Sanderson

I am a vibration novice, so I am sure there is a bridge from simple FFT analysis to the more involved FRF hammer test analysis.

I am loving the interactive nature of Pluto in exploring new subjects. A great learning tool for a self-learner .

Thanks!

1 Like

There is a nice DSP library but alas not FRF’s. However they are not that difficult to write, which is what I have done. The normal procedure is

  1. Find a trigger point on your force spectrum.
  2. Pretrigger by a number of samples
  3. Take the number of samples you want for the required frequency bin spacing.
  4. Apply an exponential window to all the channels
  5. Apply a force window to the impact channel only.
  6. Perform a spectrum on each channel.
  7. Calculate the auto power and cross power spectra.
  8. Calculate H1, H2, and coherance as appropriate. Hv and Hs are more involved.

You should be able to find the equations by doing a search on H1 FRF. Happy computing.

1 Like

Jake,

I have my googling cut out for me!
So I see the initial hit (high g) on my hammer strike. By ‘spectrum’ do you mean the g vs. time waveform? Or do I need to do an FFT on it first?
I appreciate you laying out the steps.
Can you recommend a good book or paper on doing this process without preconfigured software or equipment? Thanks!

Check out the book by Peter Avitabile.

1 Like

I just ordered it from Amazon after previewing it. It should help me understand this from first principles, and work through Jake’s steps above. Thank you!

GitHub - baggepinnen/ControlSystemIdentification.jl: System Identification toolbox for LTI systems, compatible with ControlSystems.jl has some utilities for estimating and working with FRFs.
See in particular
GitHub - baggepinnen/ControlSystemIdentification.jl: System Identification toolbox for LTI systems, compatible with ControlSystems.jl
and some example notebooks here
ControlExamples.jl/ at master · JuliaControl/ControlExamples.jl · GitHub

The implementation of the data structure for a nonparametric FRF and methods for estimating and manipulating them are contained in this file

4 Likes

Some more less standard terminology
Spectrum - FFT of time waveform (Sx)
Auto Power - Sx * conj(Sx) (Gxx))
Cross Power - Sx * conj(Sf) (Gxf)
subscript f is force
subscript x is response
H1 - FRF where all the noise is assumed on the excitation side
H2 - FRF where all the noise is assumed on the response side

For an intuitive understanding of what is going on, I still really like the HP technical note series AN 243, HP Application Notes. If you look in the Appendix of the Fundamentals, it will give the equation for H1. Another really good place to browse is Technical Papers | Vibrant Technology. Mark Richardson is a former HP person. FFTW, available in Julia is known as the best FFT algorithm around. It is also in many other packages like MATLAB and Octave. If you need more help on it than Julia gives, you may get some good hints from looking up MATLAB’s help.
The paper by Ron Potter on windows is a classic. http://papers.vibetech.com/Paper64-CompilationofTimeWindowsandLineShapesRPotter.pdf. For the exponential and force window the paper “CAUSE AND EFFECT OF APPLYING THE EXPONENTIAL WINDOW
TO AN IMPACT FORCE SIGNAL” by Fladung is excellent. There is also some good theory on using a combination of coherence and number of averages for the statistical validity of the data.
For more in depth information, the course offered by UC SDRL on signal processing and experimental modal analysis is excellent. I have not seen it, but by reputation Peter Avitable’s book is excellent. There is also the book “Introduction to Operational Modal Analysis” by Brincker and Ventura. You should be able to find many equivalent Matlab (and Python) algorithms on Matlab Central, Anders Brandt’s toolbox “AbraVibe” ABRAVIBE Toolbox and Tom Irving’s website, http://www.vibrationdata.com/.

2 Likes

Thanks! I’ll have a closer look this weekend.

Great links! Thanks @Jake.
I received the Modal Analysis Basics book by Avitabile, and I have skimmed through some of it. Very good reference for me, but it made me realize I should be looking into the stochastic inputs in operational modal analysis, because most of my testing was a roving hammer without force input monitoring. I did put an accelerometer on the experimental hammer, and I have that data, but I would be curious to see if I can derive the natural frequency of the structure from the non-instrumented hammer.

Another curiosity is that I have been keeping two notebooks - one in Mathematica, and the other in Julia/Pluto, and the exact same steps are giving different results out of the FFT portion. I suspect Julia is doing a 1D FFT and Mathematica is doing a 2D FFT? Am I correct in assuming that I would do:

[Abs[FFT[list]]]^2 in either language for the power spectrum?

I am looking through both documentations for the answer.

I am using a list of data [{time, Az},{time, Az}…] where Az is the g reading in the z-axis. I could make the time element a period (1/256 Hz), or could I simply index it as 1, 2, 3, etc.? Am I wrong in assuming that a 1D FFT only cares about the list [{Az},{Az},…], and I can eliminate {time} understanding that the x-axis should be number of samples or the period (1/256 Hz) or ~0.004s? Thanks!

I cannot speak for Mathematica as I have not used it. In Matlab or Julia your steps would be something like this:

Gyy = zeros(Float64, floor(Int, length(bs)/2))
for i = 1:numberofaverages
bs is the blocksize or number of points per block
data = timedata[1:bs] .- mean(timedata[1:bs)
data = data .* timewindow  # timewindow is the window used.  if your signal is transient you can leave it out for now.  For random data the hanning window is often used, for your data the exponential window is the required window
newfft = fft(data,bs) ./ bs   # using FFTW and Matlab help
newfft = windowcorrectionfactor .* newfft # windowcorrectionfactor = 1/mean(timewindow); in some cases this is already built into the window
newfft = 2 .*[1:floor(Int, length(newfft)/2)]  # use only the positive frequency lines and muliply by 2 to compensate
Gyy = Gyy .+ real(newfft .* conj(newfft))  # works
Gyy = Gyy .+ abs2(newfft)  # probably better 
end
Gyy = Gyy ./ numberofaverages

This code should be at least roughly Julian though I did not check for errors.

The data start is found by a trigger location and then backup a few samples.

Note that when calculating the FRF’s the window correction factor cancels so it is not a concern for this calculation. Technically the multiply by 2 factor should not be done on the zero frequency bin, but we have set it to zero by subtracting the mean so that is not a concern. If your impact is a wholly contained transient, ie has decayed by the end of the block then a uniform window (no window) is in principle valid. An exponential window is still preferred as it will get rid of noise, but that can come later.

The time axis is not important for the FFT. It only cares about the Vector going in. You will need to recreate the frequency axis. df = 1/T where T is the duration of the block of time data. In your case where fs = 256 Hz and a blocksize of 1024 time samples = 4 s. Then df = 1/4 = 0.25 Hz. This is the frequency bin spacing.
With your data you can calibrate your modal hammer using F=ma to convert the measured acceleration to a force input. Or just assume that the hammer accelerometer output is force and then all the calculations will be correct, the results will be consistent but have a scaling factor deviation. No big deal for your learning.
You may want to read a bit more of Peter’s book and perform a graphical modal analysis on your data. I would probably do a traditional modal analysis, but have seen an OMA analysis applied to test results like this.
One you have the auto power spectrum for one point, or just the spectrum you should be able to see what the resonant frequency is.

1 Like

I am working my way through Peter’s book. It’s a lot for a beginner to digest, but I am making progress. Unfortunately, most of my data is just the response output from striking the structure at measured intervals with a non-instrumented hammer. Fortunately, I did get some useable signals with the makeshift hammer to which I attached an accelerometer. My primary goal is to determine the structure’s natural frequency in the horizontal and vertical axes. I have looked at OMA techniques, but they require sufficient ambient excitation to work with, but I basically had a static structure, not in use, simply struck with my non-instrumented hammer.

I now need to process the accompanying response spectra from the structure. I am still at a mental impasse on how to actually implement the conversion of my data/graphs of FFT vs. Sample No. or Accel vs. Time into the frequency domain. I am using Pluto to work through this, and it has been great for this exercise.

It does look like you have two nice resonances.

What you are looking at here is the positive part of the FFT and then starting at the right hand side and going left the negative part of the FFT. So sample 1 is at frequency 0
sample 2 is df (delta frequency)
sample 3 is 2df

sample n/2+1 is Fnyquist
sample n/2+2 is Fnyquist - df

sample n-1 is -2
df
sample n is -df
ie there is no -0 frequency.

So you would normally not display the negative frequencies, and often it is easier to delete them but to keep the same energy multiply the positive frequencies by 2.

To find df go to your time block and find how many seconds long it is. Let us call this T seconds. Then df is 1/T. In other words, if your time block above is 4 seconds long then the frequency lines are 0.25 Hz apart giving the frequency points in your plot as 0:0.25:N where N is the number of lines after you have truncated the array to be all positive values.

A typical setup that I might use would be a 1024 point FFT with a sampling rate of 256 samples per second which gives the nice 4 second block. This gives a Nyquist frequency of 128 Hz, and most spectrum analyzers would consider this as an Fmax of 100 Hz, where the lines between 100 Hz and 128 Hz are affected by the anti-alias filter. Hopefully Peter covers these things in his book and my explanation will be helpful.

As a simple exercise, assume all your averaged hammer impacts at each location was at a constant impulse (force). Pick the values from the spectra at each resonant frequency. Then sketch your structure. Then put little arrows at each location where you collected your data. Make the length of the arrow proportional to the value you picked off the FFT. The arrow height will be a simple hand sketch depiction of the mode shape. Obviously if you have the impact force as well, then plotting H1 is more accurate, but start with the data you have.

2 Likes

Thanks so much @Jake! I’ll look this over more carefully when I get home.
I will share my Pluto notebook to clarify my study. Peter’s book goes over the modal shapes from data, which I find amazing. My primary goal is to determine the horizontal and vertical natural frequency of the structure. I am hoping it is well above 5 Hz for both. The plot in the previous post is of the excitation input I captured from my experimental hammer. I was happy to see some shape similaritites to graphs of input forces in Peter’s book. Right now I am doing side-by-side plots in Pluto of the input force and the response waveforms or spectra. I am curious to see if I can derive the natural frequency solely form the response waveforms without the input force data. If not, I will try and do the next simple step - without windows, and see if there is an obvious fundamental or natural frequency in the response spectra. Thanks for the explanation of deriving the frequencies! I’ll give it a try later.

This graph does not seems to make much sense.
Is there a symmetry wrt ~270 Hz?!

I was assuming it was because I took the absolute value of the FFT, and this caused it to mirror itself?

Anyway, here’s where I am so far in trying to do a modal analysis in Julia using a Pluto notebook. I have input data for 1 second @ 100 Hz (100 samples), and 1 second @ 256 Hz (256 samples) of the steel structure’s response at one point. I am guessing 1 second is too short a window of data, and I still need to convert the time vs. input and time versus response data to the frequency domain. Is that at the absolute value of the FFT stage I did, or before? I am trying to figure out how to do the H1 FFR, but I am not sure how to multiply and divide complex conjugates and the like for Gxx and Gxy. More reading and experimenting to do. Does anyone know of a program (free hopefully) I can use to do the modal analysis to check my work? Thanks everyone!

ModalAnalysis-1|666x500

@Jake You had summarized the steps as follows:

You had also showed me how to calculate an H1 FRF here:

Gyy = zeros(Float64, floor(Int, length(bs)/2))
for i = 1:numberofaverages
bs is the blocksize or number of points per block
data = timedata[1:bs] .- mean(timedata[1:bs)
data = data .* timewindow  # timewindow is the window used.  if your signal is transient you can leave it out for now.  For random data the hanning window is often used, for your data the exponential window is the required window
newfft = fft(data,bs) ./ bs   # using FFTW and Matlab help
newfft = windowcorrectionfactor .* newfft # windowcorrectionfactor = 1/mean(timewindow); in some cases this is already built into the window
newfft = 2 .*[1:floor(Int, length(newfft)/2)]  # use only the positive frequency lines and muliply by 2 to compensate
Gyy = Gyy .+ real(newfft .* conj(newfft))  # works
Gyy = Gyy .+ abs2(newfft)  # probably better 
end
Gyy = Gyy ./ numberofaverages

…which I almost have working.

I have two gaps in my understanding:

  1. 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.

  2. For the input-output (force/response) data sets, I have followed your lead as quoted above in defining Gxx, I have applied the same to Gyy. The Cross power spectrum is then defined as Sy * S(conjugate of X) in Peter’s book. This leads to H1 as Gyx/Gxx, or the Cross power spectrum divided by Gxx. For H2 FRF there is a Gxy. How does Gxy differ from Gyx? Could this be a typo?

Here is the code that works derived from your Gxx example:

begin
	Gxx = zeros(Float64, floor(Int, 100))
	bsx = 100
	datax = HInput[:,2] .- mean(HInput[:,2])
	newfftx = fft(datax)
	#newfft = windowcorrectionfactor .* newfft
	#newfft = 2 .*[1:floor(Int, length(newfft)/2)]
	#Gyy = Gyy .+ abs2(newfft)
	Gxx = Gxx .+ real(newfftx .* conj(newfftx))
	plot(Gxx, title="Input Power Spectra",x_lims=(0,50),label="Gxx")
end

The abs2 line bombed out, so I went with the conj method.
I didn’t know how to define “windowcorrectionfactor” for the previous line.
There was a type error with the floor function line.

All-in-all it seems to be sufficient for the example I am using. A mistake I made was not setting my input accelerometer to its maximum 200 Hz frequency (it defaulted to 100 Hz) given my output tri-axial accelerometer was set to 256 Hz. They are both durations of 1 second, except now to do the H1 FRF I need Gyx which means I am going to be doing a *. and conj on a 256 vs. 100 element matrix (vector). How can I rectify this? Normalization or some sort of upsampling or downsampling?

This has been a great Julia/Pluto learning exercise, and since most of my questions have turned to FRFs and such, I will wait until I have more Julia questions. Thanks for your help!

Perhaps I did not explain myself very well, there absolutely is symmetry and it is expected. To the left of 270th frequency bin it is the positive frequencies, and to the right the negative frequencies going from the highest frequency to the lowest, non-zero frequency. The amplitude for the positive and negative frequency is the same, thus the symmetry.

1 Like

@Jake I finally made it through this self-teaching exercise, thanks to you and everyone who helped. To give back a bit to those who may stumble upon this in need of help, here are some clarifications:

The top line worked, and the only correction to the second is to add a ‘.’ after ‘abs2’ -

Gyy = Gyy .+ abs2.(newfft)

I had to resample my experimental instrumented hammer, since both accelerometers defaulted to their default setting so 256 Hz for the tri-axial accelerometer on the structure, and 100 Hz for the one I put on the back of my soft-face hammer. I had them both at 200 Hz the day before! The resample() function did a good job of keeping the data close to the original after the resample.

Is there a way to share Pluto notebooks in posts?
Here are some screen shots of the end of my notebook showing the linear Fourier spectra, the cross power spectra, and the FRF (H1). I learned so much from @Jake and Peter Avitabile’s book Modal Testing: A Practitioner’s Guide.

image

image

So it looks like the the first highest peak in the FRF at H1[13]. Given a sample rate of 256 Hz and a sample period of 1 second, does that mean that each unit on the x-axis is 1 Hz, so H1[13] is 13 Hz as my natural frequency?

1 Like

Cool! Would you mind sharing the data for others to play with?

Not having followed the thread cannot comment on the coincidence of having 256 samples and also 256 Hz as the sampling frequency.
In any case the general formula below applies to any FFT dataset:

N = 256         # number of samples
f0 = 256        # sampling rate [Hz]
Δt = 1/f0       # sampling period [s]
Δf = f0/N       # Δf = 1/(N*Δt), frequency step in [Hz] between FFT samples

NB: the above are physical frequencies, not angular frequencies