50 views

Skip to first unread message

Apr 21, 2006, 3:55:57 PM4/21/06

to

A new Matlab user emailed me about applying what amounts to a notch

filter to his data set. He has two years' worth of 30-minute data and

wants to remove daily fluctutations from it. Assuming his daily

fluctuations were sinusoidal, I advised him as follows below.

filter to his data set. He has two years' worth of 30-minute data and

wants to remove daily fluctutations from it. Assuming his daily

fluctuations were sinusoidal, I advised him as follows below.

- Suggestions on my use of the FFT would be appreciated.

- What if the daily fluctuations are regular but not sinusoidal? I

was thinking of segmenting the data into 24-hour portions, average

out the noise, and then calling that the daily fluctuation. Is there

a more sophisticated approach?

Thanks,

Mike

----Original Message Follows----

From: mkleder @ hotmail . com

To:

Subject: for help in Matlab

Date: Fri, 21 Apr 2006

% First, I'm going to create some data to mimic what you

% described as your data stream with the following Matlab

% code:

% create time step variable of two years worth of

% half-hour data points:

t=0:35040;

% create sinusoidal daily cycle every 24 hours

% (every 48 time steps)

x=sin(t/48*2*pi)/2+1;

% superimpose gaussian fluctuations on the data with

% zero mean and stdev=.15:

x=x+randn(size(x))*.15;

% look at the first 3 days of data:

figure

plot(t,x)

xlim([0,48*3])

grid on

xlabel('Elapsed Time in 30-Minute Incrememts')

ylabel('Measurement Value')

% Now the trick is to remove the daily cycles, by which I

% mean sinusoidal wavelengths equal to 24-hours or

% 48 measurements. To do this, I'm first going to take a

% fast fourier transform (FFT) of the data, then I'm going

% to zero-out frequencies corresponding to the daily

% cycle, and then recover the measurements via an

% inverse FFT.

% compute the Fast Fourier Transform

z=fft(x);

% Look at the non-daily flucutations

y=z;

y(731)=0;

y(34312)=0;

xr=real(ifft(y)); % non-daily fluctuations only

figure

plot(t,xr)

xlim([0,48*3])

title('Other Fluctuation Data')

grid on

xlabel('Elapsed Time in 30-Minute Incrememts')

ylabel('Measurement Value')

% Look at the daily fluctuations:

y=0*z;

y(731)=z(731);

y(34312)=z(34312);

xd=real(ifft(y)); % daily fluctuations only

figure

plot(t,xd)

xlim([0,48*3])

title('Daily Sinusoidal Fluctuation Data')

grid on

xlabel('Elapsed Time in 30-Minute Incrememts')

ylabel('Measurement Value')

% So what was done is to remove the daily fluctutation

% from x and put it into xd ("x daily") and the remaining

% fluctuation into xr ("x remainder"). We can check that

% x = xd + xr via:

max(abs(x-(xd+xr)))

% which gives a tiny error due to numerical precision.

% (The "real" in front of the "ifft" function calls is

% present for the same reason: to eliminate any tiny

% imaginary quantites arising from numerical precision

% issues).

% The only remaining question is, how did I know that the

% correct elements in the Fourier transform wre #731 and

% #34312? There are two answers. First, I could calculate

% where those frequency elemtns should be, or second, I

% can just take the FFT of a sine wave of the desired

% frequency imbedded in a data stream of the same length

% as yours and just observe where it's spikes are. (I did

% the latter.)

% Note that changing the PHASE of the daily sinusoid

% doesn't change the process at all. The FFT stores the

% phase as part of the data that is being zeroed out. For

% example, you can change the creation of the original

% daily fluctuation data as follows and get a similar

% result:

x=sin(t/48*2*pi+1)/2+1;

% Keep in mind this trick works best if the daily

% fluctuation data is sinusoidal. If not, then you have

% to move to fancier techniques.

% Hope it helps,

% Mike

----Original Message Follows----

From:

To: mkleder @ hotmail . com

Subject: for help in Matlab

Date: Fri, 21 Apr 2006

hi dear Kleder,

I am a new comer to Matlab. I want to remover the daily cycles

from my heat flux data with filters in Matlab. The time series is two

years long. The dataset is monitored every 30 minutes every day. Can

you give me some suggestions?

Apr 21, 2006, 4:03:14 PM4/21/06

to

Michael Kleder skrev:

> A new Matlab user emailed me about applying what amounts to a notch

> filter to his data set. He has two years' worth of 30-minute data and

> wants to remove daily fluctutations from it. Assuming his daily

> fluctuations were sinusoidal, I advised him as follows below.

>

> - Suggestions on my use of the FFT would be appreciated.

> - What if the daily fluctuations are regular but not sinusoidal? I

> was thinking of segmenting the data into 24-hour portions, average

> out the noise, and then calling that the daily fluctuation. Is there

> a more sophisticated approach?

As far as I can tell from the code you posted, you DFTed the data and

set the relevant coefficients to 0. I suppose it does work, but it

isn't

really an ideal approach. You are touching on one problem: What if

there is some sort of drift in the variation. One solution could be to

DFT the whole data set and set the relevant bins to 0. I can't figure

out off the top of my head how many data points that is, but I

suspect we are talking about tens of thousands of points.

So, what are the alternatives.

A few days ago, I posted this example to get rid of powerline

intereference:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fs = 1000; % Sampling frequency

f = 60; % Location of notch

r = 0.95; % Control parameter. 0 < r < 1.

A = ones(3,1);

B = A;

cW = cos(2*pi*f/fs);

A(2) = -2*r*cW;

A(3) = r^2;

B(2)= -2*cW;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

What this filter does, is that it puts a zero exatly on the unit circle

at the 60 Hz line. The zero removes the line at 60 Hz, but pulls

down the surounding frequencies too. To counter that, a pole

is put just inside the unit circle to push the frequency

response back up.

The notch is put on the nominal location of the interfering line

by choosing the ratio f/fs. The width of the notch is controlled

by the parameter r.

This approach ought to give a little better control and leeway

than setting the spectrum coefficient to 0.

If the variation turns out to be large, one might cosider to use a

band-stop filter or an adaptive filter. It's hard to tell without

seeing the data in question.

Rune

Apr 21, 2006, 4:31:53 PM4/21/06

to

> - What if the daily fluctuations are regular but not sinusoidal?

Any regular (periodic) process can be represented as a series of

harmonically related sinusoids.

> I was thinking of segmenting the data into 24-hour >portions,

average out the noise, and then calling that >the daily

fluctuation. Is there a more sophisticated >approach?

That is essentially low-pass filtering.

Yes, low pass filtering the data. You sample rate is :

30 minutes -> 30/1440 = 0.0208 days

with sampling frequency Fs= 48 samples/day.

So your low-pass filter cutoff will be

Fc= 1 samples/day

So that, normalizing wrt to the sampling frequency we have

wc= 2*pi*(1/48)

This is the cutoff frequency for you digital filter (you can make it

smaller to compensate for the slow filter roll-off). This will filter

any variations within a period of 24 hrs...

In general its not a good idea to simply take the fft of the data and

set those coeffs. to zero beacause you will be getting ringing (a

form of Gibbs effect) in time domain. You will see periodicities that

do not belong to your data.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu