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

Question

Ashfaq Ahmed 2023-3-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1926585-signal-fitting-with-sine-curve-and-how-to-find-out-the-phase-shift

评论： Star Strider 2023-3-15

采纳的回答： Star Strider

SignalFitting.mat

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 个评论
显示 3更早的评论隐藏 3更早的评论

Jan 2023-3-10

在 MATLAB Online 中打开

"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 2023-3-10

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 2023-3-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1926585-signal-fitting-with-sine-curve-and-how-to-find-out-the-phase-shift#answer_1190580

在 MATLAB Online 中打开

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 个评论
显示 6更早的评论隐藏 6更早的评论

Star Strider 2023-3-13

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.

.

Star Strider 2023-3-13

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 2023-3-13

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.

Ashfaq Ahmed 2023-3-15

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 2023-3-15

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

.

请先登录，再进行评论。

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

5 个评论 显示 3更早的评论隐藏 3更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论