«

»

Jul 09

Nyquist Signal Expansion with Python

resampler.py

import numpy as np
import scipy.signal as signal

def expander(x_in, L):
xout = np.zeros(x_in.size*L)
xout[::L] = x_in
return xout

def resample(x_in, L, M):
#assume that fs_in = 1 for this function
x_e = expander(x_in, L)
h = L*signal.remez(20*L+1, [0, 0.35/float(L), 0.5/float(L), 0.5], [1, 0],[1, 1])
x_ef = signal.lfilter(h,1,x_e)
delay = (len(h)-1)/2.0
return [x_ef[0::M], delay]

if __name__ == "__main__":
import matplotlib.pyplot as plt
import sys
if len(sys.argv) > 2:
L = int(sys.argv[1])
M = int(sys.argv[2])
else:
L=3
M=2
fs_in = 6
Ts_in = 1/float(fs_in)
t_in = np.r_[0:5:Ts_in]
x_in = np.sin(2*np.pi*t_in) + 2*np.cos(4*np.pi*t_in)
t_cont = np.r_[0:5:0.01]
x_cont = np.sin(2*np.pi*t_cont) + 2*np.cos(4*np.pi*t_cont)
[x_out, delay] = resample (x_in, L, M)
Ts_out=Ts_in*M/float(L)
t_out = np.r_[0:len(x_out)]*Ts_out

delay_time = Ts_in/float(L)*delay
x_delayed = np.sin(2*np.pi*(t_cont-delay_time))+2*np.cos(4*np.pi*(t_cont-delay_time))


f2 = plt.figure()
ax2 = f2.add_subplot(111)
ax2.plot(t_cont,x_delayed,alpha=0.2)
#ax2.plot(t_out, x_out, 'g')
ax2.stem(t_out, x_out, linefmt='g--', markerfmt='go', basefmt='k--')
ax2.set_title('Resampled Data L=%d, M=%d' % (L,M))
ax2.set_xlabel('[s]')
ax2.text(delay_time,2.5, 'Filter Delay = %1.2fs' % delay_time, clip_on=True, multialignment='left', backgroundcolor='white')
ax2.grid()
f2.savefig("resampler test L%d M%d.png" % (L, M), pad_inches=0.1, dpi=72, bbox_inches='tight')
plt.show()

 

Permanent link to this article: http://blog.curioussystem.com/2013/07/nyquist-signal-expansion-with-python/

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>