MATLAB Answers

I'm calculating the harmonics of a time series. How do I determine the significance of the harmonic fit onto the data?

조회 수: 1(최근 30일)
Ikmal Rosli
Ikmal Rosli 2021년 7월 5일
댓글: Mathieu NOE 2021년 7월 6일
I am calculating the the first few harmonics of a time series. The data and code below is a replication of an example from a textbook from Wilks (2011), pg 371 - 381, Statistical Methods in the Atmospheric Sciences. I've got the amplitude and phase shift right. But I'm not sure how to evaluate the significance of the fit.
clear;clc
% fitting a harmonic
yt = [22.2 22.7 32.2 44.4 54.8 64.3 68.8 67.1 60.2 49.5 39.3 27.4]'; % data
t = (1:12)';
n = 12;
% first harmonic
omegaA1 = cosd(360*t/n);
omegaB1 = sind(360*t/n);
% second harmonic
omegaA2 = cosd(360*t*2/n);
omegaB2 = sind(360*t*2/n);
% regress
b1 = [ones(length(yt), 1) omegaA1 omegaB1 omegaA2 omegaB2]\yt;
% get first harmonic and phase shift;
A1 = hypot(b1(2), b1(3)); % first harmonic amplitude
if b1(2) > 0
phi1 = atand(b1(3)/b1(2))
elseif b1(2) < 0
phi1 = atand(b1(3)/b1(2));
if phi1 < 180
phi1 = phi1 + 180;
elseif phi1 > 180
phi1 = phi1 -180;
end
elseif b(1) == 0
phi1 = 90;
end
% second harmonic and phase shift;
A2 = hypot(b1(4), b1(5)); % second harmonic amplitude
if b1(4) > 0
phi2 = atand(b1(5)/b1(4));
elseif b1(4) < 0
phi2 = atand(b1(5)/b1(4));
if phi2 < 180
phi2 = phi2 + 180;
elseif phi2 > 180
phi2 = phi2 - 180;
end
elseif b1(4) == 0
phi2 = 90
end
% fit and plot
close all
yfit1 = b1(1) + A1*cosd(360*t/n - phi1);
plot(t, yt, 'o')
hold on
plot(t, yfit1)
hold off
[~, p] = corr(yt,yfit1); % I am not quite sure about this part

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 7월 5일
hello
why not simply compute the mean squared error between your fitted data and the input data
I also tried including the second harmonic , but it did not improve the fit...
clearvars;clc
% fitting a harmonic
yt = [22.2 22.7 32.2 44.4 54.8 64.3 68.8 67.1 60.2 49.5 39.3 27.4]'; % data
t = (1:12)';
n = 12;
% first harmonic
omegaA1 = cosd(360*t/n);
omegaB1 = sind(360*t/n);
% second harmonic
omegaA2 = cosd(360*t*2/n);
omegaB2 = sind(360*t*2/n);
% regress
b1 = [ones(length(yt), 1) omegaA1 omegaB1 omegaA2 omegaB2]\yt;
% get first harmonic and phase shift;
A1 = hypot(b1(2), b1(3)); % first harmonic amplitude
if b1(2) > 0
phi1 = atand(b1(3)/b1(2))
elseif b1(2) < 0
phi1 = atand(b1(3)/b1(2));
if phi1 < 180
phi1 = phi1 + 180;
elseif phi1 > 180
phi1 = phi1 -180;
end
elseif b(1) == 0
phi1 = 90;
end
% second harmonic and phase shift;
A2 = hypot(b1(4), b1(5)); % second harmonic amplitude
if b1(4) > 0
phi2 = atand(b1(5)/b1(4));
elseif b1(4) < 0
phi2 = atand(b1(5)/b1(4));
if phi2 < 180
phi2 = phi2 + 180;
elseif phi2 > 180
phi2 = phi2 - 180;
end
elseif b1(4) == 0
phi2 = 90
end
% fit and plot
close all
yfit1 = b1(1) + A1*cosd(360*t/n - phi1);
yfit2 = b1(1) + A1*cosd(360*t/n - phi1)+ A2*cosd(360*t/n - phi2);
mse1 = sqrt(mean((yfit1-yt).^2)) % mean squared error for fit #1
mse2 = sqrt(mean((yfit2-yt).^2)) % mean squared error for fit #2
plot(t, yt, 'o',t, yfit1,t, yfit2)
legend('data','fit 1' ,'fit 2');

추가 답변(0개)

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by