-
Notifications
You must be signed in to change notification settings - Fork 0
/
int_fourier.py
132 lines (108 loc) · 3.69 KB
/
int_fourier.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
from numpy import *
from scipy import *
from pylab import *
import numpy.random as random
import pdb
ion()
i = 1j
def my_slow_fft(f):
# a slow ifft that **CAN'T** interpolate!
N = len(f)
F = zeros(N, dtype=complex)
for k in range(N):
for n in range(N):
F[k] += f[n]*e**(-2.0*pi*i*k*n/N)
return F
def my_slow_ifft(F):
# a slow ifft that **CAN'T** interpolate!
N = len(F)
f = zeros(N, dtype=complex)
for n in range(N):
for k in range(N):
f[n] += F[k]*e**(2.0*pi*i*k*n/N)
return (1.0/N)*f
def my_slow_ifft_interp(F,t):
# a slow ifft that **CAN'T** interpolate!
N = len(F)
f = zeros(len(t), dtype=complex)
for n in range(len(t)):
for k in range(N):
f[n] += F[k]*e**(2.0*pi*i*k*t[n])
return (1.0/N)*f
def my_correct_slow_ifft_interp(F,t):
# an ifft that **CAN** interpolate
N = len(F)
f = zeros(len(t), dtype=complex)
for n in range(len(t)):
for k in range(int(N/2)):
f[n] += F[k]*e**(2.0*pi*i*k*t[n])
for k in range(int(N/2),N):
f[n] += F[k]*e**(-2.0*pi*i*(N-k)*t[n])
return (1.0/N)*f
if __name__ == "__main__":
# Sample measurements
ts = linspace(0,1,num=21)
f = zeros(len(ts))
hf = int(len(ts)/2)
f[0:hf] = ts[0:hf]
#f[hf:] = e**-ts[hf:] + 0.5 - e**-0.5
f[hf:] = 1.0/ts[hf:] - 1.5
#f = ts*(ts-1)
#f = sin(ts*(2.*pi)) + cos(ts*(2.*pi))#-1.0*ts**2#exp(ts)
F = fftshift(fft(f))
F_unshifted = fft(f)
ff = fftshift(ifftshift(f))
# plot where we are so far
# this includes F^-1{F(f)} for now
figure(1)
plot(ts,f, 'bo-',markerfacecolor='none', mec='blue')
#ylim(-0.6,0.6)
title("Original Function")
figure(2)
plot(ts,F.real,ts,F.imag)
legend(('real', 'imaginary'))
title("Shifted Fourier Transform")
figure(22)
plot(ts,F_unshifted.real,ts,F_unshifted.imag)
legend(('real', 'imaginary'))
title("Unshifted Fourier Transform")
figure(3)
plot(ts,ff)
#ylim(-0.6,0.6)
title("Getting Back The Original")
# now interpolate
# new timestamps, m samples between former samples.
Nf = len(ts)
t = linspace(0,1,num=len(ts)*3)
Fs = my_slow_fft(f)
fm = my_slow_ifft(Fs)
# using the Fourier Coefficients "sas-is"
fwi = my_slow_ifft_interp(F_unshifted,t)
fci = my_correct_slow_ifft_interp(F_unshifted,t)
figure(4)
plot(ts,Fs.real,ts,Fs.imag)
legend(('real', 'imaginary'))
title("My Unshifted Fourier Transform")
figure(5)
plot(ts,f, 'bo-',markerfacecolor='none', mec='blue')
plot(ts,fm,'g^-',markerfacecolor='none', mec='green')
#ylim(-0.6,0.6)
title("Original Using My Inverse")
figure(500)
plot(ts,f, 'bo-',markerfacecolor='none', mec='blue')
plot(ts[:-1],my_correct_slow_ifft_interp(fft(f[:-1]),ts[:-1]),'g^-',markerfacecolor='none', mec='green')
#ylim(-0.6,0.6)
title("Original Using My Inverse Interp ")
figure(6)
plot(ts,f, 'bo-',markerfacecolor='none', mec='blue')
plot(t,fwi,'go-',markerfacecolor='none', mec='green')
#ylim(-0.6,0.6)
legend(('original', 'interpolatation'))
title("Interpolating with Fourier Coefficients, No Change")
figure(7)
plot(ts,f,'bo-',markerfacecolor='none', mec='blue')
plot(t[:-1],my_correct_slow_ifft_interp(fft(f[:-1]),t[:-1]),'g^-',markerfacecolor='none', mec='green')
#ylim(-0.6,0.6)
legend(('original', 'interpolatation'))
title("Interpolating with Fourier Coefficients Accounting for Aliasing")
show()