notch filter, sinusoidal vs. non-sinusoidal
Michael Kleder
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?
Rune Allnor
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
Dapple
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.