Phase Shift Problem

[ptc-5189599]

Hi,

I'm currently working on a project which involved some oscillation testing of a scale model. The issue I'm having is that there is a phase lag in the output data, i.e. the maximum force does not correspond with a zero crossing of amplitude as it should. Rather, it is either in phase with the force, or anti-phase with the force. I was hoping there was some way to sort this in MathCAD, it isn't a straightforward phase shift as the phase lag differs with each test run. Is there some coding that would allow to me to check for the various times at which amplitude = 0, and then adjust the data in the force column so that Fmax correlates with A = 0? The frequency and thus period are the same, so by sorting the phase shift corresponding to 1 point, the rest of the data should fall in place accordingly. Hopefully, this could be applied regardless of which data file is used so that for all test runs, Fmax correlates with a zero crossing of amplitude, and the results can be readily obtained.

I've created a very basic mathcad sheet to illustrate the problem. The input values for 'Series' and 'File' can be obtained by checking the arrays on page 3 i.e. the 'ID' column in each array refers to the file input, and R1/R4 refers to series 1/4 respectively.

The non-dimensionalised Force vs Amplitude shows a tilt in the graph which highlights the problem - this graph should be an ellipsis around the y axis passing through the origin.

Any help on this is greatly appreciated. I am using MathCAD 15

Thanks

[Werner_E]

Could you show in a sheet an example how you would manually "shift" your data to achieve what you want?

Especially how you would determine the point of maximal force and how the shifting should be applied. As I understand it, shifting Fs to the right would mean getting rid of some values at the end of Fs and the same amount at the beginning at As. Shifting to the left - the other way round. Obviously ts would have to be shortened accordingly, too. So you would end up with shorter vectors not showing full periods as it seems to be now.

BTW, is vector t sorted in every data file? Then I would guess that FindRange() could be written simpler. In any case this routine should not be defined using a global assignment! Its interesting that the results differ slightly (and the first index seems not to be correct in any case) depending which way the function is defined:

Addendum: I just realized that the values t1 and t2 you pick from the tables at the end of your sheet do NOT mark begin and end of full periods. That way your way of getting rid of the constant component by subtraction of mean() does not work as intended. Better seen in the graph if you chose "Crossed" instead of "Boxed".

Maybe better to use As:=As-(max(As)+min(As))/2

Same with Fs, but here its more problematic as its very noisy.

[Werner_E]

Here is a shift routine. Not sure if its OK to apply it to Fs/As/ts rather than to A/F/t.

Some questions:

1) Whats the meaning of the values t1 and t2 you pick from the table? Is it mandatory that t2 has to be the last value in ts or may the vector be a but shorter (as is the case when using my shift routine)

2) Whats the output (to be saved) of the overall calculation? (Fs/As/ts or something else? Should it be stored)

3) Is the frequency of the signals 2*pi/5 in every case?

4) Is it really OK to change the measured data by simply shifting force values along the time axis?

Find also attached an animation for your pleasure - as you can see it wouldn't be that easy to decide which shift should be considered best. A perfect circle is not to be expected anyway.

[RichardJ]

Here's another approach that will find the shift with a very high degree of accuracy. Depending on the goal of the analysis, it may also provide parameters that would lead directly to the required solution, also with a high degree of accuracy, and possibly with no need to actually shift the data arrays.

[Werner_E]

Here is an appraoch using genfit to find the sinusoidal fit. I guess using FFT could do the job, too.

Additional an avi showing the results for all 27 data sheets.

[RichardJ]

Here is an appraoch using genfit to find the sinusoidal fit.

Personally, I never use genfit. I find minerr easier to set up for anything other than very simple fits, and much more versatile.

I guess using FFT could do the job, too.

Not with much accuracy, because the data does not have an integer number of periods, and the number of periods is small. However, I have found in the past that a least squares fit of a sinusoid can be very sensitive to the guesses (especially if the sinusoid is undersampled), in which case an FFT is a great way of getting starting guesses. Some years ago I posted a worksheet that uses this approach. See http://communities.ptc.com/message/134317#134317

[Werner_E]

I guess using FFT could do the job, too.

Not with much accuracy, because the data does not have an integer number of periods, and the number of periods is small.

Yes, think you are right.

Added: But as "A" is a very clear signal and knowing that "A" and "F" would always have the same frequency we could cut the vectors so we get data with integer multiples of the period.

