Signal Fitting with sine curve (and how to find out the phase shift?)

Hi!

I have a temperature signal (Y-axis) consisting of 764 values that looks like this -

And the time (X-axis) of this signal is a bit nonperiodic. Each of the observation points are like this -

....

I have timetable and the temperature data (attached). Can anyone please suggest me how can I fit this signal into a sine curve and get the amplitude, mean, and phase shift? Any feedback will be greatly appreciated!

[Jan]

Fitting a sine wave to a non-periodic function is an expensive random number generator.

Did you search in the forum already? You'd find e.g. https://www.mathworks.com/matlabcentral/answers/121579-curve-fitting-to-a-sinusoidal-function#answer_128539

[Ashfaq Ahmed]

What is an expensive random number generator?

[Ashfaq Ahmed]

Yes, I actually used that code already. Excellent code. But I am not sure about how is the code cslculating the phase shift.

[Jan]

"Expensive random number generator" means, that the output has not more mathematical significance than rand(1, 4).

The function, Star Strider is fitting is:

fit = @(b,x)  b(1).*(sin(2*pi*x.*b(2) + 2*pi*b(3))) + b(4)

Then 2*pi*b(3) is the phaseshift.

[John D'Errico]

That signal does not have an even remotely fixed frequency. So there is no constant phase shift that could be estimated. Trying to do so would be meaningless.

Answer 1

[Star Strider]

I initially wnated to see if I could get this to fit using fminsearch, however I wasn’t guessing the correct initial parameter estimates and it would not return an acceptable result.

So I went to ga, since it can usually solve these problems when I cannot, and it succeeded brilliantly (in my opinion, anyway).

Here are those results —

LD = load(websave('SignalFitting','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1320795/SignalFitting.mat'));
Per = LD.Period;
Temp = LD.Average_Temp;

DP = diff(Per);
[lo,hi] = bounds(DP);
[h,m,s] = hms(mean(DP));

Stats = [lo hi mean(DP) median(DP) mode(DP)]

Stats = 1×5 duration array
     24:00:00   3456:00:00    441:16:44    384:00:00    192:00:00

TT1 = timetable(Per,Temp);
LenTT1 = size(TT1,1);
Ts = hours(h);
TT1r = retime(TT1,'regular','linear', 'TimeStep',hours(h))          % Resample Data To A Uniform Sampling Frequency

TT1r = 765×1 timetable
       Per          Temp   
    __________    _________

    1984-05-02       12.879
    1984-05-20       11.992
    1984-06-07       14.559
    1984-06-26       21.892
    1984-07-14       21.111
    1984-08-01        20.33
    1984-08-20       19.549
    1984-09-07       18.768
    1984-09-26       17.175
    1984-10-14       14.069
    1984-11-01       14.388
    1984-11-20        9.125
    1984-12-08        5.261
    1984-12-26       2.6238
    1985-01-14    -0.013477
    1985-02-01      -1.5444

DOY = day(TT1r.Per,'dayofyear');
LY = eomday(year(TT1.Per),2);
LYv = LY(diff(DOY)<1)-28;

YRv = [0; cumsum(diff(DOY)<1)];
max(YRv)

ans = 38

CDOY = DOY+365.25*YRv;

[dlo,dhi] = bounds(CDOY)

dlo = 123

dhi = 1.4161e+04

dCOD = dhi - dlo

dCOD = 1.4038e+04

[dlo,dhi] = bounds(TT1r.Per)

dlo = datetime
   1984-05-02

dhi = datetime
   2022-10-08

DPer = (dhi - dlo)/24

DPer = duration
   14038:30:00

figure
plot(CDOY, TT1r.Temp)
grid
xlabel('Cumulative Days')
ylabel('Temp')

L = size(TT1r,1);
Fs = 1/h;                                                           % Sampling Frequency (cycles/hour)
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FT_Temp = fft((TT1r.Temp-mean(TT1r.Temp)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

[pks,locs] = findpeaks(abs(FT_Temp(Iv))*2, 'MinPeakProminence',2.5)

pks = 4.7925

locs = 53

pkfreq = Fv(locs)                                                   % cycles/hour

pkfreq = 1.1515e-04

pkper = 1/pkfreq                                                    % hours/cycle

pkper = 8.6843e+03

pkper_d = pkper/24                                                  % days/cycle

pkper_d = 361.8462

figure
plot(Fv, abs(FT_Temp(Iv))*2)
grid
xlabel('Frequency (cycles/hour)')
ylabel('Magnitude')

% Elapsed Time:   4.923885300000001E+01		00:00:49.238
% 
% Fitness value at convergence = 102.6832
% Generations 				 = 163
% 
% 	Rate Constants:
% 		Theta( 1) =  9.42146
% 		Theta( 2) =  5.05100
% 		Theta( 3) =  1.30257
% 		Theta( 4) = 14.16108
% 
% 
% Elapsed Time:   3.035811500000000E+01		00:00:30.358
% 
% Fitness value at convergence = 102.6837
% Generations 				 = 132
% 
% 	Rate Constants:
% 		Theta( 1) =  9.42147
% 		Theta( 2) =  1.05100
% 		Theta( 3) =  3.30268
% 		Theta( 4) = 14.16108
% 
% Elapsed Time:   5.428350700000000E+01		00:00:54.283
% 
% Fitness value at convergence = 102.6837
% Generations 				 = 177
% 
% 	Rate Constants:
% 		Theta( 1) =  9.42146
% 		Theta( 2) =  6.94900
% 		Theta( 3) = 11.19732
% 		Theta( 4) = 14.16109

objfcn = @(b,t) b(1) .* sin(2*pi*t*b(2) + 2*pi*b(3)) + b(4);
% B0 = [max(TT1r.Temp)-min(TT1r.Temp); pkfreq/(365); 1; mean(TT1r.Temp)]
% B0 = [9.4; 1; 3; 14];
% [B,fval] = fminsearch(@(b) norm(TT1r.Temp - objfcn(b,CDOY)), B0)

B = [9.42147; 1.051; 3.303; 14.161];
TempFit = objfcn(B,CDOY);

figure
plot(TT1r.Per, TT1r.Temp, ':', 'DisplayName','Data')
hold on
plot(TT1r.Per, TempFit, '-r', 'DisplayName','Regression')
hold off
grid
xlabel('Per')
ylabel('Temp')
xlim([datetime(1985,03,01) datetime(1988,03,01)])

figure
plot(TT1r.Per, TT1r.Temp, ':', 'DisplayName','Data')
hold on
plot(TT1r.Per, TempFit, '-r', 'DisplayName','Regression')
hold off
grid
xlabel('Per')
ylabel('Temp')

The ga funciton came up with several similar results, although some of the parameters (specifically the frequency parameter ‘b(2)’ and the phase parameter, ‘b(3)’) differ. This simply means that there is no unique result for them, which is perfectly acceptable.

I defer to you to interpret these results.

.

[Ashfaq Ahmed]

Dear @Star Strider,

This has been very helpful. Thank you so much. Specially, the idea of creating a uniformly distributed time interval was super clever!!

I have two question just to be confirmed with your algorithm and the relevance of the result:

a) What is the meaning of the peak (amplitude) from the Frequency vs Magnitude graph? I can see the peak magnitude is 4.79 but the final calculation is giving the ampliitude's value as 9.42.

B = [9.42147; 1.051; 3.303; 14.161];

Can you tell me the difference?

b) What is the concept of b(2) and b(3) in the sine wave equation? WHat does the phase and frequency parameter represent in this temeprature time series?

[Star Strider]

As always, my pleasure!

‘Specially, the idea of creating a uniformly distributed time interval was super clever!!’

Thank you! That is actually required for most signal processing, and while it is no longer absolutely necessary to calculate the fft (for example using nufft) it generally makes signal processing easier.

‘What is the meaning of the peak (amplitude) from the Frequency vs Magnitude graph?’

That is the amplitude at that specific frequency. The actual waveform can be made of more than one frequency (and here appears to be), so the apparent amplitude of the time-domain signal can be (and frequently is) different.

‘What is the concept of b(2) and b(3) in the sine wave equation? WHat does the phase and frequency parameter represent in this temeprature time series?’

The ‘b(2)’ parameter is the frequency in cycles/tiime unit. It is multiplied by to produce the radian frequency (radians/tiime unit) that the sin function requires. The ‘b(3)’ parameter is the phase, and while here is multiplied by as well, is actually on the interval of or about 1.9 radians. The phase simply shifts the sin curve to match the waveform. It probably would have been easier not to multiply it by anything, and just estimate it directly.

.

[Ashfaq Ahmed]

Thank you! Now I understand the answer to the second question. But for question (a),

"That is the amplitude at that specific frequency."

Is this the amplitude (at any given cycle) for the entire fourier transformed temperature series? Because we are essentially plotting -

plot(Fv, abs(FT_Temp(Iv))*2)

Also, can you please tell me why was it multiplied with 2?

[Star Strider]

As always, my pleasure!

‘Is this the amplitude (at any given cycle) for the entire fourier transformed temperature series?’

No, just at that specific frequency as described by the ‘bin width’ defining it. A closer look at the frequencies on both sides of the peak frequency shows non-zero components. These are frequencies near enough to the peak frequency to be plotted as part of it and contribute to its amplitude, and what is called ‘spectral leaking’ that is the result of the numerical computation of the fft. (I do not see harmonics at integer multiples of the peak frequency, so it is actually a relatively ‘pure’ waveform, although with some distortion revealed by the width of the peak.)

Note that this:

plot(Fv, abs(FT_Temp(Iv))*2)

plots the entire Fourier transform.

The Fourier transform is computed as a two-sided transform with symmetric negative and positive frequencies, with the total time-domain signal energy divided approximately evenly between the two halves. (This can more easily be seen after using the fftshift function on the original fft result.) In a one-sided Fourier transform, such as I plotted here, multiplying by 2 corrects fotr this, approximating the amplitude of the original signal at those frequencies.

.

[Ashfaq Ahmed]

Hi @Star Strider,

Now this makes much more sense. Brilliant!

I want to comment on the frequency spectra graph that you plotted. The peak is at ~0.1e-3 cycle/hour. So it means after around every 10,000 hours (nearly 365 days), 'something' completes one cycle. And from this data set, the 'something' is the yearly seasonal temperature change. So it means it gets high during the summer and low during winter. Thus, it continues the cycle.

I am trying to find the frequency spectra for tides that are more frequent than this seasonal cycle because it happens twice daily (every 12.25 hours). So I want to catch the frequency at 1/12.25 = 0.08 = 80e-3 cycle/hour range, in more right-side to this plot. I would love your suggestion on how I can do that by not manipulating the time interval but instead using the given time gap. Will the Lomb-Scargle method help?

Footnote: Is it alright if I add you to my paper's acknowledgment section? You have been a fantastic help from the beginning of this project.

[Star Strider]

Thank you!

These temperature data are sampled irregularly, with the most frequent interval being about every 192 hours (every 8 days) and a median interval of 384 hours. The tide data would need to be sampled much more frequently, perhaps every hour, in order to get reliable data from them. Sampling them less frequently about every 3 hours would result in ‘aliasing’ such that the derived frequency would be unreliable and quite likely wrong.

‘Will the Lomb-Scargle method help?’

It easily could, and the Signal Processing Toolbox function plomb implements it. (I did not think to use it here, since this appeared to be a curve-fitting problem, at least as I interpreted it.)

‘Is it alright if I add you to my paper's acknowledgment section?’

I would be honoured!

.

[Ashfaq Ahmed]

Dear @Star Strider,

I have tried to solve the signal with the "plomb" function but could not make much progress, especially keeping these two conditions fixed -

a) the sampling time can not be made uniform. It should be as it is.

b) the frequency spectra should go up to 0.08 cycles/hours to capture the signals coming from the tides (because that happens every 12.5 hours).

I would really appreciate it if you guide me through these two changes.

[PS.: I sent you a message!]

[Star Strider]

The current data are not appropriate for any circadian calculations, since the sampling intervals are not appropriate. They do not have to be regular, however they should ideally be every 1-3 hours. The current sampling intervals can differ by days, and so any tidal variations would be ‘aliased’ and not be correct.

If you have new data for the tides, I would be willing to see if I can help with them.

I just now replied to your mesage! (For the record, it had nothing to do with the content of this thread, so there is no hidden information. I mention this because sometimes people share information offline with one person, and that is not fair to anyone else who might want to help solve a particular problem.)

.