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