Signal Fitting with sine curve (and how to find out the phase shift?)
71 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Ashfaq Ahmed
am 10 Mär. 2023
Kommentiert: Star Strider
am 15 Mär. 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 Kommentare
Jan
am 10 Mär. 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
am 10 Mär. 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.
Akzeptierte Antwort
Star Strider
am 10 Mär. 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 Kommentare
Star Strider
am 15 Mär. 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.)
.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Parametric Spectral Estimation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!