← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/fftpack/tests/test_pseudo_diffs.orig

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/fftpack/tests/test_pseudo_diffs.orig
 
1
#!/usr/bin/env python
 
2
# Created by Pearu Peterson, September 2002
 
3
""" Test functions for fftpack.pseudo_diffs module
 
4
"""
 
5
__usage__ = """
 
6
Build fftpack:
 
7
  python setup_fftpack.py build
 
8
Run tests if scipy is installed:
 
9
  python -c 'import scipy;scipy.fftpack.test(<level>)'
 
10
Run tests if fftpack is not installed:
 
11
  python tests/test_pseudo_diffs.py [<level>]
 
12
"""
 
13
import sys
 
14
from numpy.testing import *
 
15
set_package_path()
 
16
from fftpack import diff,fft,ifft,tilbert,itilbert,hilbert,ihilbert,rfft
 
17
from fftpack import shift
 
18
from fftpack import fftfreq
 
19
restore_path()
 
20
 
 
21
from numpy import arange, add, array, sin, cos, pi,exp,tanh,sum,sign
 
22
 
 
23
def random(size):
 
24
    return rand(*size)
 
25
 
 
26
def direct_diff(x,k=1,period=None):
 
27
    fx = fft(x)
 
28
    n = len (fx)
 
29
    if period is None:
 
30
        period = 2*pi
 
31
    w = fftfreq(n)*2j*pi/period*n
 
32
    if k<0:
 
33
        w = 1 / w**k
 
34
        w[0] = 0.0
 
35
    else:
 
36
        w = w**k
 
37
    if n>2000:
 
38
        w[250:n-250] = 0.0
 
39
    return ifft(w*fx).real
 
40
 
 
41
def direct_tilbert(x,h=1,period=None):
 
42
    fx = fft(x)
 
43
    n = len (fx)
 
44
    if period is None:
 
45
        period = 2*pi
 
46
    w = fftfreq(n)*h*2*pi/period*n
 
47
    w[0] = 1
 
48
    w = 1j/tanh(w)
 
49
    w[0] = 0j
 
50
    return ifft(w*fx)
 
51
 
 
52
def direct_itilbert(x,h=1,period=None):
 
53
    fx = fft(x)
 
54
    n = len (fx)
 
55
    if period is None:
 
56
        period = 2*pi
 
57
    w = fftfreq(n)*h*2*pi/period*n
 
58
    w = -1j*tanh(w)
 
59
    return ifft(w*fx)
 
60
 
 
61
def direct_hilbert(x):
 
62
    fx = fft(x)
 
63
    n = len (fx)
 
64
    w = fftfreq(n)*n
 
65
    w = 1j*sign(w)
 
66
    return ifft(w*fx)
 
67
 
 
68
def direct_ihilbert(x):
 
69
    return -direct_hilbert(x)
 
70
 
 
71
def direct_shift(x,a,period=None):
 
72
    n = len(x)
 
73
    if period is None:
 
74
        k = fftfreq(n)*1j*n
 
75
    else:
 
76
        k = fftfreq(n)*2j*pi/period*n
 
77
    return ifft(fft(x)*exp(k*a)).real
 
78
 
 
79
 
 
80
class test_diff(ScipyTestCase):
 
81
 
 
82
    def check_definition(self):
 
83
        for n in [16,17,64,127,32]:
 
84
            x = arange(n)*2*pi/n
 
85
            assert_array_almost_equal(diff(sin(x)),direct_diff(sin(x)))
 
86
            assert_array_almost_equal(diff(sin(x),2),direct_diff(sin(x),2))
 
87
            assert_array_almost_equal(diff(sin(x),3),direct_diff(sin(x),3))
 
88
            assert_array_almost_equal(diff(sin(x),4),direct_diff(sin(x),4))
 
89
            assert_array_almost_equal(diff(sin(x),5),direct_diff(sin(x),5))
 
90
            assert_array_almost_equal(diff(sin(2*x),3),direct_diff(sin(2*x),3))
 
91
            assert_array_almost_equal(diff(sin(2*x),4),direct_diff(sin(2*x),4))
 
92
            assert_array_almost_equal(diff(cos(x)),direct_diff(cos(x)))
 
93
            assert_array_almost_equal(diff(cos(x),2),direct_diff(cos(x),2))
 
94
            assert_array_almost_equal(diff(cos(x),3),direct_diff(cos(x),3))
 
