forked from JuliaLang/julia
-
Notifications
You must be signed in to change notification settings - Fork 0
/
dsp.jl
116 lines (98 loc) · 5.1 KB
/
dsp.jl
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
# Filter
b = [1., 2., 3., 4.]
x = [1., 1., 0., 1., 1., 0., 0., 0.]
@test filt(b, 1., x) == [1., 3., 5., 8., 7., 5., 7., 4.]
@test filt(b, [1., -0.5], x) == [1., 3.5, 6.75, 11.375, 12.6875, 11.34375, 12.671875, 10.3359375]
# With ranges
@test filt(b, 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.]
@test filt(1.:4., 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.]
# Across an array is the same as channel-by-channel
@test filt(b, 1., [x 1.0:8.0]) == [filt(b, 1., x) filt(b, 1., 1.0:8.0)]
@test filt(b, [1., -0.5], [x 1.0:8.0]) == [filt(b, [1., -0.5], x) filt(b, [1., -0.5], 1.0:8.0)]
si = zeros(3)
@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, si) filt(b, 1., 1.0:8.0, si)]
@test si == zeros(3) # Will likely fail if/when arrayviews are implemented
si = [zeros(3) ones(3)]
@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, zeros(3)) filt(b, 1., 1.0:8.0, ones(3))]
# With initial conditions: a lowpass 5-pole butterworth filter with W_n = 0.25,
# and a stable initial filter condition matched to the initial value
b = [0.003279216306360201,0.016396081531801006,0.03279216306360201,0.03279216306360201,0.016396081531801006,0.003279216306360201]
a = [1.0,-2.4744161749781606,2.8110063119115782,-1.703772240915465,0.5444326948885326,-0.07231566910295834]
si = [0.9967207836936347,-1.4940914728163142,1.2841226760316475,-0.4524417279474106,0.07559488540931815]
@test_approx_eq filt(b, a, ones(10), si) ones(10) # Shouldn't affect DC offset
# Convolution
a = [1., 2., 1., 2.]
b = [1., 2., 3.]
@test_approx_eq conv(a, b) [1., 4., 8., 10., 7., 6.]
@test_approx_eq conv(complex(a, ones(4)), complex(b)) complex([1., 4., 8., 10., 7., 6.], [1., 3., 6., 6., 5., 3.])
# Discrete cosine transform (DCT) tests
if Base.fftw_vendor() != :mkl
a = rand(8,11) + im*rand(8,11)
@test norm(idct(dct(a)) - a) < 1e-8
X = reshape([1,2,7,2,1,5,9,-1,3,4,6,9],3,4)
Y = rand(17,14)
Y[3:5,9:12] = X
sX = slice(Y,3:5,9:12)
true_Xdct = [ 13.856406460551018 -3.863239728836245 2.886751345948129 -0.274551994240164; -2.828427124746190 -2.184015211898548 -4.949747468305834 3.966116180118245; 4.898979485566356 -0.194137576915510 -2.857738033247041 2.731723009609389 ]
true_Xdct_1 = [ 5.773502691896258 4.618802153517007 6.350852961085884 10.969655114602890; -4.242640687119286 -2.121320343559643 4.242640687119286 -3.535533905932738; 1.632993161855452 2.041241452319315 5.715476066494083 0.408248290463863 ]
true_Xdct_2 = [ 8. -3.854030797826254 -3.0 3.761176226848022;
4.0 -2.071929829606556 4.0 -2.388955165168770; 12. -0.765366864730179 4.0 -1.847759065022573 ]
Xdct = dct(X)
Xdct! = float(X); dct!(Xdct!)
Xdct_1 = dct(X,1)
Xdct!_1 = float(X); dct!(Xdct!_1,1)
Xdct_2 = dct(X,2)
Xdct!_2 = float(X); dct!(Xdct!_2,2)
Xidct = idct(true_Xdct)
Xidct! = copy(true_Xdct); idct!(Xidct!)
Xidct_1 = idct(true_Xdct_1,1)
Xidct!_1 = copy(true_Xdct_1); idct!(Xidct!_1,1)
Xidct_2 = idct(true_Xdct_2,2)
Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2)
pXdct = plan_dct(X)(X)
pXdct! = float(X); plan_dct!(pXdct!)(pXdct!)
pXdct_1 = plan_dct(X,1)(X)
pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)(pXdct!_1)
pXdct_2 = plan_dct(X,2)(X)
pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)(pXdct!_2)
pXidct = plan_idct(true_Xdct)(true_Xdct)
pXidct! = copy(true_Xdct); plan_idct!(pXidct!)(pXidct!)
pXidct_1 = plan_idct(true_Xdct_1,1)(true_Xdct_1)
pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)(pXidct!_1)
pXidct_2 = plan_idct(true_Xdct_2,2)(true_Xdct_2)
pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)(pXidct!_2)
sXdct = dct(sX)
psXdct = plan_dct(sX)(sX)
sYdct! = copy(Y); sXdct! = slice(sYdct!,3:5,9:12); dct!(sXdct!)
psYdct! = copy(Y); psXdct! = slice(psYdct!,3:5,9:12); plan_dct!(psXdct!)(psXdct!)
for i = 1:length(X)
@test_approx_eq Xdct[i] true_Xdct[i]
@test_approx_eq Xdct![i] true_Xdct[i]
@test_approx_eq Xdct_1[i] true_Xdct_1[i]
@test_approx_eq Xdct!_1[i] true_Xdct_1[i]
@test_approx_eq Xdct_2[i] true_Xdct_2[i]
@test_approx_eq Xdct!_2[i] true_Xdct_2[i]
@test_approx_eq pXdct[i] true_Xdct[i]
@test_approx_eq pXdct![i] true_Xdct[i]
@test_approx_eq pXdct_1[i] true_Xdct_1[i]
@test_approx_eq pXdct!_1[i] true_Xdct_1[i]
@test_approx_eq pXdct_2[i] true_Xdct_2[i]
@test_approx_eq pXdct!_2[i] true_Xdct_2[i]
@test_approx_eq Xidct[i] X[i]
@test_approx_eq Xidct![i] X[i]
@test_approx_eq Xidct_1[i] X[i]
@test_approx_eq Xidct!_1[i] X[i]
@test_approx_eq Xidct_2[i] X[i]
@test_approx_eq Xidct!_2[i] X[i]
@test_approx_eq pXidct[i] X[i]
@test_approx_eq pXidct![i] X[i]
@test_approx_eq pXidct_1[i] X[i]
@test_approx_eq pXidct!_1[i] X[i]
@test_approx_eq pXidct_2[i] X[i]
@test_approx_eq pXidct!_2[i] X[i]
@test_approx_eq sXdct[i] true_Xdct[i]
@test_approx_eq psXdct[i] true_Xdct[i]
@test_approx_eq sXdct![i] true_Xdct[i]
@test_approx_eq psXdct![i] true_Xdct[i]
end
end # fftw_vendor() != :mkl