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: 0351
MODEL NAME: NyquistFreq
SHORT DESCRIPTION: A continuous function, a sum of 4 cosines, is sampled at
delta-t =1/2 and undersampled at delta-t=1 illustrating the folding of frequencies around the Nyquist frequency.
*/
import nsrunit;
unit conversion off;
math NyquistFreq {
// INDEPENDENT VARIABLES AND SAMPLING RATES
realDomain t1 sec; t1.min=2; t1.max=128; t1.delta=1;
realDomain t2 sec; t2.min=1; t2.max=128; t2.delta=1/2;
realDomain t3 sec; t3.min=0; t3.max=128; t3.delta=0.1;
// FOUR AMPLITUDES AND FREQUENCIES
real Am1 = 1, Am2=1, Am3 = 1, Am4=1;
real wa1 = (1/8)*PI sec^(-1), wa2 = (3/8)*PI sec^(-1),
wa3 = (11/8)*PI sec^(-1), wa4 = (13/8)*PI sec^(-1);
/* THREE REALIZATIONS OF THE FUNCTION
R1 under sampled
R2 adequately sampled
R3 the "continuous" function discretely realized
*/
real R1(t1) = Am1*cos(wa1*t1)+Am2*cos(wa2*t1)+Am3*cos(wa3*t1)+Am4*cos(wa4*t1);
real R2(t2) = Am1*cos(wa1*t2)+Am2*cos(wa2*t2)+Am3*cos(wa3*t2)+Am4*cos(wa4*t2);
real R3(t3) = Am1*cos(wa1*t3)+Am2*cos(wa2*t3)+Am3*cos(wa3*t3)+Am4*cos(wa4*t3);
// FREQUENCY DOMAIN
// THE FOURIER TRANSFORM OF R1
int kmax1, odd1;
kmax1=floor(t1.ct/2);
odd1=if(2*kmax1=t1.ct) 0 else 1;
realDomain k1; k1.min=1; k1.max = kmax1; k1.delta=1;
real f1(k1), w1(k1);
f1=k1/(t1.ct*t1.delta);
w1=f1*(2*PI);
real A10, A1(k1), B1(k1); // Fourier coefficients
real Periodogram1(k1), Phase1(k1);
// CALCULATE THE FORWARD FOURIER TRANSFORM OF R1
A10 = sum(R1@t1)/t1.ct;
A1(k1)=if(k1<>k1.max+odd1) 2/(t1.ct)*sum((R1*cos(w1*(t1-t1.delta)))@t1)
else sum( ((-1)^(round(t1/t1.delta)-1)*R1(t1))@t1)/t1.ct;
B1(k1)=if(k1<>k1.max+odd1) 2/(t1.ct)*sum((R1*sin(w1*(t1-t1.delta)))@t1)
else 0;
Periodogram1=A1*A1+B1*B1;
Phase1=atan(B1(k1),A1(k1));
// THE FOURIER TRANSFORM OF R2
int kmax2, odd2;
kmax2=floor(t2.ct/2);
odd2=if(2*kmax2=t2.ct) 0 else 1;
realDomain k2; k2.min=1; k2.max = kmax2; k2.delta=1;
real f2(k2), w2(k2);
f2=k2/(t2.ct*t2.delta);
w2=f2*(2*PI);
real A20, A2(k2), B2(k2); // Fourier coefficients
real Periodogram2(k2), Phase2(k2);
// CALCULATE THE FORWARD FOURIER TRANSFORM OF R2
A20 = sum(R2@t2)/t2.ct;
A2(k2)=if(k2<>k2.max+odd2) 2/(t2.ct)*sum((R2*cos(w2*(t2-t2.delta)))@t2)
else sum( ((-1)^(round(t2/t2.delta)-1)*R2(t2))@t2)/t2.ct;
B2(k2)=if(k2<>k2.max+odd2) 2/(t2.ct)*sum((R2*sin(w2*(t2-t2.delta)))@t2)
else 0;
Periodogram2=A2*A2+B2*B2;
Phase2=atan(B2(k2),A2(k2));
// BACK TRANSFORM INTO THE TIME DOMAIN giving R1, R2
real backR1(t1), backR2(t2);
backR1(t1)=A10+sum((A1*cos(w1*(t1-t1.delta))+B1*sin(w1*(t1-t1.delta)))@k1);
backR2(t2)=A20+sum((A2*cos(w2*(t2-t2.delta))+B2*sin(w2*(t2-t2.delta)))@k2);
// LOCATE ACTUAL PEAKS AND FOLDED PEAKS
real peak(k2);
real s = abs(w2(1)-w2(2))/2;
peak = if( abs(w2-wa1)<s) 5 else
if( abs(w2-wa2)<s) 5 else
if( abs(w2-wa3)<s) 5 else
if( abs(w2-wa4)<s) 5 else
if( abs(w2-(2*PI sec^(-1))+wa3)<s) 5 else
-5;
}
/*
DETAILED DESCRIPTION:
The Fourier transform of an undersampled continuous function "folds"
higher frequencies into lower frequencies.
Let series be of length N, and the series samples (V(i), i=1..N), sampled at
(t(i), i=1..N). The intervals between consecutive samples is t.delta. The model is
initialized with an odd number of samples. t.min=t(1). t.delta is the
interval between samples. t.ct, the number of samples, equals N.
The maximum number of frequencies is
kmax = N/2, N even,
= (N-1)/2, N odd.
The frequencies are f(i) = i/(N*t.delta), i=1..kmax.
The highest frequency is 1/(2*t.delta) which is called the Nyquist frequency,
also called the folding frequency. In radians/sec it is 2*PI/(2*t.delta).
Let
odd=0 , N is even and
odd=1, N is odd.
The forward transform is calculated as
N
A0 = Sum(V(i)/N).
i=1
for j=1 to kmax,
IF ( j not equal to kmax+odd),
N
A(j)=2/N*sum(V(i)*cos(2*PI*f(j)*(t(i)-t.min-t.delta))
i=1
N
B(j)=2/N*sum(V(i)*sin(2*PI*f(j)*(t(i)-t.min-t.delta))
i=1
ELSE
N
A(j)=sum( (-1)^(round((t(i)-t.min)/t.delta)-1)*V(i)/N).
i-1
B(j)=0.
END
END
The Periodogram estimate P(j) = A(j)*A(j)+B(j)*B(j).
The back transform is given as
for i = 1 to N
kmax
V(i)=A0+sum ( A(j)*cos(2*PI*f(j)*(t(i)-t.min-t.delta)) + B(j)*sin(2*PI*f(j)*(t(i)-t.min-t.delta)) )
j=1
The original series was built using four frequencies,
[ PI/8, 3PI/8, 11PI/8 and 13PI/8 radians] and unit amplitudes.
The periodograms (lower plot) are plotted in blue for R1 (under sampled)
and R2 (adequately sampled).
R2 with a deltaT=1/2 has a cut off frequency of 2*PI radians, and hence
the periodogram for R2 captures this information.
R1 with a deltaT=1 has a cut of frequency of PI radians. Therefore
the last frequency of 13PI/8 is folded to PI-(13PI/8-PI) = 3PI/8 and the
peak for the periodogram is boosted from 1 to 4 (1+1)^2=4. The peak
that was at 11PI/8 appears as a new peak at PI-(11PI/8)-PI) = 5PI/8.
Notice how the peaks are folded around the Nyquist frequency. The
decaying envelope around each frequency peak is because the frequencies
calculated do not match the frequencies used in the function. By shifting
the starting times of the samples (Parameter set Sharp) we can make the
peaks appear at the exact frequencies.
KEY WORDS:
time series, Nyquist, folding, Fourier, undersampling, frequencies
REFERENCES:
None.
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
*/
Read the DETAILED DESCRIPTION following the program.