95
            assert_array_almost_equal(diff(cos(x),4),direct_diff(cos(x),4))
 
96
            assert_array_almost_equal(diff(cos(2*x)),direct_diff(cos(2*x)))
 
97
            assert_array_almost_equal(diff(sin(x*n/8)),direct_diff(sin(x*n/8)))
 
98
            assert_array_almost_equal(diff(cos(x*n/8)),direct_diff(cos(x*n/8)))
 
99
            for k in range(5):
 
100
                assert_array_almost_equal(diff(sin(4*x),k),direct_diff(sin(4*x),k))
 
101
                assert_array_almost_equal(diff(cos(4*x),k),direct_diff(cos(4*x),k))
 
102
 
 
103
    def check_period(self):
 
104
        for n in [17,64]:
 
105
            x = arange(n)/float(n)
 
106
            assert_array_almost_equal(diff(sin(2*pi*x),period=1),
 
107
                                      2*pi*cos(2*pi*x))
 
108
            assert_array_almost_equal(diff(sin(2*pi*x),3,period=1),
 
109
                                      -(2*pi)**3*cos(2*pi*x))
 
110
 
 
111
    def check_sin(self):
 
112
        for n in [32,64,77]:
 
113
            x = arange(n)*2*pi/n
 
114
            assert_array_almost_equal(diff(sin(x)),cos(x))
 
115
            assert_array_almost_equal(diff(cos(x)),-sin(x))
 
116
            assert_array_almost_equal(diff(sin(x),2),-sin(x))
 
117
            assert_array_almost_equal(diff(sin(x),4),sin(x))
 
118
            assert_array_almost_equal(diff(sin(4*x)),4*cos(4*x))
 
119
            assert_array_almost_equal(diff(sin(sin(x))),cos(x)*cos(sin(x)))
 
120
 
 
121
    def check_expr(self):
 
122
        for n in [64,77,100,128,256,512,1024,2048,4096,8192][:5]:
 
123
            x = arange(n)*2*pi/n
 
124
            f=sin(x)*cos(4*x)+exp(sin(3*x))
 
125
            df=cos(x)*cos(4*x)-4*sin(x)*sin(4*x)+3*cos(3*x)*exp(sin(3*x))
 
126
            ddf=-17*sin(x)*cos(4*x)-8*cos(x)*sin(4*x)\
 
127
                 -9*sin(3*x)*exp(sin(3*x))+9*cos(3*x)**2*exp(sin(3*x))
 
128
            d1 = diff(f)
 
129
            assert_array_almost_equal(d1,df)
 
130
            assert_array_almost_equal(diff(df),ddf)
 
131
            assert_array_almost_equal(diff(f,2),ddf)
 
132
            assert_array_almost_equal(diff(ddf,-1),df)
 
133
            #print max(abs(d1-df))
 
134
 
 
135
    def check_expr_large(self):
 
136
        for n in [2048,4096]:
 
137
            x = arange(n)*2*pi/n
 
138
            f=sin(x)*cos(4*x)+exp(sin(3*x))
 
139
            df=cos(x)*cos(4*x)-4*sin(x)*sin(4*x)+3*cos(3*x)*exp(sin(3*x))
 
140
            ddf=-17*sin(x)*cos(4*x)-8*cos(x)*sin(4*x)\
 
141
                 -9*sin(3*x)*exp(sin(3*x))+9*cos(3*x)**2*exp(sin(3*x))
 
142
            assert_array_almost_equal(diff(f),df)
 
143
            assert_array_almost_equal(diff(df),ddf)
 
144
            assert_array_almost_equal(diff(ddf,-1),df)
 
145
            assert_array_almost_equal(diff(f,2),ddf)
 
146
 
 
147
    def check_int(self):
 
148
        n = 64
 
149
        x = arange(n)*2*pi/n
 
150
        assert_array_almost_equal(diff(sin(x),-1),-cos(x))
 
151
        assert_array_almost_equal(diff(sin(x),-2),-sin(x))
 
152
        assert_array_almost_equal(diff(sin(x),-4),sin(x))
 
153
        assert_array_almost_equal(diff(2*cos(2*x),-1),sin(2*x))
 
154
 
 
155
    def check_random_even(self):
 
156
        for k in [0,2,4,6]:
 
157
            for n in [60,32,64,56,55]:
 
158
                f=random ((n,))
 
159
                af=sum(f)/n
 
160
                f=f-af
 
161
                # zeroing Nyquist mode:
 
