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

59 views (last 30 days)

Show older comments

Ashfaq Ahmed
on 10 Mar 2023

Commented: Star Strider
on 15 Mar 2023

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!

##### 5 Comments

Jan
on 10 Mar 2023

"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
on 10 Mar 2023

### Accepted Answer

Star Strider
on 10 Mar 2023

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)]

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

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)

CDOY = DOY+365.25*YRv;

[dlo,dhi] = bounds(CDOY)

dCOD = dhi - dlo

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

DPer = (dhi - dlo)/24

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)

pkfreq = Fv(locs) % cycles/hour

pkper = 1/pkfreq % hours/cycle

pkper_d = pkper/24 % days/cycle

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.

.

##### 8 Comments

Star Strider
on 15 Mar 2023

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

.

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!