{{{id=15|
# Copyright 2009, P. Lutus
# Released under the GPL (http://www.gnu.org/copyleft/gpl.html)
///
}}}
{{{id=4|
#auto
# reset() # commented out for now -- ticket 7255
forget()
# special equation rendering
def render(x,name = "temp.png",size = "normal"):
if(type(x) != type("")): x = latex(x)
latex.eval("\\" + size + " $" + x + "$",{},"",name)
# magnitude of 2-component cartesian vector
def mag(x):
return sqrt(x[0]^2+x[1]^2)
var('a, b, c, d, f, g, m, q, t, x, y')
# time-domain function plot
def time_plot_funct(funct,nl,w=0.25,color='#006000',labels=('Time','Amplitude')):
n = 500
list = [[t*w/n,funct(t*w/n,nl)] for t in range(n)]
return list_plot(list,rgbcolor=color,axes_labels=labels,plotjoined=True)
# time-domain FFT array plot
def time_plot_fft(fftobj,freq,line_color='#006000',labels=('Time','Amplitude')):
lfft = len(fftobj)
data = [[x*2*freq/lfft,fftobj[x][0]] for x in range(lfft/2)]
return list_plot(data,plotjoined=True,rgbcolor=line_color,axes_labels=labels)
# frequency domain plot
def fft_plot(fftobj,line_color='#006000',labels=('Frequency','Amplitude')):
size = len(fftobj)
dt = 1.0/size
list = map(lambda x:(mag(x))*dt,fftobj[:size/2])
return list_plot(list,rgbcolor=line_color,plotjoined=True,axes_labels=labels)
///
}}}
{{{id=68|
def tf(t,nl): return sin(100*pi*t) + sin(200*pi*t) + sin(300*pi*t) + (random()-0.5) * nl
def time_freq(a,nl):
tp = time_plot_funct(tf,nl,0.05,'#006000',('',''))
tlbl=text("Time, Seconds",(0.025,-1),rgbcolor='black')
show(tp+tlbl,figsize=(3.5,2),fontsize=10,ymin=-2.5,ymax=2.5)
n = 400
fft = FFT(n)
for t in range(n):
fft[t] = tf(t/n,nl)
fft.forward_transform()
fp = fft_plot(fft,'#006000',('',''))
flbl=text("Frequency, Hertz",(100,-0.15),rgbcolor='black')
show(fp+flbl,figsize=(3.5,2),fontsize=10,ymin=-.2,ymax=0.6)
time_freq(t,0)
time_freq(t,6)
///
}}}
{{{id=75|
render("$f(y) = \int_{-\infty}^\infty f(x) \, e^{-2\pi i x y} dx$","equ_ftransform.png","large")
///
}}}
{{{id=76|
render("$f(x) = \int_{-\infty}^\infty f(y) \, e^{2\pi i x y} dy$","equ_finvtransform.png","large")
///
}}}
{{{id=77|
render("$e^{2\pi i \\theta} = cos \, 2\pi\\theta + i \, sin \, 2\pi\\theta $","equ_eulers_equation.png","large")
///
}}}
{{{id=78|
lbl = text("$e^{2\pi i \\theta} = cos \, 2\pi\\theta + i \, sin \, 2\pi\\theta $"
,(1,1.2),fontsize=14,rgbcolor='black')
rp = plot(lambda x:N(e^(2*pi*i*x)).real(),x,0,2,rgbcolor='#006000')
ip = plot(lambda x:N(e^(2*pi*i*x)).imag(),x,0,2)
show(rp+ip+lbl,figsize=(4,3))
///
}}}
{{{id=84|
render("$X_k = \sum_{n=0}^{N-1} x_n e^{-\\frac{2\pi i}{N} k n} \, k=0,...,N-1$","equ_dft.png","large")
///
}}}
{{{id=86|
render("$x_n = \\frac{1}{N} \sum_{k=0}^{N-1} X_k e^{\\frac{2\pi i}{N} k n} \, n=0,...,N-1$","equ_idft.png","large")
///
}}}
{{{id=87|
ff(f,t) = sin(2*pi*f*t)
f = 5
@interact
def _(n = (5..20)):
list = [ff(f,t/(2*n)) for t in range(0,n)]
show(list_plot(list,plotjoined=True),figsize=(4,3))
///
}}}
{{{id=88|
# declare constants
epsilon = 1e-12
samples = 25
# define a waveform-generating function (wgf)
wgf(t) = 12*sin(t*3) + 11*sin(t*4) + 10*sin(t*5) + 9*sin(t*6)
# display the time-domain waveform for the wgf
show(time_plot_funct(wgf,0,10),figsize=(4,3))
# create a Fast Fourier Transform object
# and populate it with wgf data
fft = FFT(samples)
for t in range(samples):
fft[t] = wgf(2*pi*t/samples)
# convert to frequency domain
fft.forward_transform()
# locate and print spectral lines
for f in range(samples/2):
v = 2*mag(fft[f])/samples
if(v > epsilon):
print "Frequency: %3.0f, Amplitude: %3.0f" % (f,v)
///
Frequency: 3, Amplitude: 12
Frequency: 4, Amplitude: 11
Frequency: 5, Amplitude: 10
Frequency: 6, Amplitude: 9
}}}
{{{id=90|
# declare constants
epsilon = .05
samples = 1000
# define a waveform-generating function (wgf)
def wgf(t,nl): return sgn(N(cos(2*pi*t)))
# display the time-domain waveform for the wgf
show(time_plot_funct(wgf,0,2),figsize=(4,3))
# create a Fast Fourier Transform object
# and populate it with wgf data
fft = FFT(samples)
for t in range(samples):
fft[t] = wgf(t/samples,0)
# convert to frequency domain
fft.forward_transform()
# locate and print spectral linessamples = 10000
for f in range(samples/2):
v = mag(fft[f])/(samples*2/pi)
if(v > epsilon):
# print "Frequency: %3.0f, Amplitude: 1/%.0f" % (f,1/v)
print "Frequency: %3.0f, Amplitude: %.4f" % (f,v)
///
Frequency: 1, Amplitude: 1.0000
Frequency: 3, Amplitude: 0.3333
Frequency: 5, Amplitude: 0.2000
Frequency: 7, Amplitude: 0.1428
Frequency: 9, Amplitude: 0.1111
Frequency: 11, Amplitude: 0.0909
Frequency: 13, Amplitude: 0.0769
Frequency: 15, Amplitude: 0.0666
Frequency: 17, Amplitude: 0.0588
Frequency: 19, Amplitude: 0.0526
}}}
{{{id=98|
sgn(N(cos(3)))
///
-1
}}}
{{{id=93|
# declare constants
samples = 5000
@interact
def _(frequency = (1..8),harmonics = (1..64)):
# create a Fast Fourier Transform object
fft = FFT(samples)
# populate the FFT object with spectral data
fd2 = frequency*2.0
for n in range(harmonics):
a = fd2*(1+n*2)
b = -fd2*4/(a*pi)
fft[a] = (0,b)
# convert from frequency domain to time domain
fft.backward_transform()
# display the result
show(time_plot_fft(fft,1),ymin=-1.3,ymax=1.3)
///
}}}
{{{id=95|
# just to produce an article graphic
# declare constants
samples = 5000
frequency = 2
harmonics = 64
# create a Fast Fourier Transform object
fft = FFT(samples)
# populate the FFT object with spectral data
fd2 = frequency*2.0
for n in range(harmonics):
a = fd2*(1.0+n*2.0)
b = -fd2*4.0/(a*pi)
fft[a] = (0,b)
# comvert from frequency domain to time domain
fft.backward_transform()
# display the result
show(time_plot_fft(fft,1),ymin=-1.3,ymax=1.3,figsize=(4,3))
///
}}}
{{{id=96|
render("$y = \\frac{4}{\pi} \sum_{n=1}^\infty \\frac{sin(2\pi f t(2n-1))}{2n-1} $","equ_sqw.png","large")
///
}}}
{{{id=97|
///
}}}