162
                f = diff(diff(f,1),-1)
 
163
                assert_almost_equal(sum(f),0.0)
 
164
                assert_array_almost_equal(diff(diff(f,k),-k),f)
 
165
                assert_array_almost_equal(diff(diff(f,-k),k),f)
 
166
 
 
167
    def check_random_odd(self):
 
168
        for k in [0,1,2,3,4,5,6]:
 
169
            for n in [33,65,55]:
 
170
                f=random ((n,))
 
171
                af=sum(f)/n
 
172
                f=f-af
 
173
                assert_almost_equal(sum(f),0.0)
 
174
                assert_array_almost_equal(diff(diff(f,k),-k),f)
 
175
                assert_array_almost_equal(diff(diff(f,-k),k),f)
 
176
 
 
177
    def check_zero_nyquist (self):
 
178
        for k in [0,1,2,3,4,5,6]:
 
179
            for n in [32,33,64,56,55]:
 
180
                f=random ((n,))
 
181
                af=sum(f)/n
 
182
                f=f-af
 
183
                # zeroing Nyquist mode:
 
184
                f = diff(diff(f,1),-1)
 
185
                assert_almost_equal(sum(f),0.0)
 
186
                assert_array_almost_equal(diff(diff(f,k),-k),f)
 
187
                assert_array_almost_equal(diff(diff(f,-k),k),f)
 
188
 
 
189
    def bench_random(self,level=5):
 
190
        print
 
191
        print 'Differentiation of periodic functions'
 
192
        print '====================================='
 
193
        print ' size  |  convolve |    naive'
 
194
        print '-------------------------------------'
 
195
        for size,repeat in [(100,1500),(1000,300),
 
196
                            (256,1500),
 
197
                            (512,1000),
 
198
                            (1024,500),
 
199
                            (2048,200),
 
200
                            (2048*2,100),
 
201
                            (2048*4,50),
 
202
                            ]:
 
203
            print '%6s' % size,
 
204
            sys.stdout.flush()
 
205
            x = arange (size)*2*pi/size
 
206
            if size<2000:
 
207
                f = sin(x)*cos(4*x)+exp(sin(3*x))
 
208
            else:
 
209
                f = sin(x)*cos(4*x)
 
210
            assert_array_almost_equal(diff(f,1),direct_diff(f,1))
 
211
            assert_array_almost_equal(diff(f,2),direct_diff(f,2))
 
212
            print '| %9.2f' % self.measure('diff(f,3)',repeat),
 
213
            sys.stdout.flush()
 
214
            print '| %9.2f' % self.measure('direct_diff(f,3)',repeat),
 
215
            sys.stdout.flush()
 
216
            print ' (secs for %s calls)' % (repeat)
 
217
 
 
218
 
 
219
class test_tilbert(ScipyTestCase):
 
220
 
 
221
    def check_definition(self):
 
222
        for h in [0.1,0.5,1,5.5,10]:
 
223
            for n in [16,17,64,127]:
 
224
                x = arange(n)*2*pi/n
 
225
                y = tilbert(sin(x),h)
 
226
                y1 = direct_tilbert(sin(x),h)
 
227
                assert_array_almost_equal (y,y1)
 
228
                assert_array_almost_equal(tilbert(sin(x),h),
 
229
                                          direct_tilbert(sin(x),h))
 
230
                assert_array_almost_equal(tilbert(sin(2*x),h),
 
231
                                          direct_tilbert(sin(2*x),h))
 
232
 
 
233
    def check_random_even(self):
 
234
        for h in [0.1,0.5,1,5.5,10]:
 
235
            for n in [32,64,56]:
 
236
                f=random ((n,))
 
237
                af=sum(f)/n
 
238
                f=f-af
 
239
                assert_almost_equal(sum(f),0.0)
 
240
                assert_array_almost_equal(direct_tilbert(direct_itilbert(f,h),h),f)
 
241
 
 
242
    def check_random_odd(self):
 
243
        for h in [0.1,0.5,1,5.5,10]:
 
244
            for n in [33,65,55]:
 
245
                f=random ((n,))
 
246
                af=sum(f)/n
 
247
                f=f-af
 
248
                assert_almost_equal(sum(f),0.0)
 
249
                assert_array_almost_equal(itilbert(tilbert(f,h),h),f)
 
250
                assert_array_almost_equal(tilbert(itilbert(f,h),h),f)
 
251
 
 
252
    def bench_random(self,level=5):
 
253
        print
 
