2

I'm studying mechanical engineering and I would like to determine the magnitude and phase of a signal. Fourier analysis has been mentioned many times during my study and I think I know what it is about but I never had to apply it myself, which is (as usual) way more difficult than following the analysis someone else did;) I tried to find answers to my questions, but haven't been successful so far. I started to get more confused the more I read so I decided to ask somebody for help. I don't expect a final solution, like source code, but could please someone point me in the right direction?

Technical background:
I would like to analyze the bending moment of the rotor mast of a helicopter. The bending moment is calculated from the signal of a strain gauge, which rotates with the rotor. Apparently, the signal of the bending moment is overlain by the rotor revolution. What I need is the phase and the magnitude of the signal.

Approach:
I have time-discrete flight test data of the bending moment and a signal, which indicates, when rotor blade number one is located above the tailboom. I applied a fast fourier transformation to the data of one revolution and would like to determine phase and magnitude from the imaginary and real part of the fourier coefficients. The computations are done by a python script.

Questions:
As far as I understood, the first fourier coefficent equals the mean value of the signal and the next coefficients equals the coefficient of the first harmonic frequency, which would be the rotor frequency, since I analyse the data of one revolution. If I change the number of terms considered, the coefficents change. This makes perfectly sense to me, but affects phase and magnitude as well. So choosing the second coefficient returned by my fft-function can't be right. What am I missing?

Thanks in advance!
Johannes

How do I choose N? Does N equal the number of measured values?

Never heard of atan2, but I found information about it...

If the frequency of the coefficient equals k/N, I know, that the coefficient with k=1 stands for the rotor rpm. I did the calculation for 36 revolutions and in 7 cases, the coefficient with k=0 has had the larges magnitude by far. So I'll increase the sampling rate. Fortunately, the raw flight test data is sampled with 25kHz and the data I used unitl now is sampled with 1kHz.

Does it affect the quality of the result, if I calculate the phase and magnitude for every revolution seperately and calculate the mean afterwards instead of calculating the phase and magnitude for the whole flight at once? I decided to calculate the phase and magnitude for each revolution to be able to calculate the standard deviation to asses the quality.

thanks alot!

Here some final notes:

  • the time for the signal calculated from the fourier coefficients starts at t=0s no matter what time the actual signal starts. This is what leads to the phase shift in the diagramm ("Signals") I posted above.

  • calculating the mean of all phases for each revolution leads to a tremendous standard deviation (up to 180°). This is caused by the discontinuity of arctan. It can be solved by calculating the mean of the real and imaginary part of the fourier coefficients befor calculting the phase. I'm not yet sure under which circumstances calculating the mean of the fourier coefficients is valid...

  • arctan2 from python expects the arguments in a different order than posted above or as other programms like excel... I used numpy.angle(imaginary number) and numpy.abs(imaginary number) in the end anyway...

Thanks alot for your help! Johannes

Cameron Buie
  • 102,994
Johannes
  • 21
  • Couple of questions: 1. What do you mean when you say that the "signal of the bending moment is overlain by the rotor revolution"? 2. Are you assuming that the signal of the bending moment is sinusoidal? If not, how would you define the "beginning" of the signal? I don't think you can solve this problem without defining those things. – Adrian Keister Jul 19 '13 at 17:17
  • yes, the signal is sinusoidal. In forward flight, the mast is bended about the lateral axis of the helicopter. Since the strain gauge is mounted to the rotating axis, the signal isn't equal to the magnitude of the bending moment, but a sinusoidal signal oscillating with the rpm of the rotor. I don't know, what you mean with "beginning". I know the beginning of a revolution due to a second signal.

    I hope this answers your questions. It's not that easy for me to describe, since English isn't my mother tongue;)

     – Johannes Jul 19 '13 at 23:23

1 Answers1

2

Here's the algorithm I would use:

  1. Take the DFT/FFT of your sinusoidal signal. If your sinusoidal signal is $x_{n}$, then your transformed signal is $$X_{k}=\sum_{n=0}^{N-1}\left[x_{n} \cdot e^{-2\pi ikn/N}\right].$$
  2. You want the basic amplitude and phase of the strongest frequency component of your signal. So, find $\displaystyle \max_{k} |X_{k}|$, call it $X_{\max}$.
  3. The amplitude of the strongest component of your signal is given by $$\frac{|X_{\max}|}{N}=\frac{\sqrt{\text{Re}^{2}(X_{\max})+\text{Im}^{2}(X_{\max})}}{N},$$ and the phase is given by $$\arg(X_{\max})=\text{atan2}(\text{Im}(X_{\max}),\text{Re}(X_{\max})),$$ where you may need to offset this value by the signal that tells you when the rotor blade one is located above the tailboom. I don't know if you're familiar with the $\text{atan2}$ function, but it's basically the arctangent function without the $\pi$ ambiguity that it normally has. I'm pretty sure Python has this function natively.

If you happen to be interested in the frequency of the sinusoid represented by $X_{\max}$, it is $k/N$ cycles per sample, where $k$ is the index of this component.

Incidentally, if you find that the first component of your transformed signal is by far the strongest, that's probably an indication that you're not sampling fast enough.