[Werner_E]

I guess the call to the shift routine before the A/F-plot at the end should read

[vector(Fs,As,ts)]:=shift(Fs,As,ts, - round(PointShift))

I tried getting guess values using cfft, but to no avail. They are too far off. Genfit (and it applies to your solve block, too) is very sensitive especially about the guess value for frequency. I uses omega (2.pi.f) and this guess value has to be in the range 1.1 to 1.4 to make it work. cfft would suggest a value of 0.0094!

See the "FFT guess" section in the attached sheet - maybe I am doing something wrong.

I think I sticked to the sheet you mentioned, but I guess that the frequency is not just inaccurate but calculated wrong (in my sheet). Shouldn't it be multiplied by the (unknown) sample frequency? But you didn't in the sheet you referenced and it seems to work there without.

So the question is if there is a possibility to get the frequency/period of a signal if we do not know the sampling frequency and sampled data are not integer multiples of a period.

Data set 4 contains some heavily oscillating data (esp. 4/1059 and 4/1061) where the A/F-plot won't be suitable for visual judgement of the fit and shift.

[Werner_E]

OK, I should have know better. Of course we need to multiply by the sample frequency and of course we know it! What we don't know is the samples per period - think I got confused a bit. In Richard's reference there was no special "time" vector, just the indices. So sampling frequency was naturally 1.

It works pretty well now - even for noisy Fs's.

What I am still not sure about is the phase shift. I thought it would be the argument of the complex value with maximal magnitude in the cfft (without the first), but according to Richards reference its the complement of that argument to 90° (the arguments in the atan2() function he uses are switched). Anyway - either of both ways would provide nice, nearly perfect fits for one data set and awfully wrong ones (up to half a period away) for others. An not only for the vectors Fs but also for As, which is a very well behaved measured signal in all cases and so I would expect a better result.

Any explanation for that?

Edit: Seems that Richard is right with the phase angle and the complement to pi/2. At least we get perfect results if the data length is an integer multiple of the period. If the data length is not an integer multiple the phase angle can be way off. Is there any relationship we could build on? We get a really good guess for the frequency, so we could calculate how many periods we have. Any way to correct the phase angle with that knowledge?

[RichardJ]

What I am still not sure about is the phase shift. I thought it would be the argument of the complex value with maximal magnitude in the cfft (without the first), but according to Richards reference its the complement of that argument to 90° (the arguments in the atan2() function he uses are switched).

After the FT the phase of a cosine wave would be zero, but the fitting function is a sine wave.

Anyway - either of both ways would provide nice, nearly perfect fits for one data set and awfully wrong ones (up to half a period away) for others. An not only for the vectors Fs but also for As, which is a very well behaved measured signal in all cases and so I would expect a better result.

Any explanation for that?

Edit: Seems that Richard is right with the phase angle and the complement to pi/2. At least we get perfect results if the data length is an integer multiple of the period. If the data length is not an integer multiple the phase angle can be way off. Is there any relationship we could build on? We get a really good guess for the frequency, so we could calculate how many periods we have. Any way to correct the phase angle with that knowledge?

The problem is that the FFT does not know about the time axis, and the calculated phase is therefore referenced to the first point being at zero (so, in this case, not even the indices of the array, because ORIGIN=1). So there is an additional phase shift of 2*PI*(number of periods) before the data starts. I tried to correct for this, and kept getting unexpected results. I just noticed that you defined your function in an unusual way: sin(b2*(t+b3)), rather than sin(b2*t+b3), so the phase is the product of b2 and b3. No wonder everything I tried didn't work! Easier I think to switch to the function b1*cos(b2*t+b3)+b4. I have to step away from it for now though, but if you don't post a fix before I get to it I'll look at it later.

[Werner_E]

The problem is that the FFT does not know about the time axis, and the calculated phase is therefore referenced to the first point being at zero

I made a test sheet (attached, feel free to change sin to cos or whatever) and also found out that its not the integer multiple of the period which is responsible for the good phase shift result.

In contrary to the data of the initiator of this thread in that test sheet a data length different from an integer multiple of the period would give awful frequency results, no matter how many full periods would be included 😞

noticed that you defined your function in an unusual way: sin(b2*(t+b3)), rather than sin(b2*t+b3),

