##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
- Installed: sp (signal processing package for python)...
yamaro -
r18:19
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,59
1 '''
2 Created on May 29, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 #from sp import multirate
8 import cshift
9 from multirate import upfirdn
10
11 def afb(x, af):
12
13 # Analysis filter bank
14 #
15 # USAGE:
16 # [lo, hi] = afb(x, af)
17 # INPUT:
18 # x - N-point vector, where
19 # 1) N is even
20 # 2) N >= length(af)
21 # af - analysis filters
22 # af(:, 1) - lowpass filter (even length)
23 # af(:, 2) - highpass filter (even length)
24 # OUTPUT:
25 # lo - Low frequecy output
26 # hi - High frequency output
27 # EXAMPLE:
28 # [af, sf] = farras;
29 # x = rand(1,64);
30 # [lo, hi] = afb(x, af);
31 # y = sfb(lo, hi, sf);
32 # err = x - y;
33 # max(abs(err))
34 #
35 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
36 # http://taco.poly.edu/WaveletSoftware/
37
38 N = x.size;
39 L = (af).size/2;
40 x = cshift(x,-L);
41
42 # lowpass filter
43 lo = upfirdn(x, af[:,0], 1, 2);
44
45 # VERIFY THIS!!!!!!!!!!!!
46 for i in range(0, L):
47 lo[i] = lo[N/2+[i]] + lo[i];
48
49 lo = lo[1:N/2];
50
51 # highpass filter
52 hi = upfirdn(x, af[:,2], 1, 2);
53
54 for j in range(0, L):
55 hi[j] = hi(N/2+[j]) + hi[j];
56
57 hi = hi[1:N/2];
58
59 return lo, hi
Show file before | Show file after | Hide comments
@@ -0,0 +1,29
1 '''
2 Created on May 30, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import numpy as np
8
9 def cshift(x, m):
10
11 # Circular Shift
12 #
13 # USAGE:
14 # y = cshift(x, m)
15 # INPUT:
16 # x - N-point vector
17 # m - amount of shift
18 # OUTPUT:
19 # y - vector x will be shifted by m samples to the left
20 #
21 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
22 # http://taco.poly.edu/WaveletSoftware/
23
24 N = x.size;
25 n = np.arange(N-1);
26 n = np.mod(n-m, N);
27 y = x[n];
28
29 return y No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,63
1 '''
2 Created on May 29, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import numpy as np
8 import afb
9
10 def dualtree(x, J, Faf, af):
11
12 # Dual-tree Complex Discrete Wavelet Transform
13 #
14 # USAGE:
15 # w = dualtree(x, J, Faf, af)
16 # INPUT:
17 # x - N-point vector
18 # 1) N is divisible by 2^J
19 # 2) N >= 2^(J-1)*length(af)
20 # J - number of stages
21 # Faf - filters for the first stage
22 # af - filters for the remaining stages
23 # OUTPUT:
24 # w - DWT coefficients
25 # w{j}{1}, j = 1..J - real part
26 # w{j}{2}, j = 1..J - imaginary part
27 # w{J+1}{d} - lowpass coefficients, d = 1,2
28 # EXAMPLE:
29 # x = rand(1, 512);
30 # J = 4;
31 # [Faf, Fsf] = FSfarras;
32 # [af, sf] = dualfilt1;
33 # w = dualtree(x, J, Faf, af);
34 # y = idualtree(w, J, Fsf, sf);
35 # err = x - y;
36 # max(abs(err))
37 #
38 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
39 # http://taco.poly.edu/WaveletSoftware/
40
41 # normalization
42 x = x/np.sqrt(2);
43
44 w = np.zeros(shape=(J,2)) ### VERIFY THIS DEFINITION
45
46 # Tree 1
47 [x1,w[1,0]] = afb(x, Faf[0,1]); # w{1}{1}
48
49 for j in range (2,J):
50 [x1,w[j,0]] = afb(x1, af[0,1]); #check this
51
52 w[J+1,1] = x1;
53
54 # Tree 2
55 [x2,w[1,2]] = afb(x, Faf[0,1]);
56
57 for j in range (2,J):
58 [x2,w[j,2]] = afb(x2, af[0,1]);
59
60 w[J+1,2] = x2;
61
62 return w
63
Show file before | Show file after | Hide comments
@@ -0,0 +1,238
1 #!/usr/bin/env python
2 """Module providing Multirate signal processing functionality.
3 Largely based on MATLAB's Multirate signal processing toolbox with consultation
4 of Octave m-file source code.
5
6 From: https://github.com/mubeta06/python/blob/master/signal_processing/sp/multirate.py
7 """
8 #import sys
9 import fractions
10 import numpy
11 from scipy import signal
12
13
14 def downsample(s, n, phase=0):
15 """Decrease sampling rate by integer factor n with included offset phase.
16 """
17 return s[phase::n]
18
19
20 def upsample(s, n, phase=0):
21 """Increase sampling rate by integer factor n with included offset phase.
22 """
23 return numpy.roll(numpy.kron(s, numpy.r_[1, numpy.zeros(n-1)]), phase)
24
25
26 def decimate(s, r, n=None, fir=False):
27 """Decimation - decrease sampling rate by r. The decimation process filters
28 the input data s with an order n lowpass filter and then resamples the
29 resulting smoothed signal at a lower rate. By default, decimate employs an
30 eighth-order lowpass Chebyshev Type I filter with a cutoff frequency of
31 0.8/r. It filters the input sequence in both the forward and reverse
32 directions to remove all phase distortion, effectively doubling the filter
33 order. If 'fir' is set to True decimate uses an order 30 FIR filter (by
34 default otherwise n), instead of the Chebyshev IIR filter. Here decimate
35 filters the input sequence in only one direction. This technique conserves
36 memory and is useful for working with long sequences.
37 """
38 if fir:
39 if n is None:
40 n = 30
41 b = signal.firwin(n, 1.0/r)
42 a = 1
43 f = signal.lfilter(b, a, s)
44 else: #iir
45 if n is None:
46 n = 8
47 b, a = signal.cheby1(n, 0.05, 0.8/r)
48 f = signal.filtfilt(b, a, s)
49 return downsample(f, r)
50
51
52 def interp(s, r, l=4, alpha=0.5):
53 """Interpolation - increase sampling rate by integer factor r. Interpolation
54 increases the original sampling rate for a sequence to a higher rate. interp
55 performs lowpass interpolation by inserting zeros into the original sequence
56 and then applying a special lowpass filter. l specifies the filter length
57 and alpha the cut-off frequency. The length of the FIR lowpass interpolating
58 filter is 2*l*r+1. The number of original sample values used for
59 interpolation is 2*l. Ordinarily, l should be less than or equal to 10. The
60 original signal is assumed to be band limited with normalized cutoff
61 frequency 0=alpha=1, where 1 is half the original sampling frequency (the
62 Nyquist frequency). The default value for l is 4 and the default value for
63 alpha is 0.5.
64 """
65 b = signal.firwin(2*l*r+1, alpha/r);
66 a = 1
67 return r*signal.lfilter(b, a, upsample(s, r))[r*l+1:-1]
68
69
70 def resample(s, p, q, h=None):
71 """Change sampling rate by rational factor. This implementation is based on
72 the Octave implementation of the resample function. It designs the
73 anti-aliasing filter using the window approach applying a Kaiser window with
74 the beta term calculated as specified by [2].
75 Ref [1] J. G. Proakis and D. G. Manolakis,
76 Digital Signal Processing: Principles, Algorithms, and Applications,
77 4th ed., Prentice Hall, 2007. Chap. 6
78
79 Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
80 Discrete-time signal processing, Signal processing series,
81 Prentice-Hall, 1999
82 """
83 gcd = fractions.gcd(p,q)
84 if gcd>1:
85 p=p/gcd
86 q=q/gcd
87
88 if h is None: #design filter
89 #properties of the antialiasing filter
90 log10_rejection = -3.0
91 stopband_cutoff_f = 1.0/(2.0 * max(p,q))
92 roll_off_width = stopband_cutoff_f / 10.0
93
94 #determine filter length
95 #use empirical formula from [2] Chap 7, Eq. (7.63) p 476
96 rejection_db = -20.0*log10_rejection;
97 l = numpy.ceil((rejection_db-8.0) / (28.714 * roll_off_width))
98
99 #ideal sinc filter
100 t = numpy.arange(-l, l + 1)
101 ideal_filter=2*p*stopband_cutoff_f*numpy.sinc(2*stopband_cutoff_f*t)
102
103 #determine parameter of Kaiser window
104 #use empirical formula from [2] Chap 7, Eq. (7.62) p 474
105 beta = signal.kaiser_beta(rejection_db)
106
107 #apodize ideal filter response
108 h = numpy.kaiser(2*l+1, beta)*ideal_filter
109
110 ls = len(s)
111 lh = len(h)
112
113 l = (lh - 1)/2.0
114 ly = numpy.ceil(ls*p/float(q))
115
116 #pre and postpad filter response
117 nz_pre = numpy.floor(q - numpy.mod(l,q))
118 hpad = h[-lh+nz_pre:]
119
120 offset = numpy.floor((l+nz_pre)/q)
121 nz_post = 0;
122 while numpy.ceil(((ls-1)*p + nz_pre + lh + nz_post )/q ) - offset < ly:
123 nz_post += 1
124 hpad = hpad[:lh + nz_pre + nz_post]
125
126 #filtering
127 xfilt = upfirdn(s, hpad, p, q)
128
129 return xfilt[offset-1:offset-1+ly]
130
131
132 def upfirdn(s, h, p, q):
133 """Upsample signal s by p, apply FIR filter as specified by h, and
134 downsample by q. Using fftconvolve as opposed to lfilter as it does not seem
135 to do a full convolution operation (and its much faster than convolve).
136 """
137 return downsample(signal.fftconvolve(h, upsample(s, p)), q)
138
139 # def main():
140 # """Show simple use cases for functionality provided by this module. Each
141 # example below attempts to mimic the examples provided by mathworks MATLAB
142 # documentation, http://www.mathworks.com/help/toolbox/signal/
143 # """
144 # import pylab
145 # argv = sys.argv
146 # if len(argv) != 1:
147 # print >>sys.stderr, 'usage: python -m pim.sp.multirate'
148 # sys.exit(2)
149 #
150 # #Downsample
151 # x = numpy.arange(1, 11)
152 # print 'Down Sampling %s by 3' % x
153 # print downsample(x, 3)
154 # print 'Down Sampling %s by 3 with phase offset 2' % x
155 # print downsample(x, 3, phase=2)
156 #
157 # #Upsample
158 # x = numpy.arange(1, 5)
159 # print 'Up Sampling %s by 3' % x
160 # print upsample(x, 3)
161 # print 'Up Sampling %s by 3 with phase offset 2' % x
162 # print upsample(x, 3, 2)
163 #
164 # #Decimate
165 # t = numpy.arange(0, 1, 0.00025)
166 # x = numpy.sin(2*numpy.pi*30*t) + numpy.sin(2*numpy.pi*60*t)
167 # y = decimate(x,4)
168 # pylab.figure()
169 # pylab.subplot(2, 1, 1)
170 # pylab.title('Original Signal')
171 # pylab.stem(numpy.arange(len(x[0:120])), x[0:120])
172 # pylab.subplot(2, 1, 2)
173 # pylab.title('Decimated Signal')
174 # pylab.stem(numpy.arange(len(y[0:30])), y[0:30])
175 #
176 # #Interp
177 # t = numpy.arange(0, 1, 0.001)
178 # x = numpy.sin(2*numpy.pi*30*t) + numpy.sin(2*numpy.pi*60*t)
179 # y = interp(x,4)
180 # pylab.figure()
181 # pylab.subplot(2, 1, 1)
182 # pylab.title('Original Signal')
183 # pylab.stem(numpy.arange(len(x[0:30])), x[0:30])
184 # pylab.subplot(2, 1, 2)
185 # pylab.title('Interpolated Signal')
186 # pylab.stem(numpy.arange(len(y[0:120])), y[0:120])
187 #
188 # #upfirdn
189 # L = 147.0
190 # M = 160.0
191 # N = 24.0*L
192 # h = signal.firwin(N-1, 1/M, window=('kaiser', 7.8562))
193 # h = L*h
194 # Fs = 48000.0
195 # n = numpy.arange(0, 10239)
196 # x = numpy.sin(2*numpy.pi*1000/Fs*n)
197 # y = upfirdn(x, h, L, M)
198 # pylab.figure()
199 # pylab.stem(n[1:49]/Fs, x[1:49])
200 # pylab.stem(n[1:45]/(Fs*L/M), y[13:57], 'r', markerfmt='ro',)
201 # pylab.xlabel('Time (sec)')
202 # pylab.ylabel('Signal value')
203 #
204 # #resample
205 # fs1 = 10.0
206 # t1 = numpy.arange(0, 1 + 1.0/fs1, 1.0/fs1)
207 # x = t1
208 # y = resample(x, 3, 2)
209 # t2 = numpy.arange(0,(len(y)))*2.0/(3.0*fs1)
210 # pylab.figure()
211 # pylab.plot(t1, x, '*')
212 # pylab.plot(t2, y, 'o')
213 # pylab.plot(numpy.arange(-0.5,1.5, 0.01), numpy.arange(-0.5,1.5, 0.01), ':')
214 # pylab.legend(('original','resampled'))
215 # pylab.xlabel('Time')
216 #
217 # x = numpy.hstack([numpy.arange(1,11), numpy.arange(9,0,-1)])
218 # y = resample(x,3,2)
219 # pylab.figure()
220 # pylab.subplot(2, 1, 1)
221 # pylab.title('Edge Effects Not Noticeable')
222 # pylab.plot(numpy.arange(19)+1, x, '*')
223 # pylab.plot(numpy.arange(29)*2/3.0 + 1, y, 'o')
224 # pylab.legend(('original', 'resampled'))
225 # x = numpy.hstack([numpy.arange(10, 0, -1), numpy.arange(2,11)])
226 # y = resample(x,3,2)
227 # pylab.subplot(2, 1, 2)
228 # pylab.plot(numpy.arange(19)+1, x, '*')
229 # pylab.plot(numpy.arange(29)*2/3.0 + 1, y, 'o')
230 # pylab.title('Edge Effects Very Noticeable')
231 # pylab.legend(('original', 'resampled'))
232 #
233 # pylab.show()
234 # return 0
235 #
236 # if __name__ == '__main__':
237 # sys.exit(main())
238 No newline at end of file
Show file before | Show file after | Hide comments
@@ -16,10 +16,10
16 J = 4;
16 J = 4;
17 [Faf, Fsf] = FSfarras;
17 [Faf, Fsf] = FSfarras;
18 [af, sf] = dualfilt1;
18 [af, sf] = dualfilt1;
19 # #
19
20 # # # compute transform of zero vector
20 # compute transform of zero vector
21 # # x = zeros(1,N);
21 x = np.zeros(shape=(1,N));
22 # # w = dualtree(x, J, Faf, af);
22 #w = dualtree(x, J, Faf, af);
23 # #
23 # #
24 # # # Uses both real and imaginary wavelets
24 # # # Uses both real and imaginary wavelets
25 # # for i in range (1, J+1):
25 # # for i in range (1, J+1):
Show file before | Show file after | Hide comments
@@ -30,8 +30,10
30 u = Q*c; # minimum 2-norm solution
30 u = Q*c; # minimum 2-norm solution
31 I = sps.eye(M);
31 I = sps.eye(M);
32
32
33 #---------- not needed, defined above--------------
33 # Spacing of floating point numbers
34 # Spacing of floating point numbers
34 eps = np.spacing(1)
35 #eps = np.spacing(1)
36 #--------------------------------------------------
35
37
36 # Loop until damping parameter is small enough
38 # Loop until damping parameter is small enough
37 while (eps > 1e-7):
39 while (eps > 1e-7):
Show file before | Show file after | Hide comments
@@ -1,10 +1,11
1 '''
1 '''
2 Created on May 27, 2014
2 Created on May 30, 2014
3
3
4 @author: Yolian Amaro
4 @author: Yolian Amaro
5 '''
5 '''
6
6
7 from irls_dn import *
7 from irls_dn import *
8 from scipy.optimize import fsolve
8
9
9 def irls_dn2(A,b,p,G):
10 def irls_dn2(A,b,p,G):
10
11
@@ -17,34 +18,38
17 # Start with a large lambda1, and do a line search until fidelity <= G.
18 # Start with a large lambda1, and do a line search until fidelity <= G.
18 # (Inversions with large lambda1 are really fast anyway).
19 # (Inversions with large lambda1 are really fast anyway).
19
20
20 # Then spin up fzero to localize the root even better
21 # Then spin up fsolve to localize the root even better
21
22
22 # Line Search
23 # Line Search
23
24
24 alpha = 2; # Line search parameter
25 alpha = 2; # Line search parameter
25 lambda1 = 1e5; # What's a reasonable but safe initial guess?
26 lambda1 = 1e5; # What's a reasonable but safe initial guess?
26 u = irls_dn(A,b,p,lambda1);
27 u = irls_dn(A,b,p,lambda1);
27 # fid = np.norm(A*u-b)^2;
28 fid = np.norm(A*u-b)**2;
28 #
29
29 # print '----------------------------------\n';
30 print '----------------------------------\n';
30 #
31
31 # while (fid >= G)
32 while (fid >= G):
32 # lambda1 = lambda1 / alpha; # Balance between speed and accuracy
33 lambda1 = lambda1 / alpha; # Balance between speed and accuracy
33 # u = irls_dn(A,b,p,lambda1);
34 u = irls_dn(A,b,p,lambda1);
34 # fid = np.norm(A*u-b)^2;
35 fid = np.norm(A*u-b)**2;
35 # print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid);
36 print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid;
36 #
37
37 # # Refinement using fzero
38 # Refinement using fzero
38 # lambda10 = [lambda1 lambda1*alpha]; # interval with zero-crossing
39 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
39 # f = @(lambda1) np.norm(A*irls_dn(A,b,p,lambda1) - b)^2 - G;
40
40 # opts = optimset('fzero');
41 f = lambda lambda1: np.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
42
43
44 # opts = optimset('fzero');
41 # # opts.Display = 'iter';
45 # # opts.Display = 'iter';
42 # opts.Display = 'none';
46 # opts.Display = 'none';
43 # opts.TolX = 0.01*lambda1;
47 # opts.TolX = 0.01*lambda1;
44 # lambda1 = fzero(f,lambda10,opts);
45 #
46 # u = irls_dn(A,b,p,lambda1);
47 #
48 #
49 # return u;
50
48
49 lambda1 = fsolve(f,lambda0); # FALTA OPTIMIZE ESTO
50
51 u = irls_dn(A,b,p,lambda1);
52
53
54 return u;
55
General Comments 0
You need to be logged in to leave comments. Login now