254
        print ' Tilbert transform of periodic functions'
 
255
        print '========================================='
 
256
        print ' size  | optimized |    naive'
 
257
        print '-----------------------------------------'
 
258
        for size,repeat in [(100,1500),(1000,300),
 
259
                            (256,1500),
 
260
                            (512,1000),
 
261
                            (1024,500),
 
262
                            (2048,200),
 
263
                            (2048*2,100),
 
264
                            (2048*4,50),
 
265
                            ]:
 
266
            print '%6s' % size,
 
267
            sys.stdout.flush()
 
268
            x = arange (size)*2*pi/size
 
269
            if size<2000:
 
270
                f = sin(x)*cos(4*x)+exp(sin(3*x))
 
271
            else:
 
272
                f = sin(x)*cos(4*x)
 
273
            assert_array_almost_equal(tilbert(f,1),direct_tilbert(f,1))
 
274
            print '| %9.2f' % self.measure('tilbert(f,1)',repeat),
 
275
            sys.stdout.flush()
 
276
            print '| %9.2f' % self.measure('direct_tilbert(f,1)',repeat),
 
277
            sys.stdout.flush()
 
278
            print ' (secs for %s calls)' % (repeat)
 
279
 
 
280
class test_itilbert(ScipyTestCase):
 
281
 
 
282
    def check_definition(self):
 
283
        for h in [0.1,0.5,1,5.5,10]:
 
284
            for n in [16,17,64,127]:
 
285
                x = arange(n)*2*pi/n
 
286
                y = itilbert(sin(x),h)
 
287
                y1 = direct_itilbert(sin(x),h)
 
288
                assert_array_almost_equal (y,y1)
 
289
                assert_array_almost_equal(itilbert(sin(x),h),
 
290
                                          direct_itilbert(sin(x),h))
 
291
                assert_array_almost_equal(itilbert(sin(2*x),h),
 
292
                                          direct_itilbert(sin(2*x),h))
 
293
 
 
294
class test_hilbert(ScipyTestCase):
 
295
 
 
296
    def check_definition(self):
 
297
        for n in [16,17,64,127]:
 
298
            x = arange(n)*2*pi/n
 
299
            y = hilbert(sin(x))
 
300
            y1 = direct_hilbert(sin(x))
 
301
            assert_array_almost_equal (y,y1)
 
302
            assert_array_almost_equal(hilbert(sin(2*x)),
 
303
                                      direct_hilbert(sin(2*x)))
 
304
 
 
305
    def check_tilbert_relation(self):
 
306
        for n in [16,17,64,127]:
 
307
            x = arange(n)*2*pi/n
 
308
            f = sin (x)+cos (2*x)*sin(x)
 
309
            y = hilbert(f)
 
310
            y1 = direct_hilbert(f)
 
311
            assert_array_almost_equal (y,y1)
 
312
            y2 = tilbert(f,h=10)
 
313
            assert_array_almost_equal (y,y2)
 
314
 
 
315
    def check_random_odd(self):
 
316
        for n in [33,65,55]:
 
317
            f=random ((n,))
 
318
            af=sum(f)/n
 
319
            f=f-af
 
320
            assert_almost_equal(sum(f),0.0)
 
321
            assert_array_almost_equal(ihilbert(hilbert(f)),f)
 
322
            assert_array_almost_equal(hilbert(ihilbert(f)),f)
 
323
 
 
324
    def check_random_even(self):
 
325
        for n in [32,64,56]:
 
326
            f=random ((n,))
 
327
            af=sum(f)/n
 
328
            f=f-af
 
329
            # zeroing Nyquist mode:
 
330
            f = diff(diff(f,1),-1)
 
331
            assert_almost_equal(sum(f),0.0)
 
332
            assert_array_almost_equal(direct_hilbert(direct_ihilbert(f)),f)
 
333
            assert_array_almost_equal(hilbert(ihilbert(f)),f)
 
334
 
 
335
    def bench_random(self,level=5):
 
336
        print
 
337
        print ' Hilbert transform of periodic functions'
 
338
        print '========================================='
 
339
        print ' size  | optimized |    naive'
 
340
        print '-----------------------------------------'
 
341
        for size,repeat in [(100,1500),(1000,300),
 
342
                            (256,1500),
 
343
                            (512,1000),
 
344
                            (1024,500),
 
345
                            (2048,200),
 
346
                            (2048*2,100),
 
347
                            (2048*4,50),
 
348
                            ]:
 
349
            print '%6s' % size,
 
350
            sys.stdout.flush()
 
351
            x = arange (size)*2*pi/size
 
352
            if size<2000:
 
353
                f = sin(x)*cos(4*x)+exp(sin(3*x))
 
354
            else:
 
355
                f = sin(x)*cos(4*x)
 
356
            assert_array_almost_equal(hilbert(f),direct_hilbert(f))
 
357
            print '| %9.2f' % self.measure('hilbert(f)',repeat),
 
358
            sys.stdout.flush()
 
359
            print '| %9.2f' % self.measure('direct_hilbert(f)',repeat),
 
360
            sys.stdout.flush()
 
361
            print ' (secs for %s calls)' % (repeat)
 
362
 
 
363
class test_ihilbert(ScipyTestCase):
 
364
 
 
365
    def check_definition(self):
 
366
        for n in [16,17,64,127]:
 
367
            x = arange(n)*2*pi/n
 
368
            y = ihilbert(sin(x))
 
369
            y1 = direct_ihilbert(sin(x))
 
370
            assert_array_almost_equal (y,y1)
 
371
            assert_array_almost_equal(ihilbert(sin(2*x)),
 
372
                                      direct_ihilbert(sin(2*x)))
 
373
 
 
374
    def check_itilbert_relation(self):
 
375
        for n in [16,17,64,127]:
 
376
            x = arange(n)*2*pi/n
 
377
            f = sin (x)+cos (2*x)*sin(x)
 
378
            y = ihilbert(f)
 
379
            y1 = direct_ihilbert(f)
 
380
            assert_array_almost_equal (y,y1)
 
381
            y2 = itilbert(f,h=10)
 
382
            assert_array_almost_equal (y,y2)
 
383
 
 
384
class test_shift(ScipyTestCase):
 
385
 
 
386
    def check_definition(self):
 
387
        for n in [18,17,64,127,32,2048,256]:
 
388
            x = arange(n)*2*pi/n
 
389
            for a in [0.1,3]:
 
390
                assert_array_almost_equal(shift(sin(x),a),direct_shift(sin(x),a))
 
391
                assert_array_almost_equal(shift(sin(x),a),sin(x+a))
 
392
                assert_array_almost_equal(shift(cos(x),a),cos(x+a))
 
393
                assert_array_almost_equal(shift(cos(2*x)+sin(x),a),
 
394
                                          cos(2*(x+a))+sin(x+a))
 
395
                assert_array_almost_equal(shift(exp(sin(x)),a),exp(sin(x+a)))
 
396
            assert_array_almost_equal(shift(sin(x),2*pi),sin(x))
 
397
            assert_array_almost_equal(shift(sin(x),pi),-sin(x))
 
398
            assert_array_almost_equal(shift(sin(x),pi/2),cos(x))
 
399
 
 
400
    def bench_random(self,level=5):
 
401
        print
 
402
        print ' Shifting periodic functions'
 
403
        print '=============================='
 
404
        print ' size  | optimized |    naive'
 
405
        print '------------------------------'
 
406
        for size,repeat in [(100,1500),(1000,300),
 
407
                            (256,1500),
 
408
                            (512,1000),
 
409
                            (1024,500),
 
410
                            (2048,200),
 
411
                            (2048*2,100),
 
412
                            (2048*4,50),
 
413
                            ]:
 
414
            print '%6s' % size,
 
415
            sys.stdout.flush()
 
416
            x = arange (size)*2*pi/size
 
417
            a = 1
 
418
            if size<2000:
 
419
                f = sin(x)*cos(4*x)+exp(sin(3*x))
 
420
                sf = sin(x+a)*cos(4*(x+a))+exp(sin(3*(x+a)))
 
421
            else:
 
422
                f = sin(x)*cos(4*x)
 
423
                sf = sin(x+a)*cos(4*(x+a))
 
424
            assert_array_almost_equal(direct_shift(f,1),sf)
 
425
            assert_array_almost_equal(shift(f,1),sf)
 
426
            print '| %9.2f' % self.measure('shift(f,a)',repeat),
 
427
            sys.stdout.flush()
 
428
            print '| %9.2f' % self.measure('direct_shift(f,a)',repeat),
 
429
            sys.stdout.flush()
 
430
            print ' (secs for %s calls)' % (repeat)
 
431
 
 
432
if __name__ == "__main__":
 
433
    ScipyTest('fftpack.pseudo_diffs').run()