Adrian Keister
  • 10,099
  • 13
  • 30
  • 43
  • 1
    How do I choose N? Does N equal the number of measured values?

    Never heard of atan2, but I found information about it...

    If the frequency of the coefficient equals k/N, I know, that the coefficient with k=1 stands for the rotor rpm. I did the calculation for 36 revolutions and in 7 cases, the coefficient with k=0 has had the larges magnitude by far. So I'll increase the sampling rate. Fortunately, the raw flight test data is sampled with 25kHz and the data I used unitl now is sampled with 1kHz.

     – Johannes Jul 20 '13 at 16:04
  • 1
    Does it affect the quality of the result, if I calculate the phase and magnitude for every revolution seperately and calculate the mean afterwards instead of calculating the phase and magnitude for the whole flight at once? I decided to calculate the phase and magnitude for each revolution to be able to calculate the standard deviation to asses the quality.

    thanks alot!

     – Johannes Jul 20 '13 at 16:05
  • $N$ is the number of samples. If you do an FFT, then you may have to zero-pad in order to get $N$ to be a power of $2$. How fast are the rotor blades rotating? Generally, you need to sample at least ten times as fast as the frequency of the phenomenon in which you're interested. I know Nyquist says twice, but that's a theoretical frequency, without taking noise into account. Can you post a sample FFT plotted against $k$? My step # 2 might not actually be correct. If you have a well-defined peak somewhere out there, that might be what you seek. There are ways to find those. – Adrian Keister Jul 20 '13 at 18:25
  • As for your second comment, I would probably analyze each rotation separately, like you're doing, for the reasons you gave. – Adrian Keister Jul 20 '13 at 18:28
  • 1
    The rpm of the rotor is adapted to e.g. speed and air density. But I have the signal of the position of one blade, so I can calculate the rpm for each revolution. The rpm varies between 5.4 Hz and 6.77Hz. So the minimum sample rate equals 67.7Hz. The sample rate of the data I analyse equals 1kHz. But there are still 7 out of 36 revolutions with a way larger first component than the others... – Johannes Jul 20 '13 at 19:49
  • 1
    What do you mean with "sample FFT plotted against k"? The coefficients in the complex plane?

    Until now, I selected the data of each revolution, without considering, that N should be a power of 2. Should I crop the interval or take values of the next revolution into account?

     – Johannes Jul 20 '13 at 19:56
  • 1
    What do you mean with "well defined peak"? The sinusoidal character of the values of the beinding moment? Unfortunately, it is not always distinct. How ever, this hardly happens, so I could simply exclude those revolutions. Thanks for your support! What you wrote helped me way more, than pretty much everything I read in the last couple of days in various books and web pages... – Johannes Jul 20 '13 at 20:03
  • Ok, sounds like your sampling rate is adequate. I would still, if your computing resources can handle it, use all the data in your analysis. By the plot, I mean the absolute value of each coefficient plotted against its index $k$. As for $N$ being a power of two, your FFT might handle that already. Check the documentation. You can compute a DFT without using an FFT algorithm. If your sample size is small, you might want to do that anyway. By "well-defined peak", I mean a peak in the plot I just mentioned: $|X_{k}|$ versus $k$, that is not at the far left-hand end of the plot. Glad I can help! – Adrian Keister Jul 20 '13 at 20:24
  • Is this the type of plot, you asked for? frequency spectrum

    here a detail of the low frequencies frequency spectrum detail low frequency

    I export the data with the sampling rate of 1kHz with a software tool to csv. The files have 3GBs already and analyzing the data takes about 20 minutes. I can crosscheck the magnitude and phase with another tool. If the deviation is significant, I'll increase the sampling rate.

     – Johannes Jul 20 '13 at 21:49
  • The python function does modify the data, so N equals a power of two. I don't know how, but the documentation says it does;)
    Why would I do a DFT without using a FFT, if there is a function, I can import. Do you mean I write the algorithm myselfe?     – Johannes Jul 20 '13 at 21:57
  • I almost forgot: The data I plotted is from a interval of 6 seconds of lengths. I wasn't sure, if you need information about the characteristics of the signal in general or of one revolution... – Johannes Jul 20 '13 at 22:06
  • Yes, those are the kinds of plots I was asking for. I'm guessing the second plot is a zoom in on the first? Looks like the first major peak is around 6Hz, and I think that would be the one in which you're interested, based on your rotor speed. It actually is the maximum peak, and it's separated enough from the first bin to be a distinct peak. More data points would help out, for sure. There's no particular reason to do a bare DFT minus the FFT algorithm, most of the time: zero padding does not hurt your data. – Adrian Keister Jul 22 '13 at 14:18
  • Almost there;) I still got some phase displacement:(

    Frequency spectrum,

    Raw data,

    low frequency spectrum,

    signals.

    Yes, the second plot is a zoom in on the first. I figured out, that I can use a rea FFT instead of a complex FFT. The complex FFT calculated those peaks on the left hand side, which are complex conjugated coefficienc only.

     – Johannes Jul 23 '13 at 09:43
  • And I normalized the coefficients by dividing by the number of data points considered and multiplied by two, so the magnitude of the fourier signal equals the magnitude of the fourier signal. – Johannes Jul 23 '13 at 09:51