How to estimate annual and semi-annual signals?

4 visualizaciones (últimos 30 días)
hanif hamden
hanif hamden el 25 de Nov. de 2020
Comentada: Mathieu NOE el 26 de Nov. de 2020
Hi everyone, I would like to estimate annual ans semi-annual signal from my data. Attached is my data as well as my coding, I managed to estimate bias and linear trend. However, I was stuck in estimating annual and semi-annual signal. I hope anyone can help me on this matter. Thank you.
A = load('FPT1.txt');
t = A(:,1);
slv = A(:,5);
dt1 = num2str(t,'%.2f');
dtm1 = datetime(dt1,'InputFormat','yyyyMMddHHmmss.SS');
dtn1 = datenum(dtm1);
dtmFill1 = fillmissing(dtm1,'linear');
%%PARAMETER FIT
% hobs = h0 + h1T + h2*sin(2*pi*T) + h3*cos(2*pi*T) + h4*sin(4*pi*T) + h5*cos(4*pi*T)
% Estimate bias/offset
h0 = mean(slv);
% Estimate linear trend
b = polyfit(dtn1,slv,1);
h1 = polyval(b, dtn1);
%% Stuck here
% Estimate Annual Signals [h2*sin(2*pi*T) + h3*cos(2*pi*T)]
% However, the time for this data is every 9.9156 days
% f1=1/(365*24); %frequency of 1 yr oscillations
h2 = sin(2*pi*T)/slv;
h3 = cos(2*pi*T)/slv;
% Estimate Semi-Annual Signals [h4*sin(4*pi*T) + h5*cos(4*pi*T)]
% f2=1/((365*24))/2; %frequency of 1/2 year oscillations
h4 = slv*sin(4*pi*T);
h5 = slv*cos(4*pi*T);

Respuesta aceptada

Mathieu NOE
Mathieu NOE el 25 de Nov. de 2020
hello
you're not stuck anymore
below my code (the h parameters computation is the same as DFT coefficients)
if you're stisfied, pls accept my answer
all the best
A = load('FPT1.txt');
t = A(:,1);
slv = A(:,5);
dt1 = num2str(t,'%.2f');
dtm1 = datetime(dt1,'InputFormat','yyyyMMddHHmmss.SS');
dtn1 = datenum(dtm1);
% resample with constant time sampling
% remember time is in days
dt = 7; % 7 days interval
Fs = 1/dt; % NB not in Hz = 1/s but in 1/day
new_t = min(dtn1):dt:max(dtn1);
new_data = interp1(dtn1,slv,new_t);
% figure(1),plot(dtn1,slv,'b+-',new_t,new_data,'r*-');
%%PARAMETER FIT
% hobs = h0 + h1T + h2*sin(2*pi*T) + h3*cos(2*pi*T) + h4*sin(4*pi*T) + h5*cos(4*pi*T)
% Estimate bias/offset
% h0 = mean(new_data);
% Estimate linear trend
b = polyfit(new_t,new_data,1);
% h1 = polyval(b, new_t);
%% Not anymore Stuck here
% Estimate Annual Signals [h2*sin(2*pi*T) + h3*cos(2*pi*T)]
% However, the time for this data is every 9.9156 days
f1=1/(365); %frequency of 1 yr oscillations time base in days)
h2 = 2/length(new_t)*sum(sin(2*pi*f1*new_t).*new_data);
h3 = 2/length(new_t)*sum(cos(2*pi*f1*new_t).*new_data);
% Estimate Semi-Annual Signals [h4*sin(4*pi*T) + h5*cos(4*pi*T)]
f2=2*f1; %frequency of 1/2 year oscillations
h4 = 2/length(new_t)*sum(sin(2*pi*f2*new_t).*new_data);
h5 = 2/length(new_t)*sum(cos(2*pi*f2*new_t).*new_data);
% reconstructed signal
hobs = h2*sin(2*pi*f1*new_t) + h3*cos(2*pi*f1*new_t) + h4*sin(2*pi*f2*new_t) + h5*cos(2*pi*f2*new_t);
% use the delta = new_data-hobs, for better linear trend estimation
delta = new_data-hobs;
% Estimate linear trend
b = polyfit(new_t,delta,1);
h0 = b(2);
h1 = b(1);
hobs = h0 + h1*new_t + h2*sin(2*pi*f1*new_t) + h3*cos(2*pi*f1*new_t) + h4*sin(2*pi*f2*new_t) + h5*cos(2*pi*f2*new_t);
figure(2),plot(dtn1,slv,'b+-',new_t,hobs,'r-*');
  2 comentarios
hanif hamden
hanif hamden el 26 de Nov. de 2020
Thank you sir. the code works well. Just one more question. I wanted to subtract the hobs model with the original data, since the matrix is not the same. Is it correct if I interpolate the hobs based on time?
% -----Additional---
new_t2 = new_t';
hobs2 = hobs';
dtm2 = datetime(new_t2,'ConvertFrom','datenum');
Ts_hobs = [new_t2 hobs2];
new_Ts_hobs = interp1(Ts_hobs(:,1),Ts_hobs(:,2),dtn1,'makima');
Nslv = slv-new_Ts_hobs;
figure(3),plot(dtn1,slv,'b+-',dtn1,Nslv,'r-*');
Mathieu NOE
Mathieu NOE el 26 de Nov. de 2020
yes it's ok
just simplified a bit your code;
this is doing the same with less :
% -----Additional---
dtm2 = datetime(new_t,'ConvertFrom','datenum');
new_Ts_hobs = interp1(new_t,hobs,dtn1,'makima');
Nslv = slv-new_Ts_hobs;
figure(3),plot(dtn1,slv,'b+-',dtn1,Nslv,'r-*');

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Linear Model Identification en Help Center y File Exchange.

Productos


Versión

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by