81 views (last 30 days)

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.

Image Analyst
on 23 Jun 2013

Edited: Image Analyst
on 23 Jun 2013

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.