-
Notifications
You must be signed in to change notification settings - Fork 4
/
test_pade.py
66 lines (49 loc) · 1.89 KB
/
test_pade.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
import numpy as np
from pade import pade
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
def process_signal(t,signal):
# Scipy FFT
SIGMA = 5 # control amount of lineshape broadening (bigger = more narrow)
dt = t[1] - t[0]
damp = np.exp(-(dt*np.arange(len(signal)))/float(SIGMA))
fw_fft = fft(signal*damp)
freq_fft = fftfreq(t.size,d=dt)*2*np.pi # 2pi is implicit in Pade, so we add it here
# grab the positive frequency components
fw_fft = fw_fft[:t.size//2]
freq_fft = freq_fft[:t.size//2]
# Now do Pade appx of Fourier Transform; line broadening via SIGMA
# we will read in frequencies from SciPy FFT to evaluate analytic function over (read_freq)
fw_pade, freq_pade = pade(t,signal,sigma=SIGMA,read_freq=freq_fft)
# grab just imaginary components
fw_pade = np.imag(fw_pade)
fw_fft = np.imag(fw_fft)
return fw_pade, freq_pade, fw_fft, freq_fft
def test_pade():
w = 2.0
dt = 0.02
N = 5000
t = np.linspace(0,dt*N, N, endpoint=False)
signal = np.sin(w*t) + np.sin(2*w*t) + np.sin(4*w*t)
fw_pade, freq_pade, fw_fft, freq_fft = process_signal(t,signal)
# FFT and Pade agree to within some tolerance
assert np.allclose(freq_pade,freq_fft)
assert np.linalg.norm(fw_fft - fw_pade) < 1e-1 # "closeness" also depends on signal length, SIGMA, etc.
#print(np.linalg.norm(fw_fft - fw))
if __name__ == '__main__':
w = 2.0
dt = 0.02
N = 5000
t = np.linspace(0,dt*N, N, endpoint=False)
signal = np.sin(w*t) + np.sin(2*w*t) + np.sin(4*w*t)
fw_pade, freq_pade, fw_fft, freq_fft = process_signal(t,signal)
plt.plot(t,signal)
plt.savefig('signal.png')
plt.clf()
plt.cla()
plt.plot(freq_pade,fw_pade,label='pade')
plt.plot(freq_fft,fw_fft,label='FFT',ls='--')
plt.xlim([0,10])
plt.legend()
plt.savefig('fsignal.png')
plt.show()