yes, I often prefer that representations as b3 that way would correspond w/o further calculation to the point of zero-crossing and is exactly the time dist sin(b2*t) would have to be shifted.

[RichardJ]

yes, I often prefer that representations as b3 that way would correspond w/o further calculation to the point of zero-crossing and is exactly the time dist sin(b2*t) would have to be shifted.

True. I guess it's six of one and half a dozen of the other. It's like origin 0 or origin 1: whatever you get used to is easier.

[RichardJ]

This fixes the guesses, but now there's something broken in the shift calcuation, and I haven't figured out what yet.

[Werner_E]

Great, now there isn't much to do for genfit anymore. apart from some slight frequency correction in 1/1076.

I think the "base" value is wrong now. It should not be doubled and absolute valued. Thats the reason I grabbed the value from the spectrum, not the Mag.

Also a (very very seldom) possible problem in the calculation of hi: If the constant component Mag0 or Mag1 in that sheet has by chance the same value as the absolute value of the max wave, you will get index 1 instead of the correct one. As already said - not supposed to happen that often, but that was the reason I spent variable Mag1.

I'll have a look at the now wrong shifting tomorrow - it might be because of the different function definition. But the switch from sin to cos should not change anything and tghe necessary replacing of cA3 by cA3/cA2 and the same for F did not do the job.

[Werner_E]

I'll have a look at the now wrong shifting tomorrow - it might be because of the different function definition. But the switch from sin to cos should not change anything and tghe necessary replacing of cA3 by cA3/cA2 and the same for F did not do the job.

OK, its tomorrow and apart from the change described above it was also necessary for a change in sign in shift_i because you defined cos(b2 - b3).

Think it should work now - of course only if the data length is an integer multiple of the period (as can be easily demonstrated with the cfft-test sheet).

Fortunately Martins data are all nearly perfect integer multiples of the period, otherwise the FFT results would not be of much help as of the sensitivity with respect to the guess value for the frequency. Probably it would be a better approach to use a brute force wade through the As values and look for two consecutives zero-crossings.

[RichardJ]

Cool

Fortunately Martins data are all nearly perfect integer multiples of the period, otherwise the FFT results would not be of much help as of the sensitivity with respect to the guess value for the frequency.

If the number of periods of an otherwise perfect sinusoid is non integer, after the FFT the dominant frequency will still be the underlying frequency of the sinusoid (>90% sure of that, but not 100%). In this case I noticed that if you look at the spectrum the peak is always remarkably symmetric. I guess the sampling was designed based on the exitation frequency. If that's not the case then the maximum point may be a poor estimate of the frquency, but interpolation can help to get a better frequency estimate. A simple approach to that is to zero pad the data before the FFT. That can lead to big oscillations either side of the peak (the FFT of a boxcar function being a sinc function), but those can be suppressed by windowing the data. I've done a lot of this type of analysis, so trust me that this approach works (and also never trust anyone that says "trust me" ).

A little tidying up to do, and this one's a keeper. There are a lot of experiments where one applies a sinusoidal perturbation and monitors the response (hopefully also sinusoidal, although not necessarily so). This worksheet contains a lot of the basic elements required to analyze such data (more so than the MTF worksheet I pointed to, which does not have a time scale). So it's destined for my library

[Werner_E]

Yes, its a handy utility.

But the data length being an integer multiple of the period would be mandatory, otherwise none of the four calculated values would be of any value - as for example in the attached sheet.

[ptc-5189599]

Hi guys,

Sorry about my lack of correspondance, I've recently moved house and struggling to get internet when out of office. I'll have a look through all the files just now and reply to your messages.

Thanks very much for your input, it's all greatly appreciated.

Martin

[RichardJ]

Fortunately Martins data are all nearly perfect integer multiples of the period, otherwise the FFT results would not be of much help as of the sensitivity with respect to the guess value for the frequency. Probably it would be a better approach to use a brute force wade through the As values and look for two consecutives zero-crossings.

No problem. Just needs a little modification. This is based on your earlier sheet, which contains an error. You fixed that in the later sheet. I fixed it a different way in this sheet.

[Werner_E · Accepted solution]

Thats great! Zero padding really works like a charm.

I implemented it in Martin's sheet to get better guess values even if his other data will be not that good behaved as the ones he posted.

[Fred_Kohlhepp]

A phase shift between force and amplitude may indicate damping, not faulty data.