/
test_interpolation.py
144 lines (119 loc) · 5.12 KB
/
test_interpolation.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
133
134
135
136
137
138
139
140
141
142
143
144
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
The interpolation test suite for ObsPy.
"""
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from future.builtins import * # NOQA
import os
import unittest
import numpy as np
import matplotlib.pyplot as plt
from obspy.core.util.testing import ImageComparison
from obspy.signal.interpolation import (lanczos_interpolation,
calculate_lanczos_kernel,
plot_lanczos_windows)
class InterpolationTestCase(unittest.TestCase):
"""
Interpolation test case
"""
def setUp(self):
# directory where the test files are located
self.path = os.path.join(os.path.dirname(__file__), 'data')
self.path_images = os.path.join(os.path.dirname(__file__), 'images')
def test_calculate_lanczos_kernel(self):
"""
Tests the kernels implemented in C against their numpy counterpart.
"""
x = np.linspace(-5, 5, 11)
values = calculate_lanczos_kernel(x, 5, "hanning")
np.testing.assert_allclose(
values["only_sinc"], np.sinc(x), atol=1E-9)
np.testing.assert_allclose(
values["only_taper"], np.hanning(len(x)), atol=1E-9)
np.testing.assert_allclose(
values["full_kernel"], np.sinc(x) * np.hanning(len(x)),
atol=1E-9)
values = calculate_lanczos_kernel(x, 5, "blackman")
np.testing.assert_allclose(
values["only_sinc"], np.sinc(x), atol=1E-9)
np.testing.assert_allclose(
values["only_taper"], np.blackman(len(x)), atol=1E-9)
np.testing.assert_allclose(
values["full_kernel"], np.sinc(x) * np.blackman(len(x)),
atol=1E-9)
values = calculate_lanczos_kernel(x, 5, "lanczos")
np.testing.assert_allclose(
values["only_sinc"], np.sinc(x), atol=1E-9)
np.testing.assert_allclose(
values["only_taper"], np.sinc(x / 5.0), atol=1E-9)
np.testing.assert_allclose(
values["full_kernel"], np.sinc(x) * np.sinc(x / 5.0),
atol=1E-9)
def test_lanczos_interpolation(self):
"""
Tests against the instaseis implementation which should work well
enough.
"""
data = np.array([0.92961609, 0.31637555, 0.18391881, 0.20456028,
0.56772503, 0.5955447, 0.96451452, 0.6531771,
0.74890664, 0.65356987])
dt = 1.0
# Scenario 1.
new_dt = 0.45
a = 1
expected_output = np.array([
0.92961609, 0.55712768, 0.31720733, 0.24275977, 0.17825931,
0.16750234, 0.17561933, 0.20626905, 0.37726064, 0.5647072,
0.47145546, 0.59222238, 0.58665834, 0.91241347, 0.79909224,
0.61631275, 0.61258393, 0.61611633, 0.73239733, 0.56371682,
0.65356987])
output = lanczos_interpolation(
data, old_dt=dt, new_start=0.0, old_start=0.0, new_dt=new_dt,
new_npts=21, a=a)
np.testing.assert_allclose(output, expected_output, atol=1E-9)
# Scenario 2.
new_dt = 0.72
a = 12
expected_output = np.array([
0.92961609, 0.54632548, 0.14335148, 0.19675436, 0.19030867,
0.41722415, 0.60644459, 0.6018648, 0.88751628, 0.90970863,
0.58602723, 0.71521445, 0.83288791])
output = lanczos_interpolation(
data, old_dt=dt, new_start=0.0, old_start=0.0, new_dt=new_dt,
new_npts=13, a=a)
np.testing.assert_allclose(output, expected_output, atol=1E-9)
def test_lanczos_interpolation_units(self):
"""
Regression test for a bug that manifested when the original sampling
rate is not 1 Hertz and new and old start times are not identical.
"""
# Generate a highly oversampled sine curve. Upsample and downsample
# it again and it should not change a lot except at the boundaries.
# Also shift a bit to trigger the bug.
original_dt = 13.333
new_dt = 17.23
data = np.sin(np.linspace(0, 2 * np.pi, 1000))
output = lanczos_interpolation(
data, old_dt=original_dt, new_start=10 * original_dt,
old_start=0.0, a=20,
new_dt=new_dt, new_npts=int(990 * original_dt / new_dt))
output = lanczos_interpolation(
output, old_dt=new_dt, new_start=10 * original_dt,
old_start=0.0, a=20,
new_dt=original_dt, new_npts=int(980 * original_dt / new_dt) - 1)
np.testing.assert_allclose(data[220:620], output[200:600], atol=1E-4,
rtol=1E-4)
def test_plot_lanczos_window(self):
"""
Tests the plot_lanczos_window function.
"""
with ImageComparison(self.path_images,
'plot_lanczos_window.png') as ic:
plt.figure(figsize=(8, 12))
plot_lanczos_windows(a=20, filename=ic.name)
def suite():
return unittest.makeSuite(InterpolationTestCase, 'test')
if __name__ == '__main__':
unittest.main(defaultTest='suite')