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

59 views (last 30 days)
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
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
John D'Errico on 10 Mar 2023
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.

Sign in to comment.

Accepted Answer

Star Strider
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)]
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.
.
  8 Comments
Ashfaq Ahmed
Ashfaq Ahmed on 15 Mar 2023
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
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.)
.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!