How to remove the periodic oscillations from a signal?

33 次查看(过去 30 天)
The task that I have is to remove the annual and semiannual oscillation from a set of temperature measurements, taken over several years, by means of least squares method.
I found the method described here and since I have to wait for my colleague to first make some analysis on the data and pass them on to me, I decided to make a test run in Matlab, using a signal with two simple sine and cosine functions and random noise.
%%Creating the data
g=2.5;
t=linspace(0,g,g*356*24*60); %n-years of 1min data
o=rand(size(t)); %random noise
p=2*sin(2*pi*t)+1.5*cos(pi*t); %we create two sinusoids
x=p+o; %and add noise to them
%%The method
N=length(x); %no of data
n=(1:N);
T=g*365*24*60; %time that the data covers
f1=1/(365*24*60); %frequency of 1 yr oscillations
f2=1/((365*24*60)/2); %frequency of 1/2 year oscillations
a1=f1*T;
a2=f2*T;
A1=(2/N)*sum(x.*cos((2*pi*a1*n)/N));
A2=(2/N)*sum(x.*cos((2*pi*a2*n)/N));
B1=(2/N)*sum(x.*sin((2*pi*a1*n)/N));
B2=(2/N)*sum(x.*sin((2*pi*a2*n)/N));
C1=sqrt(A1^2+B1^2); %amplitude
C2=sqrt(A2^2+B2^2);
fi1=atan(B1/A1); %phase shift
fi2=atan(B2/A2);
v=(C1*cos(2*pi*t-fi1))+(C2*cos(2*pi*t-fi2)); %the reconstructed set of data
for these two frequencies
r=x-v; %I wasn't sure what was meant by removing, so i just subtracted
%%Visualization
figure %now lets plot the data
plot(t,x)
xlabel('t[min]')
ylabel('Arb. units :)')
title('Data')
figure %and, original signal and noise separately
plot(t,p,'b')
hold on
plot(t,o,'-r')
xlabel('t[min]')
ylabel('Arb. units :)')
title('Data separated')
figure %and once again data (blue), reconstructed signal
plot(t,x,'b-') %(green), and the final result afer subtraction
hold on %(red)
plot(t,v,'g--')
plot(t,r,'r:')
xlabel('t[min]')
ylabel('w/e')
The problem now is, that the reconstructed signal (v) does not resemble the original signal so much. Even when I remove the noise and work only with the pure sinusoid signal, what I get is not really a good representation of it. So, what I'm most interested is:
-Is this at all the correct use of this method? Am I using the correct equations?
-Is there an error somewhere in the code? By error I mostly consider the wrong definition, or calculation of some of the variables.
-In this "ideal case" (sinusoid without noise) should the method be able to reproduce the original curve (almost) exactly?
I hope I've laid out my problem clearly and understandably, if anyone has additional questions, please ask. I know this is rather basic stuff, but my knowledge and experience in this area is very limited (up to non-exsistant), so please bear with me.
Thank you in advance for any help.
  4 个评论
Image Analyst
Image Analyst 2022-12-22
@Tsehaye Negash you are correct. They must have made a typo. They did get it correct later though:
T=g*365*24*60; %time that the data covers

请先登录,再进行评论。

回答(1 个)

Image Analyst
Image Analyst 2013-6-23
编辑:Image Analyst 2013-6-23
Try the Curve Fitting Toolbox. I don't have it so I can't help you much more. Otherwise, try taking the FFT and zeroing out low frequencies until you get a signal without much seasonal oscillations.
You could also try John D'Errico's curve fitting app: SLM

Community Treasure Hunt

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

Start Hunting!

Translated by