Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

It is not always obvious there is a peak frequency, what to do about those trials? #40

Open
moorepants opened this issue Jan 12, 2025 · 12 comments
Assignees

Comments

@moorepants
Copy link
Contributor

For example, session 006 klinkers:

image

There are basically two dominate frequencies with similar amplitude peaks.

Session 007 klinkers is similar:

image

Session 007 tarmac:

image

@moorepants
Copy link
Contributor Author

moorepants commented Jan 12, 2025

One thing to try: smooth the spectrum more.

Edit: this has been done.

@moorepants
Copy link
Contributor Author

We could calculate the median (basically the 50% area mark). In fact, I think that is what I've often reported in the past for such an analysis.

@moorepants
Copy link
Contributor Author

This function actually calculates a box plot type values of a given curve: https://dynamicisttoolkit.readthedocs.io/en/latest/dtk.html#dtk.process.curve_area_stats

@moorepants
Copy link
Contributor Author

moorepants commented Jan 14, 2025

Here is an example of calculating a box plot from the area under the frequency spectrum:

>>> freq, amp, _, _ = tr.calculate_frequency_spectrum("SeatBotacc_ver", 400, smooth=True)
>>> from dtk.process import curve_area_stats
>>> curve_area_stats(freq, np.atleast_2d(amp).T)
{'2p': array([3.44238281]),
 'lq': array([8.66699219]),
 'median': array([27.19726562]),
 'uq': array([83.203125]),
 '98p': array([182.39746094])}
>>> fig, ax = plt.subplots()
>>> ax.bxp([{'med': 27.197, 'q1': 8.67, 'q3': 83.20, 'whislo': 3.44, 'whishi': 182.40}], vert=False, showfliers=False)
{'whiskers': [<matplotlib.lines.Line2D at 0x7dbb0fe35a90>,
  <matplotlib.lines.Line2D at 0x7dbb0fe35bd0>],
 'caps': [<matplotlib.lines.Line2D at 0x7dbb0fe35d10>,
  <matplotlib.lines.Line2D at 0x7dbb0fe35e50>],
 'boxes': [<matplotlib.lines.Line2D at 0x7dbb0fe35810>],
 'medians': [<matplotlib.lines.Line2D at 0x7dbb0fe35f90>],
 'fliers': [],
 'means': []}
>>> plt.show()

image
image

@moorepants
Copy link
Contributor Author

Riender suggested finding multiple peaks for some trials (this could help https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html).

@moorepants
Copy link
Contributor Author

From Riender:

In most plots the peak frequency is credible.
However some cases show 2 peaks of similar magnitude (see below).
Actually the second frequency seems to be twice the first frequency, and could present is first harmonic frequency.
Here it will make sense to report at least 2 peak frequencies.
In the statistics this will be complicated but we will find a trick.
This happens only for a few cases – and we can manually select the most relevant frequency.
I guess this will be the lowest of the 2 frequencies.

In the old rusty the peak is often not clear at all, especially on klinkers.
It may be a waste of time to analyse peak frequency for old rusty further.
Overplotting power as function of frequency for old rusty with new strollers will be much better.

@moorepants
Copy link
Contributor Author

Riender just came in an said that for the trials that have 2 peaks, we should pick the first one even if the magnitude is a bit lower than the second peak. For example:

Image

Those two should be updated to find the first peak. He said that would satisfy his concern.

@moorepants
Copy link
Contributor Author

I've switched to ISO weighted spectrums so this changes the frequency spectrum plots and thus the peak frequency and bandwidth. Need to review new plots.

@ForensicBMechE
Copy link

  1. What's the reasoning behind taking the first peak? And in that case, if you don't take the highest, how do you define a peak? This must be very clearly defined.
  2. Amplitude units are confusing. If we normalize them, shouldn't they have no units?
  3. What "threshold frequency" does the green line indicate?
  4. There are two obvious peaks standing out, but I would even say there's a third, much smaller one. And if I measure where they are from zero in pixels, I get about 25, 50 and 75... so aren't we just looking at a fundamental frequency and to harmonics?

@moorepants
Copy link
Contributor Author

What's the reasoning behind taking the first peak? And in that case, if you don't take the highest, how do you define a peak? This must be very clearly defined.

The reasoning is that the first peaks likely correspond to the same mechanical cause and that following peaks may be harmonics.

Amplitude units are confusing. If we normalize them, shouldn't they have no units?

Frequency spectrums have units, raw FFT will be (physical unit)/Hz. I've normalazid wrt to Hz so we get the unit of the signal.

This has a good explanation if you want to read about it: https://docs.scipy.org/doc/scipy-1.15.0/tutorial/signal.html#continuous-time-sine-signal I present that "Amplitude Spectrum" shown on that page.

What "threshold frequency" does the green line indicate?

This is explained in the paper and my prior presentations. I choose 80% of the area under the smoothed spectrum. It is just a way to show where most (being 80%) of the power is in the frequency spectrum. I will change the name later.

There are two obvious peaks standing out, but I would even say there's a third, much smaller one. And if I measure where they are from zero in pixels, I get about 25, 50 and 75... so aren't we just looking at a fundamental frequency and to harmonics?

For the January 31st deadline, I don't intend to dig in fine detail about the peak analysis.

@moorepants moorepants self-assigned this Jan 20, 2025
@moorepants
Copy link
Contributor Author

After applying the ISO filters this is less of an issue. Here are ones that may not be getting the peak we desire:

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

@moorepants moorepants removed the jan31 label Jan 31, 2025
@moorepants
Copy link
Contributor Author

For the Jan 31 report, I just said we sometimes pick the second peak. Still needs addressing for the journal submission.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants