In [1]:
# edapy_06_00 clear all variables and import various modules
# History
# 2022/06/26 - checked with most recent Python & etc. software
# 2023/07/13 - fixed some issues, incl colormap, loglof plt, eda_cvec()
# 2024/05/26 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched dot product to return scalar value
%reset -f
import os
from datetime import date
from math import exp, pi, sin, sqrt, floor, ceil
import numpy as np
import scipy.linalg as la
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.colors import ListedColormap
# eda_draw function makes a "pictorial matrix equation"
# arguments are vectors, matrices and strings
# which are plotted in the order that the appear
# except that strings starting with 'title ' are plotted
# under the subseqeunt matrix or vector
# always returns a status of 1
def eda_draw(*argv):
bw = np.zeros((256,4));
v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
bw[:,0] = v;
bw[:,1] = v;
bw[:,2] = v;
bw[:,3] = np.ones(256);
bwcmap = ListedColormap(bw);
# size of plot
W = 16;
H = 4;
fig1 = plt.figure(1);
# figsize width and height in inches
fig1.set_size_inches(W,H);
ax1 = plt.subplot(1,1,1);
plt.axis([0, W, -H/2, H/2]);
plt.axis('off');
LM = W/6; # matrix width and heoght
LV = W/40; # vector width
FS = 0.12; # character width
TO = 0.4; # title vertical offset
SP = 0.2; # space between objects
LS = 0.2; # leading space
p = LS; # starting x-position
istitle=0; # flags presence of a title
for a in argv:
if isinstance(a,np.ndarray):
sh = np.shape(a);
if len(sh) == 1: # conversion to nx1 array
n = sh[0];
m = 1;
ap = a;
a = np.zeros((n,1));
a[:,0] = ap;
else:
n = sh[0];
m = sh[1];
if m==1:
pold=p;
left=p;
right=p+LV;
bottom=-LM/2;
top=LM/2;
plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
p = p+LV;
pm = (p+pold)/2;
if istitle:
plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
istitle=0;
p = p+SP;
else:
pold=p;
left=p;
right=p+LM;
bottom=-LM/2;
top=LM/2;
plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
p = p+LM;
pm = (p+pold)/2;
if istitle:
plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
istitle=0;
p = p+SP;
elif isinstance(a,str):
ns = len(a);
istitle=0;
if( ns>=6 ):
if 'title ' in a[0:6]:
istitle=1;
titlestr=a[6:];
if( istitle != 1):
plt.text(p,0,a);
p = p + ns*FS + SP;
plt.show();
return 1;
# function to make a numpy N-by-1 column vector
# c=eda_cvec(v1, v2, ...) from a list of several
# array-like entities v1, v2, including a number
# a list of numbers, a tuple of numbers, an N-by-0 np array
# and a N-by-1 np array. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it
# It promotes integers to floats, and integers
# and floats to complex, by context.
# This version concatenates many argments,
# whereas c=eda_cvec1(v1) takes just one argiment.
# I recommend always using eda_cvec(v1, v2, ...)
def eda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = eda_cvec1(a);
N,M = np.shape(v);
Nt = Nt + N;
if( N==0 ):
continue; # skip vector of zero length
if (t==int) and isinstance(v[0,0],float):
t=float;
elif isinstance(v[0,0],complex):
t=complex;
w = np.zeros((Nt,1),dtype=t);
Nt = 0;
for a in argv:
v = eda_cvec1(a);
N,M = np.shape(v);
w[Nt:Nt+N,0] = v[0:N,0];
Nt = Nt + N;
return w;
# function to make a numpy N-by-1 column vector
# c=eda_cvec1(v) from entity v that is array-like,
# including a number, a list of numbers, a tuple
# of numbers, an N-by-0 np array and a N-by1 np array.
# It promotes integers to floats, and integers
# and floats to complex, by context. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it.
# This version takes just one input argment.
# whereas c=eda_cvec(v1,v2,...) concatenates
# many argiments.
def eda_cvec1(v):
if isinstance(v, int) or isinstance(v, np.int32):
w = np.zeros((1,1),dtype=int);
w[0,0] = v;
return w;
elif isinstance(v, float):
w = np.zeros((1,1),dtype=float);
w[0,0] = v;
return w;
elif isinstance(v, complex):
w = np.zeros((1,1),dtype=complex);
w[0,0] = v;
return w;
elif isinstance(v, np.ndarray):
s = np.shape(v);
if len(s) == 1:
return np.copy(np.reshape(v,(s[0],1)));
else:
[r,c]=s;
if( c==1 ):
return(np.copy(v));
elif(r==1):
return(np.copy(v.T));
else:
raise TypeError("eda_cvec: %d by %d ndarray not allowed" % (r, c));
elif isinstance(v, list):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
else:
raise TypeError("eda_cvec: %s not supported" % type(v));
In [2]:
# edapy_06_01
# cosine wave
# make the auxillary variable t
N=100;
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = eda_cvec( np.linspace(tmin,tmax,N) );
# a cosine of ampmitude C, time offset t0 and period T
t0 = 10.0;
T = 50;
C= 2.0;
d = C*np.cos( 2*pi*(t-t0)/T );
# plot data
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(1,1,1);
plt.axis([tmin, tmax, 1.1*np.min(d), 1.1*np.max(d)])
plt.plot(t,d,'k-');
plt.xlabel('time t');
plt.ylabel('d(t)');
plt.title('Delayed cosine wave');
plt.show();
In [3]:
# edapy_06_02
# plot columns of G
# (not normalized - see eda06_13 for properly normalized version)
# square data kernel
N=32; # must be an even number
M=N;
# make the auxillary variable t
N=100;
Dt = 1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
Dt = (tmax - tmin)/(N-1);
t = eda_cvec( np.linspace(tmin,tmax,N) );
# standard definition of nyquist frequency fmax
# frequency increment Df, frequency vector f,
# angular frequency vector w
Nf=floor(N/2+1);
fmax = 1/(2*Dt);
Df = fmax/(N/2);
f = eda_cvec( np.linspace(0,fmax,Nf) );
Nw=Nf;
wmax = 2*pi*fmax;
Dw = wmax/(N/2);
w = eda_cvec( np.linspace(0,wmax,Nw) );
# set up data kernel G
G = np.zeros((N,M));
# zero frequency column
G[:,0] = np.ones((N,1)).ravel();
# interior M/2-1 columns. Note indice arithmetic
# to keep track of pairs of columns, one for cos(),
# one for sin()
Mo2 = floor(M/2);
for i in range(1,Mo2):
j = 2*i-1;
k = j+1;
G[:,j]=np.cos(w[i,0]*t).ravel();
G[:,k]=np.sin(w[i,0]*t).ravel();
# nyquist column
G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();
# plot selected columns of dat kernel G
fig1 = plt.figure(1,figsize=(8,4));
plt.rcParams['xtick.bottom'] = False;
plt.rcParams['xtick.labelbottom'] = False;
plt.rcParams['xtick.top'] = True;
plt.rcParams['xtick.labeltop'] = True;
ip = 1;
for i in [0,1,2,3,4,M-1]:
ax = plt.subplot(1,6,ip);
ax.invert_yaxis();
plt.axis([-2, 2, tmin, tmax ]);
plt.plot(G[:,i],t,'k-');
plt.plot(G[:,i],t,'k.');
mystr = "col %d" % (i+1);
plt.title(mystr);
if (i==0):
plt.ylabel('time t');
ip=ip+1;
plt.show();
In [4]:
# edapy_06_03: example of aliasing
# make the auxillary variable t
N=100;
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
Dt1 = (tmax - tmin)/(N-1);
t = eda_cvec( np.linspace(tmin,tmax,N) );
# standard set-up of Fourier parameters
# note: N required to be an even integer
Nf=floor(N/2+1); # number of non-negative frequencies
fmax = 1/(2*Dt1); # Nyquist frequency
Df = fmax/(N/2); # frequency increment
# non-negative frequency axis
f = eda_cvec( np.linspace(0,fmax,Nf) );
Nw=Nf; # number of angular frequencie
wmax = 2*pi*fmax; # Nyquist angukar frequency
Dw = wmax/(N/2); # angular frequency incremnt
# non-negative angular frequency axis
w = eda_cvec( np.linspace(0,wmax,Nw) );
# low frequency oscillation
w1= 2*Dw;
print('low frequency', w1);
d1 = np.cos( w1*t );
# high frequency oscillation
w2= (2+N)*Dw;
print('high frequency', w2);
d2 = np.cos( w2*t );
# set up higher sample rate time
N=500;
Dt2=0.1;
tmin=0.0;
tmax = tmin+Dt2*(N-1);
t3 = eda_cvec( np.linspace(tmin,tmax,N) );
d3 = np.cos( w2*t3 );
# plot signals
fig1 = plt.figure(1,figsize=(16,4));
plt.rcParams['xtick.bottom'] = False;
plt.rcParams['xtick.labelbottom'] = False;
plt.rcParams['xtick.top'] = True;
plt.rcParams['xtick.labeltop'] = True;
ax1 = plt.subplot(2,1,1);
plt.axis([tmin, tmax, -2, 2 ]);
plt.plot(t,d1,'k-');
plt.plot(t,d2,'ko');
plt.xlabel('time t');
plt.ylabel('d(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([tmin, tmax, -2, 2 ]);
plt.plot(t,d2,'k-');
plt.plot(t,d2,'ko');
plt.plot(t3,d3,'k:');
plt.xlabel('time t');
plt.ylabel('d(t)');
plt.show();
low frequency 0.12566370614359174 high frequency 6.408849013323178
In [5]:
# edapy_06_04: another example of aliasing
# make the auxillary variable t
N=50;
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
Dt1 = (tmax - tmin)/(N-1);
t = eda_cvec( np.linspace(tmin,tmax,N) );
# standard set-up of Fourier parameters
# note: N required to be an even integer
Nf=floor(N/2+1); # number of non-negative frequencies
fmax = 1/(2*Dt1); # Nyquist frequency
Df = fmax/(N/2); # frequency increment
# non-negative frequency axis
f = eda_cvec( np.linspace(0,fmax,Nf) );
Nw=Nf; # number of angular frequencie
wmax = 2*pi*fmax; # Nyquist angukar frequency
Dw = wmax/(N/2); # angular frequency incremnt
# non-negative angular frequency axis
w = eda_cvec( np.linspace(0,wmax,Nw) );
# low frequency oscillation
w1= -2*Dw;
print('low frequency', w1);
d1 = np.sin( w1*t );
# high frequency oscillation
w2= (N-2)*Dw;
print('high frequency', w2);
d2 = np.sin( w2*t );
# set up higher sample rate time
N=500;
Dt2=0.1;
tmin=0.0;
tmax = tmin+Dt2*(N-1);
t3 = eda_cvec( np.linspace(tmin,tmax,N) );
d3 = np.cos( w2*t3 );
# plot signals
fig1 = plt.figure(1,figsize=(8,4));
plt.rcParams['xtick.bottom'] = False;
plt.rcParams['xtick.labelbottom'] = False;
plt.rcParams['xtick.top'] = True;
plt.rcParams['xtick.labeltop'] = True;
ax1 = plt.subplot(2,1,1);
plt.axis([tmin, tmax, -2, 2 ]);
plt.plot(t,d1,'k-');
plt.plot(t,d2,'ko');
plt.xlabel('time t');
plt.ylabel('d(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([tmin, tmax, -2, 2 ]);
plt.plot(t,d2,'k-');
plt.plot(t,d2,'ko');
plt.plot(t3,d3,'k:');
plt.xlabel('time t');
plt.ylabel('d(t)');
plt.show();
low frequency -0.25132741228718347 high frequency 6.031857894892403
In [6]:
# edapy_06_05: a.s.d. of the Neuse River hydrograph
# using least-squares
# (not normalized - see edapy_06_14 for properly normalized version)
# load data from file
D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
[Nraw, K]=D.shape;
traw = eda_cvec( D[:,0] ); # traw for "t raw"
Dt = traw[1,0]-traw[0,0]; # sampling interval
draw = eda_cvec( D[:,1] ); # draw for "d raw"
# round off to even number of points
No2 = floor(Nraw/2);
N = 2*No2;
d = eda_cvec( draw[0:N,0] );
t = eda_cvec( traw[0:N,0] );
# plot data
fig1 = plt.figure(1,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);
plt.plot(t,d,'k-');
plt.xlabel('time (days)');
plt.ylabel('discharge (cfs)');
plt.title('Neuse River Hydrograph');
plt.show();
# nyquist frequencies
Nf=floor(N/2+1);
fmax = 1/(2*Dt);
Df = fmax/(N/2);
f = eda_cvec(np.linspace(0,fmax,Nf));
Nw=Nf;
wmax = 2*pi*fmax;
Dw = wmax/(N/2);
w = eda_cvec(np.linspace(0,wmax,Nw));
fpos = eda_cvec(Df*np.linspace(0,No2,Nf));
wpos = eda_cvec(Dw*np.linspace(0,No2,Nf));
# set up data kernel G
M=N;
Mo2=floor(M/2);
G = np.zeros((N,M));
G[0:N,0:1] = np.ones((N,1)); # zero frequency column
# interior M/2-1 columns.
for i in range(1,Mo2):
j = 2*i-1;
k = j+1;
G[0:N,j:j+1]=np.cos(w[i,0]*t);
G[0:N,k:k+1]=np.sin(w[i,0]*t);
G[0:N,M-1:M] = np.cos(w[Nw-1,0]*t); # nyquist column
# solve least squares problem using
# analytic formula for inv(GT*G)
GTd = np.matmul( G.T, d );
gtgi = 2*np.ones((M,1))/N;
gtgi[0,0]=1/N;
gtgi[M-1,0]=1/N;
mest = np.multiply( gtgi, GTd );
# compute power spectral density s
s=np.zeros((Nw,1));
s[0,0]=np.sqrt(mest[0,0]**2); # zero frequency
for i in range(1,Mo2): # interior points
j = 2*i-1;
k = j+1;
s[i,0] = sqrt(mest[j,0]**2 + mest[k,0]**2);
s[Nw-1,0]= sqrt(mest[M-1,0]**2); # Nyquist frequency
# plot spectrum
fig2 = plt.figure(2,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,3,1);
plt.axis([0, fmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
plt.plot(fpos[1:Nf,0],s[1:Nf,0],'k-');
plt.xlabel('frequency f (cycles per day)');
plt.ylabel('amplitude spectrum');
plt.title('linear in frequency');
ax1 = plt.subplot(1,3,2);
Tmax = 1000;
plt.axis([0, Tmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
plt.plot(np.reciprocal(fpos[1:Nf,0]),s[1:Nf,0],'k-');
plt.xlabel('periot T (days)');
plt.ylabel('amplitude spectrum');
plt.title('linear in period');
ax1 = plt.subplot(1,3,3);
plt.semilogx(fpos[1:Nf,0], s[1:Nf,0], 'k-');
plt.axis([0.001, 1, 0, 1.1*np.max(s[1:Nf,0])])
plt.xlabel('frequency f (cycles per day)');
plt.ylabel('amplitude spectrum');
plt.title('logarithmic in frequency');
plt.show();
# check results
dpre =np.matmul( G, mest );
e = d - dpre; # individual errors
E = np.matmul( e.T, e ); E=E[0,0]; # total error
print('Error', E);
Error [[2.87461543e-15]]
In [7]:
# edapy_06_06: a.s.d. of the Neuse River hydrograph
# using fft command
# (unnormalized spectral density - see eda06_14 for proper normalization)
# load data from file
D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
[Nraw, K]=D.shape;
traw = eda_cvec( D[:,0] ); # traw for "t raw"
Dt = traw[1,0]-traw[0,0]; # sampling interval
draw = eda_cvec( D[:,1] ); # draw fro "d raw"
# round off to even number of points
No2 = floor(Nraw/2);
N=2*No2;
d = eda_cvec( draw[0:N,0] ); # truncate data vector
t = eda_cvec( traw[0:N,0] ); # truncate time vector
# plot data
fig1 = plt.figure(1,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);
plt.plot(t,d,'k-');
plt.xlabel('time (days)');
plt.ylabel('discharge (cfs)');
plt.title('Neuse River Hydrograph');
plt.show();
No2 = floor(N/2); # N/2 rounded off to an integer
M=N; # number of model parameters equals number of data
# set up time axis
tmin = 0;
tmax = Dt*(N-1);
t= eda_cvec( np.linspace(tmin,tmax,N) );
# generic frequency setup
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# f is full freqeuncy vector:
# f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
f=eda_cvec(
Df*np.concatenate(
(np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),
axis=0 ));
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full angular frequency vector
Nw=Nf; # number of non-negative angular f's
fpos = eda_cvec(Df * np.linspace(0,No2,Nf)); # non-neg f's
wpos=2*pi*fpos; # non-negative angular frequencies
# Compute fourier coefficients. Be careful not to forget
# the axis=0 (fft over elements of column). The default
# axis=1 (fft over elements of rows) won't work!
mest = eda_cvec( np.fft.fft(d, axis=0) );
# ampltude spectrum (note is un-normalized)
mestpos = eda_cvec( mest[0:Nf,0] ); # pos f's only
s = np.abs(mestpos); # a.s.d.
# plot spectrum
fig2 = plt.figure(2,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
# linear in frequency plot
ax1 = plt.subplot(1,3,1);
plt.axis([0, fmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
plt.plot(fpos[1:Nf,0],s[1:Nf,0],'k-');
plt.xlabel('frequency f (cycles per day)');
plt.ylabel('amplitude spectrum');
plt.title('linear in frequency');
# linear in perior plot
ax1 = plt.subplot(1,3,2);
Tmax = 1000;
plt.axis([0, Tmax, 0, 1.1*np.max(s[1:Nf,0]) ]);
plt.plot(np.reciprocal(fpos[1:Nf,0]),s[1:Nf,0],'k-');
plt.xlabel('periot T (days)');
plt.ylabel('amplitude spectrum');
plt.title('linear in period');
# semilog plot
ax1 = plt.subplot(1,3,3);
plt.semilogx(fpos[1:Nf,0], s[1:Nf,0], 'k-');
plt.axis([0.001, 1, 0, 1.1*np.max(s[1:Nf,0])])
plt.xlabel('frequency f (cycles per day)');
plt.ylabel('amplitude spectrum');
plt.title('logarithmic in frequency');
plt.show();
# compute Error of d->mest->dnew
# Inverse FFT: Don't forget the axis=0!
dnew = eda_cvec( np.real(np.fft.ifft(mest,axis=0)));
e = d - dnew; # individual errors
E = np.matmul( e.T, e ); E=E[0,0]; # total error
print('Error of d->mest->dnew:', E );
Error of d->mest->dnew: 1.9267235886294418e-20
In [8]:
# edapy_06_07
# spectra of Normal curves of varying widths
# time series of length N=256 with sampling interval of Dt=1
# generic time set up
N=256; # must be an even number
Dt=1.0;
tmin = 0;
tmax = Dt*(N-1);
t = eda_cvec( np.linspace(tmin,tmax,N) );
# generic frequency set up
No2 = floor(N/2); # N divided by two as an integer
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency increment
Nf=No2+1; # number of non-negative frequencies
# full vector of frequencies
f=eda_cvec(Df*np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full vector of angular fequencies
Nw=Nf; # number of angular frequencies
fpos=eda_cvec(Df*np.linspace(0,No2,Nf)); # non-negative frequencies
wpos=2*pi*fpos; # non-negative angular frequencies
# stamdard deviations, chosen to make an aestheic plot
Nsd=6;
sd=eda_cvec( [1.25, 2.5, 5, 10, 20, 40] );
# matrix D with each column a time series, one column for each sd
D = np.zeros((N,Nsd));
# matrix Ds with each column an amplitude spectral density
Ds = np.zeros((Nf,Nsd));
# populate D and Ds
for j in range(Nsd):
t1 = tmax/2; # mean
sd2 = sd[j]**2; # variance
d = np.exp( -np.power(t-t1,2) / (2*sd2) ); # normal function
D[:,j] = d.ravel();
# compute fourier coefficients
dfft = eda_cvec( np.fft.fft(d, axis=0) );
# delete negative frequencies
dfftpos = eda_cvec( dfft[0:Nf,0] );
# ampltude spectrum
Ds[:,j]=np.abs(dfftpos).ravel();
eda_draw(D[:,0],D[:,1],D[:,2],D[:,3],D[:,4],D[:,5],' ',Ds[:,0],Ds[:,1],Ds[:,2],Ds[:,3],Ds[:,4],Ds[:,5]);
In [9]:
# edapy_06_08 spike function
# generic time set up
N=256; # length of time sereis; must be even
Dt=1.0; # samppling interval
tmin = 0;
tmax = Dt*(N-1);
t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector
# spike time series of length N=256 with sampling interval of Dt=1
d = np.zeros((N,1));
d[0,0]=1.0;
# generic frequency setup
No2 = floor(N/2);
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# full (pos and neg) frequency vector
f = eda_cvec( Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full (pos and neg) angular frequencies
Nw=Nf; # number rof non-negative angular frequencies
fpos = eda_cvec(Df * np.linspace(0,No2,Nf) ); # non-negative frequencies
wpos= 2*pi*fpos; # non-negathive angular frequencies
# take fourier transform
# Note: Don't forget axis=0 option!
dt = eda_cvec( np.fft.fft(d, axis=0) );
# plot
fig1 = plt.figure(1,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(2,1,1);
plt.axis( [-tmax/10, tmax, -0.5, 1.5 ]);
plt.plot(t,d,'k-');
plt.plot(t,d,'k.');
plt.xlabel('time t');
plt.ylabel('d(t)');
plt.title('time series');
ax1 = plt.subplot(2,1,2);
plt.axis( [-fmax, fmax, 0, 2] );
plt.plot(f,np.real(dt),'k-'); # the transform should be purely real, so take real part
plt.plot(f,np.real(dt),'k.');
plt.xlabel('frequency f');
plt.ylabel('d-tilde(f)');
plt.title('time series');
plt.show();
In [11]:
# edapy_06_09: area under curve using FFT
# time series of length N=256 with sampling interval of Dt=1
# generic timeset up
N=256; # length of time series, must be even
Dt=1.0; # sampling interval
tmin = 0; # minimum time
tmax = Dt*(N-1); # maximum time
t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector
# generic frequency up
No2 = floor(N/2); # integer version of N/2
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency smaping
Nf=No2+1; # number on non-negative frequencies
f=eda_cvec(Df*np.concatenate(( np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));
Dw=2*pi*Df; # angular frequeny increment
w=2*pi*f; # full (pos and neg) ana
Nw=Nf; # number of non-negative angular frequencies
fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies
wpos=2*pi*fpos; # non-negative angular frequencies
# exemplary function
w0 = 2*pi*fmax/20; # frequency of oscilation
d= np.multiply(np.cos(w0*t),np.exp(-0.25*w0*t)); # product of cos(0 and delining exp()
# Fourier transform; note normalization factor of Dt
# and use of axis=0 option.
dt = eda_cvec( Dt*np.fft.fft(d, axis=0) );
# area by fft
area2 = np.real(dt[0,0]);
# area by Reimann summation
area1 = Dt*np.sum(d);
# plot
fig1 = plt.figure(1,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,1,1);
plt.axis( [0, tmax, -1, 1] );
plt.plot(t,d,'k-');
plt.xlabel('time t');
plt.ylabel('d(t)');
titlestr = "sum area %f fft area %f" % (area1, area2);
plt.title(titlestr);
plt.show();
In [12]:
# edapy_06_10: time shift using FFT
# generic time= set up
N=256;
Dt=1.0;
tmin = 0.0;
tmax = Dt*(N-1);
t=eda_cvec(np.linspace(tmin,tmax,N));
# generic time/frequency set up
No2 = floor(N/2);
fmax=1/(2.0*Dt);
Df=fmax/No2;
Nf=No2+1;
f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df;
w=2*pi*f;
Nw=Nf;
fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
wpos=2*pi*fpos;
# exemplary function
d = np.zeros((N,1));
w0 = 2*pi*fmax/20;
d[:,0] = np.multiply(np.cos(w0*t),np.exp(-0.25*w0*t)).ravel();
# time shift
t0 = t[15,0];
dt = (np.multiply(
np.exp((complex(0,-1)*w*t0)), np.fft.fft(d, axis=0)));
# inverse Fourier transform
# result should be purely real, so take real part
# to get rid of zero imaginary part
ds = eda_cvec( np.real((np.fft.ifft(dt,axis=0))));
# plot
fig1 = plt.figure(1,figsize=(10,7));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(2,1,1);
plt.axis( [0, tmax, -1, 1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
plt.title('time series before shift');
ax1 = plt.subplot(2,1,2);
plt.axis( [0, tmax, -1, 1] );
plt.plot(t,ds,'k-');
plt.xlabel('time t');
plt.ylabel('d(t)');
titlestr = "time series after rshift of %.2f" % (t0);
plt.title(titlestr);
plt.show();
In [13]:
# eda06_11: derivative using FFT
# time series of length N=256 with sampling interval of Dt=1
N=256;
Dt=1.0;
tmin = 0;
tmax = Dt*(N-1);
t=eda_cvec(np.linspace(tmin,tmax,N));
# generic time/frequency set up
No2 = floor(N/2);
fmax=1/(2.0*Dt);
Df=fmax/No2;
Nf=No2+1;
f=eda_cvec(Df * np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));
Dw=2*pi*Df;
w=2*pi*f;
Nw=Nf;
fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
wpos=2*pi*fpos;
# exemplary function
d = np.zeros((N,1));
t0 = tmax/2;
sd = 50;
d[:,0] = np.exp(-np.power((t-t0),2)/(2*sd**2)).ravel();
# take Fourier transform and multiply by iw
# except for zero freuency
dt = np.multiply(complex(0,1)*w , np.fft.fft(d, axis=0));
# inverse Fourier transform
dddt = np.real((np.fft.ifft(dt,axis=0)));
# derivative by finite differences
dddt2 = np.zeros((N,1));
dddt2[0:N-1,0]=(d[1:N,0]-d[0:N-1,0])/Dt;
# plot
fig1 = plt.figure(1,figsize=(10,8));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(3,1,1);
plt.axis( [0, tmax, -1.1*np.max(d), 1.1*np.max(d)] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
ax1 = plt.subplot(3,1,2);
plt.axis( [0, tmax, -1.1*np.max(dddt), 1.1*np.max(dddt)] );
plt.plot(t,dddt,'k-');
plt.xlabel('time t');
plt.ylabel('dd/dt by fft');
ax1 = plt.subplot(3,1,3);
plt.axis( [0, tmax, -1.1*np.max(dddt2), 1.1*np.max(dddt2)] );
plt.plot(t,dddt2,'k-');
plt.xlabel('time t');
plt.ylabel('dd/dt by sum');
plt.show();
In [14]:
# edapy_06_12: integral using FFT
# generic time/y set up
N=256;
Dt=1.0;
tmin = 0;
tmax = Dt*(N-1);
t=eda_cvec(np.linspace(tmin,tmax,N));
# generic frequency set up
No2 = floor(N/2);
fmax=1/(2.0*Dt);
Df=fmax/No2;
Nf=No2+1;
f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df;
w=2*pi*f;
Nw=Nf;
fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
wpos=2*pi*fpos;
# exemplary function
sd = 50;
t1 = tmax/4;
t2 = 3*tmax/4;
sd = 25;
d = (np.exp(-np.power((t-t1),2)/(2*sd**2)) - np.exp(-np.power((t-t2),2)/(2*sd**2)) );
# take Fourier transform and divide by iw
# except set zero frequency to zero
wr=np.concatenate(
(np.zeros((1,1)),np.reciprocal(complex(0,1)*w[1:N,0:1])) );
dt = np.multiply(wr, np.fft.fft(d, axis=0));
# inverse Fourier transform, don't forget the axis=0
# since result should be purely real, delete imaginary part
integral1 = eda_cvec( np.real((np.fft.ifft(dt,axis=0))));
# integeral by Reimann sum
integral2 = eda_cvec((Dt*np.cumsum(d)));
# adjust so that integrals have the same integration constant (zero)
integral1 = integral1 - integral1[0,0];
integral2 = integral2 - integral2[0,0];
# plot
fig1 = plt.figure(1,figsize=(10,8));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(3,1,1);
plt.axis( [0, tmax, -1.1*np.max(d), 1.1*np.max(d)] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
ax1 = plt.subplot(3,1,2);
plt.axis( [0, tmax, -1.1*np.max(integral1), 1.1*np.max(integral1)] );
plt.plot(t,integral1,'k-');
plt.xlabel('time t');
plt.ylabel('integral by fft');
ax1 = plt.subplot(3,1,3);
plt.axis( [0, tmax, -1.1*np.max(integral2), 1.1*np.max(integral2)] );
plt.plot(t,integral2,'k-');
plt.xlabel('time t');
plt.ylabel('dd/dt by diff');
plt.show();
In [15]:
# edapy_06_13: Crib Sheet 6.2
# makes exemplary time series d(t),
# amplitude spectral density s(f) pairs for
# time series of length N=256 with sampling interval of Dt=1
N=128;
Dt=0.5;
tmin = 0;
tmax = Dt*(N-1);
t=eda_cvec(np.linspace(tmin,tmax,N));
# generic frequency set up
No2 = floor(N/2);
fmax=1/(2.0*Dt);
Df=fmax/No2;
Nf=No2+1;
f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df;
w=2*pi*f;
Nw=Nf;
fpos=eda_cvec(Df * np.linspace(0,No2,Nf));
wpos=2*pi*fpos;
# duration of time series
T = tmax-tmin;
# plot 1
fig1 = plt.figure(1,figsize=(10,8));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
p=1;
# narrow function
sigmad = 1.0;
d = (np.exp(-np.power(t-tmax/2,2)/(2*(sigmad**2)))+np.random.normal(0,0.1,(N,1)));
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform, note normalization of Dt
dfftpos = eda_cvec( dfft[0:Nf,0] ) ; # positive fequencies, only
ds=sqrt(2/T)*np.abs(dfftpos); # properly normalized amplitude spectral density
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.ylabel('s(f)');
p=p+1;
# zero area wavelet
sigmad1 = 2.5;
sigmad2 = 1.0;
d = np.multiply( np.exp(-np.power(t-T/2,2)/(2*(sigmad1**2))) , np.cos((t-T/2)/sigmad2));
d = d + np.random.normal(0,0.1,(N,1));
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );
dfftpos = eda_cvec( dfft[0:Nf,0] );
ds=sqrt(2/T)*np.abs(dfftpos);
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.ylabel('s(f)');
p=p+1;
# wide function
sigmad = 7.0;
d = np.exp(-np.power(t-tmax/2,2)/(2*(sigmad**2))) + np.random.normal(0,0.1,(N,1));
# fourier transform
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # note Dt normalization
dfftpos = eda_cvec( dfft[0:Nf,0] ); # keep non-negative frequencies
ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral density
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.ylabel('s(f)');
p=p+1;
# triangle
No4 = floor(N/4);
d = np.zeros((N,1));
d[No4:2*No4,0] = np.linspace(0,1,No4);
d[2*No4:3*No4,0] = np.linspace(1,0,No4);
d[:,0] = (d + np.random.normal(0,0.01,(N,1))).ravel();
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normalization
dfftpos = eda_cvec( dfft[0:Nf,0] ); # keep positive frequencies
ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral density
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.xlabel('time t');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.xlabel('frequency f');
plt.ylabel('s(f)');
p=p+1;
plt.show();
print("\n\n\n\n");
# plot 2
fig2 = plt.figure(2,figsize=(10,8));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
p=1;
# approximately sinusoidal, with time varying frequency w0
w0 = eda_cvec( (0.25*2*pi*(0.95+0.1*np.linspace(0,N-1,N)/(N-1))) );
d = np.sin(np.multiply(w0,t)) + np.random.normal(0,0.2,(N,1));
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );
dfftpos = eda_cvec( dfft[0:Nf,0] );
ds = sqrt(2/T)*np.abs(dfftpos);
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.ylabel('s(f)');
p=p+1;
# beating, two sinusoids with frequencies in 2:1 ratio
w1 = 0.25*2*pi;
w2 = w1/2;
d = np.sin(w1*t) + 0.5*np.sin(w2*t) + np.random.normal(0,0.2,(N,1));
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier trasnform with Dt normailzation
dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative fequencies only
ds =sqrt(2/T)*np.abs(dfftpos).ravel(); # amplitude spectral density
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.ylabel('s(f)');
p=p+1;
# chirp, frequency that increases with time
w0 = eda_cvec(0.25*2*pi*(2*np.linspace(0,N-1,N)/(N-1)));
d = np.sin(np.multiply(w0,t)) + np.random.normal(0,0.05,(N,1));
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normaization
dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative frequencies only
ds = sqrt(2/T)*np.abs(dfftpos);
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.ylabel('s(f)');
p=p+1;
# sawtooth
d = ((np.mod(t,N/10)/(N/10)-0.5)+np.random.normal(0.0,0.005,(N,1)));
dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normalization
dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative frequencies
ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral ratio
ax1 = plt.subplot(4,2,p);
plt.axis( [0, T, -1.1, 1.1] );
plt.plot(t,d,'k-');
plt.xlabel('time t');
plt.ylabel('d(t)');
p=p+1;
ax1 = plt.subplot(4,2,p);
plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );
plt.plot(fpos,ds,'k-');
plt.xlabel('frequency f');
plt.ylabel('s(f)');
p=p+1;
plt.show();
In [18]:
# edapy_06_14: power spectral density
# comparison of variances between time series and its power spectral density
# random time series with zero mean and variance sigmad^2
N=256; # length of time series, must be even
Dt=0.01; # sampling interval
tmin = 0;
tmax = Dt*(N-1);
t=eda_cvec(Dt*np.linspace(0,N-1,N));
sigmad1=10; # standard deviation
sigmad1sq=sigmad1**2; # variance
d = np.random.normal(0,sigmad1,(N,1)); # random timeseroes
# this code will also work for a timeseries with non-zero mean,
# but then the "variance" is really "power"
d = d-np.mean(d); # this statement can be commented out for "power"
# variance estimated from time series
dsq = np.matmul(d.T,d); dsq=dsq[0,0];
sigmad2sq = dsq/N;
# standard frequency setup for cos() and sin() representation of data
Nf=floor(N/2+1); # number o fnon-negative frequencies
fmax = 1/(2*Dt); # nyquist frequency
Df = fmax/(N/2); # frequency interval
f = eda_cvec( np.linspace(0,fmax,Nf) ); # positive frequencies
Nw=Nf; # number of non-negative angular frequencies
wmax = 2*pi*fmax; # maximum angular frequency
Dw = 2*pi*Df; # angular frequency interval
w = 2*pi*f; # non-negative angular frequencies
# set up data kernel G for fourier coefficients
M=N;
G = np.zeros((N,M));
# zero frequency column
G[:,0] = np.ones((N,1)).ravel();
# interior M/2-1 columns
# note index arithmetic due to columns occuring in pairs
Mo2 = floor(M/2);
for i in range(1,Mo2):
j = 2*i-1;
k = j+1;
G[:,j]=np.cos(w[i,0]*t).ravel();
G[:,k]=np.sin(w[i,0]*t).ravel();
# nyquist column
G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();
# solve least squares problem using
# analytic formula for inv(GT*G)
GTd = np.matmul( G.T, d );
# since inv( GT G ) is known analytically, use shortcut
gtgi = 2*np.ones((M,1))/N;
gtgi[0,0]=1/N;
gtgi[M-1,0]=1/N;
mest = np.multiply( gtgi , GTd );
e = d - np.matmul(G,mest); # individual errors
E = np.matmul( e.T , e ); E=E[0,0]; # overall erros
# print('error', E); # uncomment to check that error is zero
# generic frequency setup for fourier transform
No2 = floor(N/2); # integer version of N/2
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# full (pos, neg) frequencies
f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df; # angular freuqncy samppling
w=2*pi*f; # angular frequencies
Nw=Nf; # number of non-negative angular frequencies
fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies
wpos=2*pi*fpos; # non-negative angular frequencies
# discrete fourier transform, nno normalization
C = eda_cvec( np.fft.fft(d, axis=0) );
# inverse Fourier transform; result should be purely real, so take real part
drecon = eda_cvec( np.real((np.fft.ifft(C,axis=0))) );
e = d - drecon; # individual errors
E = np.matmul( e.T, e ); E=E[0,0]; # total error
# print('error', E); # uncomment to check that error is zero
# comparison of coefficients
jf=10;
kf=2*jf-1;
printstr = "check agreement between sines/cosine and fft at frequency %d" % (jf);
print(printstr);
printstr = "A2 B2 %.4f %.4f" % (mest[kf,0], mest[kf+1,0]);
print(printstr);
printstr = "(2/N)real(C) -(2/N)imag(C) %.4f %.4f" % ((2/N)*np.real(C[jf,0]), -(2/N)*np.imag(C[jf,0]));
print(printstr);
printstr = "variance estimated from time series %.4f" % (sigmad2sq);
print(printstr);
sigmad3sq = 0.5*np.matmul(mest.T,mest); sigmad3sq=sigmad3sq[0,0];
# normally one doesn't bother to correct the first and last frequency, but ...
sigmad3sq = sigmad3sq + 0.5*(mest[0,0]**2 + mest[N-1,0]**2)
printstr = "variance estimated from sine/cosines %.4f" % (sigmad3sq);
print(printstr);
# variance estimated from fft coefficients
C2 = np.zeros((Nf,1))
C2[:,0] = np.real(np.multiply(C[0:Nf,0].conj(),C[0:Nf,0])).ravel();
sigmad4sq = (2/(N**2)) * np.sum(C2);
# normally one doesn't bother to correct the first and last frequency, but ...
sigmad4sq = sigmad4sq - 0.5*(2/(N**2))*C2[0,0] - 0.5*(2/(N**2))*C2[Nf-1,0];
printstr = "variance estimated from fft coefficients %.4f" % (sigmad4sq);
print(printstr);
# approximation to integral transform
dbar=Dt*C;
# variance estimated from fourier transform
dbar2 = eda_cvec( np.real(np.multiply(dbar[0:Nf,0].conj(),dbar[0:Nf,0])));
sigmad4sq = (2/((Dt**2)*(N**2))) * np.sum( dbar2 );
# normally one doesn't bother to correct the first and last frequency, but ...
sigmad4sq = sigmad4sq - 0.5*(2/((Dt**2)*(N**2)))*dbar2[0,0] - 0.5*(2/((Dt**2)*(N**2)))*dbar2[Nf-1,0];
printstr = "variance estimated from fourier transform %.4f" % (sigmad4sq);
print(printstr);
# integral from 0 to fny of |dbar(f)|2
theintegral = Df * np.sum( dbar2 );
sigmad5sq = (2/(Df*(Dt**2)*(N**2))) * theintegral;
# normally one doesn't bother to correct the first and last frequency, but ...
sigmad5sq = sigmad5sq - 0.5*Df*(2/(Df*(Dt**2)*(N**2)))*dbar2[0,0] - 0.5*Df*(2/(Df*(Dt**2)*(N**2)))*dbar2[Nf-1,0];
printstr = "variance estimated from frequency integral %.4f" % (sigmad5sq);
print(printstr);
# definition of power spectral density
T = tmax-tmin+Dt; # note the correction of the length of the window by Dt
thepsd = (2/T) * dbar2;
# integral of power spectral density
sigmad6sq = Df * np.sum(thepsd);
# normally one doesn't bother to correct the first and last frequency, but ...
sigmad6sq = sigmad6sq - 0.5*Df*thepsd[0,0] - 0.5*Df*thepsd[Nf-1,0];
printstr = "variance estimated from power spectral density %.4f" % (sigmad6sq);
print(printstr);
check agreement between sines/cosine and fft at frequency 10 A2 B2 -0.0177 0.7437 (2/N)real(C) -(2/N)imag(C) -0.0177 0.7437 variance estimated from time series 110.1368 variance estimated from sine/cosines 110.1368 variance estimated from fft coefficients 110.1368 variance estimated from fourier transform 110.1368 variance estimated from frequency integral 110.1368 variance estimated from power spectral density 110.1368
In [19]:
# edapy_06_15: p.s.d. of the Neuse River hydrograph
# properly normalized amplitude spectral density of
# Neuse River hydrograph using FFT
# load data from file
D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
[Nraw, K]=D.shape;
traw = eda_cvec( D[:,0] ); # time
Dt = traw[1,0]-traw[0,0]; # sampling interval
draw = eda_cvec( D[:,1] ); # data
# round off to even number of points
No2 = floor(Nraw/2);
N=2*No2;
d = eda_cvec( draw[0:N,0] );
t = eda_cvec( traw[0:N,0] );
# generic time set up
tmin = 0;
tmax = Dt*(N-1);
t=eda_cvec(np.linspace(tmin,tmax,N));
# generic frequency set up
No2 = floor(N/2); # integer version of N/2
M=N; # number of model parameters equal number of data
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequenices
# full (pos, neg) frequency vector
f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full (pos, neg) anguar frequency vector
Nw=Nf; # number of non-negative angular frequencies
fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies
wpos=2*pi*fpos; # non-negative angular frequencies
# compute fourier coefficients. Note normalization of Dt
dbar = Dt*np.fft.fft(d, axis=0); # fourier transform
dbarpos = dbar[0:Nf,0:1]; # non-negative frequencies
T = N*Dt; # time series length
# p.s.d.
s2 = (2/T) * np.real(np.multiply(dbarpos.conj(), dbarpos));
# plot data
fig1 = plt.figure(1,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);
plt.plot(t,d,'k-');
plt.xlabel('time (days)');
plt.ylabel('discharge (cfs)');
plt.title('Neuse River Hydrograph');
plt.show();
# plot power spectral density
fig2 = plt.figure(2,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,3,1);
plt.axis([0, fmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);
plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');
plt.xlabel('frequency f (cycles per day)');
plt.ylabel('power spectral density');
plt.title('linear in frequency');
ax1 = plt.subplot(1,3,2);
Tmax = 1000;
plt.axis([0, Tmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);
plt.plot(np.reciprocal(fpos[1:Nf,0]),s2[1:Nf,0],'k-');
plt.xlabel('periot T (days)');
plt.ylabel('power spectral density');
plt.title('linear in period');
ax1 = plt.subplot(1,3,3);
plt.axis([0.001, 1, 0, 1.1*np.max(s2[1:Nf,0])]);
plt.semilogx(fpos[1:Nf,0], s2[1:Nf,0], 'k-');
plt.xlabel('frequency f (cycles per day)');
plt.ylabel('power spectral density');
plt.title('logarithmic in frequency');
plt.show();
# variance calculation
v1 = np.std(d)**2;
printstr = "variance estimated from time series %.6e" % (v1);
print(printstr);
v2 = Df * np.sum(s2[1:Nf,0]);
printstr = "variance estimated from power spectral density %.6e" % (v2);
print(printstr);
variance estimated from time series 8.396180e+06 variance estimated from power spectral density 8.396180e+06
In [22]:
# edapy_06_16
# make a random noise dataset
# Note that he noise is saved to mynoise.txt, not noise.txt
# to avoid overwriting the standard distribution. So change
# the name if you want to overwrite noise.txt
# generic time setup
N=1024; # length of time series, must be even
Dt=1.0; # sampling interval
tmin = 0.0; # start time
tmax = Dt*(N-1); # end time
t = eda_cvec(np.linspace(tmin,tmax,N)); # time vector
# make vector d of random noise
dbar=0.0; # mean
sigmad=1.0; # standard deviation
d = np.random.normal( dbar, sigmad, (N,1) );
# plot data
fig1 = plt.figure(1,figsize=(10,4));
plt.rcParams['xtick.bottom'] = True;
plt.rcParams['xtick.labelbottom'] = True;
plt.rcParams['xtick.top'] = False;
plt.rcParams['xtick.labeltop'] = False;
ax1 = plt.subplot(1,1,1);
r = dbar+3*sigmad; # this range encloses 99% of data
plt.axis([tmin, tmax, -r, r ]);
plt.plot(t,d,'k-');
plt.plot(t,d,'k.');
plt.plot([tmin,tmax],[0,0],'k-');
plt.xlabel('time');
plt.ylabel('d(t)');
titlestr = "Normally-distributed random timeseries witm mean %.2f and sigma %.2f" % (dbar, sigmad);
plt.title(titlestr);
plt.show();
# save matrix to file of tab-separated values
# the file name has been changed from "noise.txt" to "noise2.txt"
# so as not to overwrite the existing file in the distributuon
Dm = np.concatenate( (t,d),axis=1);
np.savetxt("../Data/mynoise.txt",Dm, delimiter="\t");
In [ ]: