the SFT (slow Fourier Transform)
Number of points can be odd or even. Program corrects last
Fourier coefficients so that
if
ValueFirstPoint NOT EQUAL ValueLastPoint,
the back Transform returns the correct value
instead of the average of the first and
last values.
If you don't need the back transform,
comment it out.
/* MODEL NUMBER: 0352
MODEL NAME: TimeSeriesFilters
SHORT DESCRIPTION: Simple low pass, high pass, and bandpass
filters are illustrated by four examples.
*/
source procedure lopass(nlopass, length, input@t; output@t) {
// do nlopass averages over length number of points
// wrapping the data at the endpoints.
language="java";
maincode={{
RegularGridData t = (RegularGridData) input.grid(0);
int n = t.ct();
int NLO = (int)nlopass.realVal();
double hold[] = new double[n];
double avg[] = new double[n];
for( int i=0; i<n; i++) {
hold[i] = input.realVal(i);
}
for (int ilo =1; ilo<=NLO; ilo++) {
int L = (int)length.realVal();
for (int i=-L/2; i<n-L/2; i++) {
int inew = i+L/2;
avg[inew]=0;
for (int j=1; j<=L; j++) {
int k=i+j;
if(k<0) {k=k+n;}
if(k>=n) {k=k-n;}
avg[inew]+=hold[k];
}
avg[inew]=avg[inew]/(double)L;
}
for (int i=0; i<n; i++) {
hold[i]=avg[i];
}
}
for (int i=0; i<n; i++) {
double out = hold[i];
output.set(i,out);
}
}};
}
source procedure hipass(nhi, input@t; output@t) {
// do n two point differences on data wrapping back to first point
language="java";
maincode={{
RegularGridData t = (RegularGridData) input.grid(0);
int n = t.ct();
double ave[] = new double[n];
double hold[] = new double[n];
for (int i=0; i<n; i++) {
hold[i]=input.realVal(i);
}
int NHI = (int)nhi.realVal();
if(NHI>0) {
for (int jhi=1; jhi<=NHI; jhi++) {
double save=hold[0];
for (int i=0; i<n; i++) {
if(i==n-1) {ave[i] = -hold[i]/2+save/2;
} else {
ave[i] = -hold[i]/2+hold[i+1]/2;
}
hold[i] = ave[i];
}
}
}
for (int i=0; i<n; i++) {
double out = hold[i];
output.set(i,out);
}
}};
}
import nsrunit;
unit conversion off;
math TimeSeriesFilters {
real nlopass=1;
real Lengthlopass = 10;
real nhipass=0;
real L=nlopass+nhipass;
// TIME DOMAIN
realDomain t sec; t.min=0; t.max=63; t.delta=1;
real R1(t); // R1 is uniform random noise between
// -1 and 1
R1= 2*(random()-0.5);
real R3(t); // R1 convolved with low pass filter becomes R3
real R2(t); // R3 convolved with high pass filter becomes R2
lopass(nlopass,Lengthlopass,R1@t,R3@t);
hipass(nhipass,R3@t,R2@t);
//real ts = 1 sec;
real bR(t); // The R1 series from the back transform;
// FREQUENCY DOMAIN
/* Check if series has an odd or even number of points and set
add to modify calculation of last two Fourier Estimates. */
int kmax, odd;
kmax=floor(t.ct/2);
odd=if(2*kmax=t.ct) 0 else 1;
realDomain k; k.min=1; k.max = kmax; k.delta=1;
real f(k), w(k);
f=k/(t.ct*t.delta);
w=f*(2*PI);
real A10, A1(k), B1(k); // Fourier coefficients
real Periodogram1(k), Phase1(k);
real A20, A2(k), B2(k); // Fourier coefficients
real Periodogram2(k), Phase2(k);
/* NOTA BENE: From the Fourier transform of a time series
the periodogram can be calculated. The periodogram is
an "estimate of the spectrum," NOT the spectrum.
It is a biased estimator. */
// CALCULATE THE FORWARD FOURIER TRANSFORM OF R1
A10 = sum(R1@t)/t.ct;
A1(k)=if(k<>k.max+odd) 2/(t.ct)*sum((R1*cos(w*(t-t.min-t.delta)))@t)
else sum( ((-1)^(round((t-t.min)/t.delta)-1)*R1(t))@t)/t.ct;
B1(k)=if(k<>k.max+odd) 2/(t.ct)*sum((R1*sin(w*(t-t.min-t.delta)))@t)
else 0;
Periodogram1=A1*A1+B1*B1;
Phase1=atan(B1(k),A1(k));
// CALCULATE THE FORWARD FOURIER TRANSFORM OF R2
A20 = sum(R2@t)/t.ct;
A2(k)=if(k<>k.max+odd) 2/(t.ct)*sum((R2*cos(w*(t-t.min-t.delta)))@t)
else sum( ((-1)^(round((t-t.min)/t.delta)-1)*R2(t))@t)/t.ct;
B2(k)=if(k<>k.max+odd) 2/(t.ct)*sum((R2*sin(w*(t-t.min-t.delta)))@t)
else 0;
Periodogram2=A2*A2+B2*B2;
Phase2=atan(B2(k),A2(k));
// CALCULATE THE FOURIER TRANSFORM OF THE BANDPASS FILTER
real Filter(k);
Filter=Periodogram2/Periodogram1;
// BACK TRANSFORM INTO THE TIME DOMAIN giving R2
bR(t)=A20+sum((A2*cos(w*(t-t.delta))+B2*sin(w*(t-t.delta)))@k);
/*real R4(t);
R4 = R1-R2;
real A30, A3(k), B3(k); // Fourier coefficients
real Periodogram3(k), Phase3(k);
// CALCULATE THE FORWARD FOURIER TRANSFORM OF R2
A30 = sum(R4@t)/t.ct;
A3(k)=if(k<>k.max+odd) 2/(t.ct)*sum((R4*cos(w*(t-t.min-t.delta)))@t)
else sum( ((-1)^(round((t-t.min)/t.delta)-1)*R4(t))@t)/t.ct;
B3(k)=if(k<>k.max+odd) 2/(t.ct)*sum((R4*sin(w*(t-t.min-t.delta)))@t)
else 0;
Periodogram3=A3*A3+B3*B3;
Phase3=atan(B3(k),A3(k));*/
}
/*
DETAILED DESCRIPTION:
This program convolves low pass and high pass filters with a time series
of Gaussian white noise. The ends of the series are "wrapped."
An example will clarify this.
Let A3 be the filter (1/3, 1/3, 1/3). We will convolve it with series
Y=(y0, y1,y2,y3,y4,y5):
A3(t) convolved with Y(t)=H(t) =[ (y5+y0+y1)/3, (y0+y1+y2)/3, ... (y4+y5+y0)/3 ].
Note how y5 and y0 are wrapped. It is as if the series was
[y5, y0, y1, y2, y3, y4, y5, y0].
A high pass filter can be generated by differencing adjacent points.
By combining low and high pass filters, bandpass filters can be constructed.
Filtering operations can introduce phase shifts in the estimated Fourier
coefficients.
A Fourier transform of a time series is a least squares fit of the
time series with a set of sine and cosine functions. Since the chosen
frequencies used in the transform are orthogonal functions, their coefficients
can be found indpendently.
The periodogram at a give frequency is constructed by adding the squares
of the Fourier coefficients, i.e.,
P(f)=A(f)*A(f)+B(f)*B(f).
This is called a Power Spectrum and is usually plotted log(P) vs. log(f).
It is important to note that the Power Spectrum is not the spectrum of
the time series. It is a biased estimate of the spectrum.
The Fourier transform can be used to compute the periodogram of the
original series, P1, and the filtered series, P2. The ratio P2/P1
is the periodogram of the Filter.
In Example 1, a series of 365 Gaussian random numbers has been
generated and low passed with a 7 point low pass filter. It should be noted
that although the periodograms of the time series before and after
filtering are noisy, the frequency response of the filter is smooth.
In Example 2, the same series is low passed twice with the 7 point
filter. Note that the Frequency response of the Filter appears to be
the same, but notice the scaling of the amplitude--it doubles.
In Example 3, the same series is high passed with a two point difference
filter.
In Example 4, the same series is band passed.
These are simple examples and more complex kinds of filters have been
designed.
KEY WORDS:
time series, filter, filtering,low pass, high pass, bandpass, band pass,
Fourier
REFERENCES:
http://en.wikipedia.org/wiki/Filter_design
http://en.wikipedia.org/wiki/Window_function
http://en.wikipedia.org/wiki/Fourier_transform
REVISION HISTORY:
JSim SOFTWARE COPYRIGHT AND REQUEST FOR ACKNOWLEDGMENT OF USE:
JSim software was developed with support from NIH grants HL088516, and HL073598. Please cite these grants in any publication for which
this software is used and send one reprint of published abstracts or
articles to the address given below. Academic use is unrestricted.
Software may be copied so long as this copyright notice is included.
Copyright (C) 1999-2012 University of Washington.
Contact Information:
The National Simulation Resource,
Director J. B. Bassingthwaighte,
Department of Bioengineering,
University of Washington, Seattle, WA
98195-5061
*/
Plot Page: TimeSeries
Series R1 is generated and plotted. It is then low-passed nlopass times
and hi-passed nhipass times. The resulting series is R2, and it
is plotted with a horizontal offset. This is more apparent when
the number of points, t.max is set to a smaller value (try 15).
This model does not use a Fast Fourier Transform (FFT) routine. For
short series (<2,000 points, it is fast enough for our purposes.