signature:
def page_posts ( self , page_id , after = '' , post_type = "posts" , include_hidden = False , fields = None , ** params ) :
    """: param page _ id :
 : param after :
 : param post _ type : Can be ' posts ' , ' feed ' , ' tagged ' , ' promotable _ posts '
 : param include _ hidden :
 : param fields :
 : param params :
 : return :"""

implementation:
if fields :
    fields = "," . join ( fields )
parameters = { "access_token" : self . key , "after" : after , "fields" : fields , "include_hidden" : include_hidden }
parameters = self . merge_params ( parameters , params )
return self . api_call ( '%s/%s' % ( page_id , post_type ) , parameters )


signature:
def post_comments ( self , post_id , after = '' , order = "chronological" , filter = "stream" , fields = None , ** params ) :
    """: param post _ id :
 : param after :
 : param order : Can be ' ranked ' , ' chronological ' , ' reverse _ chronological '
 : param filter : Can be ' stream ' , ' toplevel '
 : param fields : Can be ' id ' , ' application ' , ' attachment ' , ' can _ comment ' ,
 ' can _ remove ' , ' can _ hide ' , ' can _ like ' , ' can _ reply _ privately ' , ' comments ' ,
 ' comment _ count ' , ' created _ time ' , ' from ' , ' likes ' , ' like _ count ' ,
 ' live _ broadcast _ timestamp ' , ' message ' , ' message _ tags ' , ' object ' ,
 ' parent ' , ' private _ reply _ conversation ' , ' user _ likes '
 : param params :
 : return :"""

implementation:
if fields :
    fields = "," . join ( fields )
parameters = { "access_token" : self . key , "after" : after , "order" : order , "fields" : fields , "filter" : filter }
parameters = self . merge_params ( parameters , params )
return self . api_call ( '%s/comments' % post_id , parameters )


signature:
def flatten ( dictionary , parent_key = False , separator = '.' ) :
    """Turn a nested dictionary into a flattened dictionary
 : param dictionary : The dictionary to flatten
 : param parent _ key : The string to prepend to dictionary ' s keys
 : param separator : The string used to separate flattened keys
 : return : A flattened dictionary"""

implementation:
items = [ ]
for key , value in dictionary . items ( ) :
    new_key = str ( parent_key ) + separator + key if parent_key else key
    if isinstance ( value , collections . MutableMapping ) :
        items . extend ( flatten ( value , new_key , separator ) . items ( ) )
    elif isinstance ( value , list ) :
        for k , v in enumerate ( value ) :
            items . extend ( flatten ( { str ( k ) : v } , new_key ) . items ( ) )
    else :
        items . append ( ( new_key , value ) )
return dict ( items )


signature:
def fill_gaps ( list_dicts ) :
    """Fill gaps in a list of dictionaries . Add empty keys to dictionaries in
 the list that don ' t contain other entries ' keys
 : param list _ dicts : A list of dictionaries
 : return : A list of field names , a list of dictionaries with identical keys"""

implementation:
field_names = [ ]
# ! = set bc . preserving order is better for output
for datum in list_dicts :
    for key in datum . keys ( ) :
        if key not in field_names :
            field_names . append ( key )
for datum in list_dicts :
    for key in field_names :
        if key not in datum :
            datum [ key ] = ''
return list ( field_names ) , list_dicts


signature:
def to_csv ( data , field_names = None , filename = 'data.csv' , overwrite = True , write_headers = True , append = False , flat = True , primary_fields = None , sort_fields = True ) :
    """DEPRECATED Write a list of dicts to a csv file
 : param data : List of dicts
 : param field _ names : The list column names
 : param filename : The name of the file
 : param overwrite : Overwrite the file if exists
 : param write _ headers : Write the headers to the csv file
 : param append : Write new rows if the file exists
 : param flat : Flatten the dictionary before saving
 : param primary _ fields : The first columns of the csv file
 : param sort _ fields : Sort the field names alphabetically
 : return : None"""

implementation:
# Don ' t overwrite if not specified
if not overwrite and path . isfile ( filename ) :
    raise FileExistsError ( 'The file already exists' )
# Replace file if append not specified
write_type = 'w' if not append else 'a'
# Flatten if flat is specified , or there are no predefined field names
if flat or not field_names :
    data = [ flatten ( datum ) for datum in data ]
# Fill in gaps between dicts with empty string
if not field_names :
    field_names , data = fill_gaps ( data )
# Sort fields if specified
if sort_fields :
    field_names . sort ( )
# If there are primary fields , move the field names to the front and sort
# based on first field
if primary_fields :
    for key in primary_fields [ : : - 1 ] :
        field_names . insert ( 0 , field_names . pop ( field_names . index ( key ) ) )
    data = sorted ( data , key = lambda k : k [ field_names [ 0 ] ] , reverse = True )
# Write the file
with open ( filename , write_type , encoding = 'utf-8' ) as f :
    writer = csv . DictWriter ( f , fieldnames = field_names , lineterminator = '\n' )
    if not append or write_headers :
        writer . writeheader ( )
    # Write rows containing fields in field names
    for datum in data :
        for key in list ( datum . keys ( ) ) :
            if key not in field_names :
                del datum [ key ]
            elif type ( datum [ key ] ) is str :
                datum [ key ] = datum [ key ] . strip ( )
            datum [ key ] = str ( datum [ key ] )
        writer . writerow ( datum )


signature:
def to_json ( data , filename = 'data.json' , indent = 4 ) :
    """Write an object to a json file
 : param data : The object
 : param filename : The name of the file
 : param indent : The indentation of the file
 : return : None"""

implementation:
with open ( filename , 'w' ) as f :
    f . write ( json . dumps ( data , indent = indent ) )


signature:
def save_file ( filename , source , folder = "Downloads" ) :
    """Download and save a file at path
 : param filename : The name of the file
 : param source : The location of the resource online
 : param folder : The directory the file will be saved in
 : return : None"""

implementation:
r = requests . get ( source , stream = True )
if r . status_code == 200 :
    if not path . isdir ( folder ) :
        makedirs ( folder , exist_ok = True )
    with open ( "%s/%s" % ( folder , filename ) , 'wb' ) as f :
        for chunk in r :
            f . write ( chunk )


signature:
def convert_frames_to_video ( tar_file_path , output_path = "output.mp4" , framerate = 60 , overwrite = False ) :
    """Try to convert a tar file containing a sequence of frames saved by the
 meshcat viewer into a single video file .
 This relies on having ` ffmpeg ` installed on your system ."""

implementation:
output_path = os . path . abspath ( output_path )
if os . path . isfile ( output_path ) and not overwrite :
    raise ValueError ( "The output path {:s} already exists. To overwrite that file, you can pass overwrite=True to this function." . format ( output_path ) )
with tempfile . TemporaryDirectory ( ) as tmp_dir :
    with tarfile . open ( tar_file_path ) as tar :
        tar . extractall ( tmp_dir )
    args = [ "ffmpeg" , "-r" , str ( framerate ) , "-i" , r"%07d.png" , "-vcodec" , "libx264" , "-preset" , "slow" , "-crf" , "18" ]
    if overwrite :
        args . append ( "-y" )
    args . append ( output_path )
    try :
        subprocess . check_call ( args , cwd = tmp_dir )
    except subprocess . CalledProcessError as e :
        print ( """
Could not call `ffmpeg` to convert your frames into a video.
If you want to convert the frames manually, you can extract the
.tar archive into a directory, cd to that directory, and run:
ffmpeg -r 60 -i %07d.png \\\n\t -vcodec libx264 \\\n\t -preset slow \\\n\t -crf 18 \\\n\t output.mp4
                """ )
        raise
print ( "Saved output as {:s}" . format ( output_path ) )
return output_path


signature:
def toJSON ( self ) :
    """Get a json dict of the attributes of this object ."""

implementation: return { "id" : self . id , "compile" : self . compile , "position" : self . position , "version" : self . version }


signature:
def from_string ( cls , content , position = 1 , file_id = None ) :
    """Convenience method to create a file from a string .
 This file object ' s metadata will have the id ' inlined _ input ' .
 Inputs
 content - - the content of the file ( a string ) .
 position - - ( default 1 ) rank among all files of the model while parsing
 see FileMetadata
 file _ id - - ( default ' inlined _ input ' ) the file _ id that will be used by
 kappa ."""

implementation:
if file_id is None :
    file_id = 'inlined_input'
return cls ( FileMetadata ( file_id , position ) , content )


signature:
def from_file ( cls , fpath , position = 1 , file_id = None ) :
    """Convience method to create a kappa file object from a file on disk
 Inputs
 fpath - - path to the file on disk
 position - - ( default 1 ) rank among all files of the model while parsing
 see FileMetadata
 file _ id - - ( default = fpath ) the file _ id that will be used by kappa ."""

implementation:
if file_id is None :
    file_id = fpath
with open ( fpath ) as f :
    code = f . read ( )
    file_content = str ( code )
    file_metadata = FileMetadata ( file_id , position )
    return cls ( file_metadata , file_content )


signature:
def _fix_docs ( this_abc , child_class ) :
    """Make api method docs inheritted .
 Specifically , insepect . getdoc will return values inheritted from this
 abc for standardized api methods ."""

implementation:
# After python 3.5 , this is basically handled automatically
if sys . version_info >= ( 3 , 5 ) :
    return child_class
if not issubclass ( child_class , this_abc ) :
    raise KappaError ( 'Cannot fix docs of class that is not decendent.' )
# This method is modified from solution given in
# https : / / stackoverflow . com / a / 8101598/8863865
for name , child_func in vars ( child_class ) . items ( ) :
    if callable ( child_func ) and not child_func . __doc__ :
        if name in this_abc . __abstractmethods__ :
            parent_func = getattr ( this_abc , name )
            child_func . __doc__ = parent_func . __doc__
return child_class


signature:
def add_model_string ( self , model_str , position = 1 , file_id = None ) :
    """Add a kappa model given in a string to the project ."""

implementation:
if file_id is None :
    file_id = self . make_unique_id ( 'inlined_input' )
ret_data = self . file_create ( File . from_string ( model_str , position , file_id ) )
return ret_data


signature:
def add_model_file ( self , model_fpath , position = 1 , file_id = None ) :
    """Add a kappa model from a file at given path to the project ."""

implementation:
if file_id is None :
    file_id = self . make_unique_id ( 'file_input' )
ret_data = self . file_create ( File . from_file ( model_fpath , position , file_id ) )
return ret_data


signature:
def set_default_sim_param ( self , * args , ** kwargs ) :
    """Set the simulation default simulation parameters .
 You can pass one of two things in as input :
 - a kappa _ common . SimulationParameter instance
 - the arguments and keyword argument to create such an instance .
 The parameters you specify will be used by default in simulations run
 by this client ."""

implementation:
if len ( args ) is 1 and isinstance ( args [ 0 ] , SimulationParameter ) :
    self . __default_param = args [ 0 ]
else :
    self . __default_param = SimulationParameter ( * args , ** kwargs )
return


signature:
def get_is_sim_running ( self ) :
    """Check if the current simulation is running ."""

implementation:
sim_info = self . simulation_info ( )
try :
    progress_info = sim_info [ 'simulation_info_progress' ]
    ret = progress_info [ 'simulation_progress_is_running' ]
except KeyError : # Simulation has not been created .
    ret = False
return ret


signature:
def wait_for_simulation_stop ( self , timeout = None ) :
    """Block until the simulation is done or timeout seconds exceeded .
 If the simulation stops before timeout , siminfo is returned ."""

implementation:
start = datetime . now ( )
while self . get_is_sim_running ( ) :
    sleep ( 0.5 )
    if timeout is not None :
        if ( datetime . now ( ) - start ) . seconds >= timeout :
            ret = None
            break
else :
    ret = self . simulation_info ( )
return ret


signature:
def available_devices ( ) :
    """Display available input and output audio devices along with their
 port indices .
 : return : Dictionary whose keys are the device index , the number of inputs and outputs , and their names .
 : rtype : dict"""

implementation:
devices = { }
pA = pyaudio . PyAudio ( )
device_string = str ( )
for k in range ( pA . get_device_count ( ) ) :
    dev = pA . get_device_info_by_index ( k )
    devices [ k ] = { 'name' : dev [ 'name' ] , 'inputs' : dev [ 'maxInputChannels' ] , 'outputs' : dev [ 'maxOutputChannels' ] }
    device_string += 'Index %d device name = %s, inputs = %d, outputs = %d\n' % ( k , dev [ 'name' ] , dev [ 'maxInputChannels' ] , dev [ 'maxOutputChannels' ] )
logger . debug ( device_string )
return devices


signature:
def in_out_check ( self ) :
    """Checks the input and output to see if they are valid"""

implementation:
devices = available_devices ( )
if not self . in_idx in devices :
    raise OSError ( "Input device is unavailable" )
in_check = devices [ self . in_idx ]
if not self . out_idx in devices :
    raise OSError ( "Output device is unavailable" )
out_check = devices [ self . out_idx ]
if ( ( in_check [ 'inputs' ] == 0 ) and ( out_check [ 'outputs' ] == 0 ) ) :
    raise StandardError ( 'Invalid input and output devices' )
elif ( in_check [ 'inputs' ] == 0 ) :
    raise ValueError ( 'Selected input device has no inputs' )
elif ( out_check [ 'outputs' ] == 0 ) :
    raise ValueError ( 'Selected output device has no outputs' )
return True


signature:
def interactive_stream ( self , Tsec = 2 , numChan = 1 ) :
    """Stream audio with start and stop radio buttons
 Interactive stream is designed for streaming audio through this object using
 a callback function . This stream is threaded , so it can be used with ipywidgets .
 Click on the " Start Streaming " button to start streaming and click on " Stop Streaming "
 button to stop streaming .
 Parameters
 Tsec : stream time in seconds if Tsec > 0 . If Tsec = 0 , then stream goes to infinite
 mode . When in infinite mode , the " Stop Streaming " radio button or Tsec . stop ( ) can be
 used to stop the stream .
 numChan : number of channels . Use 1 for mono and 2 for stereo ."""

implementation:
self . Tsec = Tsec
self . numChan = numChan
self . interactiveFG = 1
self . play = interactive ( self . interaction , Stream = ToggleButtons ( options = [ 'Start Streaming' , 'Stop Streaming' ] , description = ' ' , value = 'Stop Streaming' ) )
display ( self . play )


signature:
def thread_stream ( self , Tsec = 2 , numChan = 1 ) :
    """Stream audio in a thread using callback . The stream is threaded , so widgets can be
 used simultaneously during stream .
 Parameters
 Tsec : stream time in seconds if Tsec > 0 . If Tsec = 0 , then stream goes to infinite
 mode . When in infinite mode , Tsec . stop ( ) can be used to stop the stream .
 numChan : number of channels . Use 1 for mono and 2 for stereo ."""

implementation:
def stream_thread ( time , channel ) :
    self . stream ( Tsec = time , numChan = channel )
# Thread the streaming function
t = Thread ( target = stream_thread , args = ( Tsec , numChan , ) )
# Start the stream
t . start ( )


signature:
def stream ( self , Tsec = 2 , numChan = 1 ) :
    """Stream audio using callback
 Parameters
 Tsec : stream time in seconds if Tsec > 0 . If Tsec = 0 , then stream goes to infinite
 mode . When in infinite mode , Tsec . stop ( ) can be used to stop the stream .
 numChan : number of channels . Use 1 for mono and 2 for stereo ."""

implementation:
self . Tsec = Tsec
self . numChan = numChan
self . N_samples = int ( self . fs * Tsec )
self . data_capture = [ ]
self . data_capture_left = [ ]
self . data_capture_right = [ ]
self . capture_sample_count = 0
self . DSP_tic = [ ]
self . DSP_toc = [ ]
self . start_time = time . time ( )
self . stop_stream = False
# open stream using callback ( 3)
stream = self . p . open ( format = pyaudio . paInt16 , channels = numChan , rate = self . fs , input = True , output = True , input_device_index = self . in_idx , output_device_index = self . out_idx , frames_per_buffer = self . frame_length , stream_callback = self . stream_callback )
# start the stream ( 4)
stream . start_stream ( )
# infinite mode
if ( Tsec == 0 ) :
    while stream . is_active ( ) :
        if self . stop_stream :
            stream . stop_stream ( )
        time . sleep ( self . sleep_time )
else : # wait for stream to finish ( 5)
    while stream . is_active ( ) :
        if self . capture_sample_count >= self . N_samples :
            stream . stop_stream ( )
        if self . stop_stream :
            stream . stop_stream ( )
        time . sleep ( self . sleep_time )
# stop stream ( 6)
stream . stop_stream ( )
stream . close ( )
# close PyAudio ( 7)
self . p . terminate ( )
self . stream_data = True
# print ( ' Audio input / output streaming session complete ! ' )
if ( self . interactiveFG ) : # Move radio button back to ' Stop Streaming '
    self . play . children [ 0 ] . value = 'Stop Streaming'
else :
    if ( self . print_when_done == 1 ) :
        print ( 'Completed' )


signature:
def DSP_capture_add_samples ( self , new_data ) :
    """Append new samples to the data _ capture array and increment the sample counter
 If length reaches Tcapture , then the newest samples will be kept . If Tcapture = 0
 then new values are not appended to the data _ capture array ."""

implementation:
self . capture_sample_count += len ( new_data )
if self . Tcapture > 0 :
    self . data_capture = np . hstack ( ( self . data_capture , new_data ) )
    if ( self . Tcapture > 0 ) and ( len ( self . data_capture ) > self . Ncapture ) :
        self . data_capture = self . data_capture [ - self . Ncapture : ]


signature:
def DSP_capture_add_samples_stereo ( self , new_data_left , new_data_right ) :
    """Append new samples to the data _ capture _ left array and the data _ capture _ right
 array and increment the sample counter . If length reaches Tcapture , then the
 newest samples will be kept . If Tcapture = 0 then new values are not appended
 to the data _ capture array ."""

implementation:
self . capture_sample_count = self . capture_sample_count + len ( new_data_left ) + len ( new_data_right )
if self . Tcapture > 0 :
    self . data_capture_left = np . hstack ( ( self . data_capture_left , new_data_left ) )
    self . data_capture_right = np . hstack ( ( self . data_capture_right , new_data_right ) )
    if ( len ( self . data_capture_left ) > self . Ncapture ) :
        self . data_capture_left = self . data_capture_left [ - self . Ncapture : ]
    if ( len ( self . data_capture_right ) > self . Ncapture ) :
        self . data_capture_right = self . data_capture_right [ - self . Ncapture : ]


signature:
def DSP_callback_tic ( self ) :
    """Add new tic time to the DSP _ tic list . Will not be called if
 Tcapture = 0."""

implementation:
if self . Tcapture > 0 :
    self . DSP_tic . append ( time . time ( ) - self . start_time )


signature:
def DSP_callback_toc ( self ) :
    """Add new toc time to the DSP _ toc list . Will not be called if
 Tcapture = 0."""

implementation:
if self . Tcapture > 0 :
    self . DSP_toc . append ( time . time ( ) - self . start_time )


signature:
def stream_stats ( self ) :
    """Display basic statistics of callback execution : ideal period
 between callbacks , average measured period between callbacks ,
 and average time spent in the callback ."""

implementation:
Tp = self . frame_length / float ( self . fs ) * 1000
print ( 'Delay (latency) in Entering the Callback the First Time = %6.2f (ms)' % ( self . DSP_tic [ 0 ] * 1000 , ) )
print ( 'Ideal Callback period = %1.2f (ms)' % Tp )
Tmp_mean = np . mean ( np . diff ( np . array ( self . DSP_tic ) ) [ 1 : ] * 1000 )
print ( 'Average Callback Period = %1.2f (ms)' % Tmp_mean )
Tprocess_mean = np . mean ( np . array ( self . DSP_toc ) - np . array ( self . DSP_tic ) ) * 1000
print ( 'Average Callback process time = %1.2f (ms)' % Tprocess_mean )


signature:
def cb_active_plot ( self , start_ms , stop_ms , line_color = 'b' ) :
    """Plot timing information of time spent in the callback . This is similar
 to what a logic analyzer provides when probing an interrupt .
 cb _ active _ plot ( start _ ms , stop _ ms , line _ color = ' b ' )"""

implementation:
# Find bounding k values that contain the [ start _ ms , stop _ ms ]
k_min_idx = np . nonzero ( np . ravel ( np . array ( self . DSP_tic ) * 1000 < start_ms ) ) [ 0 ]
if len ( k_min_idx ) < 1 :
    k_min = 0
else :
    k_min = k_min_idx [ - 1 ]
k_max_idx = np . nonzero ( np . ravel ( np . array ( self . DSP_tic ) * 1000 > stop_ms ) ) [ 0 ]
if len ( k_min_idx ) < 1 :
    k_max = len ( self . DSP_tic )
else :
    k_max = k_max_idx [ 0 ]
for k in range ( k_min , k_max ) :
    if k == 0 :
        plt . plot ( [ 0 , self . DSP_tic [ k ] * 1000 , self . DSP_tic [ k ] * 1000 , self . DSP_toc [ k ] * 1000 , self . DSP_toc [ k ] * 1000 ] , [ 0 , 0 , 1 , 1 , 0 ] , 'b' )
    else :
        plt . plot ( [ self . DSP_toc [ k - 1 ] * 1000 , self . DSP_tic [ k ] * 1000 , self . DSP_tic [ k ] * 1000 , self . DSP_toc [ k ] * 1000 , self . DSP_toc [ k ] * 1000 ] , [ 0 , 0 , 1 , 1 , 0 ] , 'b' )
plt . plot ( [ self . DSP_toc [ k_max - 1 ] * 1000 , stop_ms ] , [ 0 , 0 ] , 'b' )
plt . xlim ( [ start_ms , stop_ms ] )
plt . title ( r'Time Spent in the callback' )
plt . ylabel ( r'Timing' )
plt . xlabel ( r'Time (ms)' )
plt . grid ( ) ;


signature:
def get_LR ( self , in_data ) :
    """Splits incoming packed stereo data into separate left and right channels
 and returns an array of left samples and an array of right samples
 Parameters
 in _ data : input data from the streaming object in the callback function .
 Returns
 left _ in : array of incoming left channel samples
 right _ in : array of incoming right channel samples"""

implementation:
for i in range ( 0 , self . frame_length * 2 ) :
    if i % 2 :
        self . right_in [ ( int ) ( i / 2 ) ] = in_data [ i ]
    else :
        self . left_in [ ( int ) ( i / 2 ) ] = in_data [ i ]
return self . left_in , self . right_in


signature:
def pack_LR ( self , left_out , right_out ) :
    """Packs separate left and right channel data into one array to output
 and returns the output .
 Parameters
 left _ out : left channel array of samples going to output
 right _ out : right channel array of samples going to output
 Returns
 out : packed left and right channel array of samples"""

implementation:
for i in range ( 0 , self . frame_length * 2 ) :
    if i % 2 :
        self . out [ i ] = right_out [ ( int ) ( i / 2 ) ]
    else :
        self . out [ i ] = left_out [ ( int ) ( i / 2 ) ]
return self . out


signature:
def IIR_lpf ( f_pass , f_stop , Ripple_pass , Atten_stop , fs = 1.00 , ftype = 'butter' ) :
    """Design an IIR lowpass filter using scipy . signal . iirdesign .
 The filter order is determined based on
 f _ pass Hz , f _ stop Hz , and the desired stopband attenuation
 d _ stop in dB , all relative to a sampling rate of fs Hz .
 Parameters
 f _ pass : Passband critical frequency in Hz
 f _ stop : Stopband critical frequency in Hz
 Ripple _ pass : Filter gain in dB at f _ pass
 Atten _ stop : Filter attenuation in dB at f _ stop
 fs : Sampling rate in Hz
 ftype : Analog prototype from ' butter ' ' cheby1 ' , ' cheby2 ' ,
 ' ellip ' , and ' bessel '
 Returns
 b : ndarray of the numerator coefficients
 a : ndarray of the denominator coefficients
 sos : 2D ndarray of second - order section coefficients
 Notes
 Additionally a text string telling the user the filter order is
 written to the console , e . g . , IIR cheby1 order = 8.
 Examples
 > > > fs = 48000
 > > > f _ pass = 5000
 > > > f _ stop = 8000
 > > > b _ but , a _ but , sos _ but = IIR _ lpf ( f _ pass , f _ stop , 0.5,60 , fs , ' butter ' )
 > > > b _ cheb1 , a _ cheb1 , sos _ cheb1 = IIR _ lpf ( f _ pass , f _ stop , 0.5,60 , fs , ' cheby1 ' )
 > > > b _ cheb2 , a _ cheb2 , sos _ cheb2 = IIR _ lpf ( f _ pass , f _ stop , 0.5,60 , fs , ' cheby2 ' )
 > > > b _ elli , a _ elli , sos _ elli = IIR _ lpf ( f _ pass , f _ stop , 0.5,60 , fs , ' ellip ' )
 Mark Wickert October 2016"""

implementation:
b , a = signal . iirdesign ( 2 * float ( f_pass ) / fs , 2 * float ( f_stop ) / fs , Ripple_pass , Atten_stop , ftype = ftype , output = 'ba' )
sos = signal . iirdesign ( 2 * float ( f_pass ) / fs , 2 * float ( f_stop ) / fs , Ripple_pass , Atten_stop , ftype = ftype , output = 'sos' )
tag = 'IIR ' + ftype + ' order'
print ( '%s = %d.' % ( tag , len ( a ) - 1 ) )
return b , a , sos


signature:
def IIR_bsf ( f_pass1 , f_stop1 , f_stop2 , f_pass2 , Ripple_pass , Atten_stop , fs = 1.00 , ftype = 'butter' ) :
    """Design an IIR bandstop filter using scipy . signal . iirdesign .
 The filter order is determined based on
 f _ pass Hz , f _ stop Hz , and the desired stopband attenuation
 d _ stop in dB , all relative to a sampling rate of fs Hz .
 Mark Wickert October 2016"""

implementation:
b , a = signal . iirdesign ( [ 2 * float ( f_pass1 ) / fs , 2 * float ( f_pass2 ) / fs ] , [ 2 * float ( f_stop1 ) / fs , 2 * float ( f_stop2 ) / fs ] , Ripple_pass , Atten_stop , ftype = ftype , output = 'ba' )
sos = signal . iirdesign ( [ 2 * float ( f_pass1 ) / fs , 2 * float ( f_pass2 ) / fs ] , [ 2 * float ( f_stop1 ) / fs , 2 * float ( f_stop2 ) / fs ] , Ripple_pass , Atten_stop , ftype = ftype , output = 'sos' )
tag = 'IIR ' + ftype + ' order'
print ( '%s = %d.' % ( tag , len ( a ) - 1 ) )
return b , a , sos


signature:
def freqz_cas ( sos , w ) :
    """Cascade frequency response
 Mark Wickert October 2016"""

implementation:
Ns , Mcol = sos . shape
w , Hcas = signal . freqz ( sos [ 0 , : 3 ] , sos [ 0 , 3 : ] , w )
for k in range ( 1 , Ns ) :
    w , Htemp = signal . freqz ( sos [ k , : 3 ] , sos [ k , 3 : ] , w )
    Hcas *= Htemp
return w , Hcas


signature:
def unique_cpx_roots ( rlist , tol = 0.001 ) :
    """The average of the root values is used when multiplicity
 is greater than one .
 Mark Wickert October 2016"""

implementation:
uniq = [ rlist [ 0 ] ]
mult = [ 1 ]
for k in range ( 1 , len ( rlist ) ) :
    N_uniq = len ( uniq )
    for m in range ( N_uniq ) :
        if abs ( rlist [ k ] - uniq [ m ] ) <= tol :
            mult [ m ] += 1
            uniq [ m ] = ( uniq [ m ] * ( mult [ m ] - 1 ) + rlist [ k ] ) / float ( mult [ m ] )
            break
    uniq = np . hstack ( ( uniq , rlist [ k ] ) )
    mult = np . hstack ( ( mult , [ 1 ] ) )
return np . array ( uniq ) , np . array ( mult )


signature:
def sos_zplane ( sos , auto_scale = True , size = 2 , tol = 0.001 ) :
    """Create an z - plane pole - zero plot .
 Create an z - plane pole - zero plot using the numerator
 and denominator z - domain system function coefficient
 ndarrays b and a respectively . Assume descending powers of z .
 Parameters
 sos : ndarray of the sos coefficients
 auto _ scale : bool ( default True )
 size : plot radius maximum when scale = False
 Returns
 ( M , N ) : tuple of zero and pole counts + plot window
 Notes
 This function tries to identify repeated poles and zeros and will
 place the multiplicity number above and to the right of the pole or zero .
 The difficulty is setting the tolerance for this detection . Currently it
 is set at 1e - 3 via the function signal . unique _ roots .
 Examples
 > > > # Here the plot is generated using auto _ scale
 > > > sos _ zplane ( sos )
 > > > # Here the plot is generated using manual scaling
 > > > sos _ zplane ( sos , False , 1.5)"""

implementation:
Ns , Mcol = sos . shape
# Extract roots from sos num and den removing z = 0
# roots due to first - order sections
N_roots = [ ]
for k in range ( Ns ) :
    N_roots_tmp = np . roots ( sos [ k , : 3 ] )
    if N_roots_tmp [ 1 ] == 0. :
        N_roots = np . hstack ( ( N_roots , N_roots_tmp [ 0 ] ) )
    else :
        N_roots = np . hstack ( ( N_roots , N_roots_tmp ) )
D_roots = [ ]
for k in range ( Ns ) :
    D_roots_tmp = np . roots ( sos [ k , 3 : ] )
    if D_roots_tmp [ 1 ] == 0. :
        D_roots = np . hstack ( ( D_roots , D_roots_tmp [ 0 ] ) )
    else :
        D_roots = np . hstack ( ( D_roots , D_roots_tmp ) )
# Plot labels if multiplicity greater than 1
x_scale = 1.5 * size
y_scale = 1.5 * size
x_off = 0.02
y_off = 0.01
M = len ( N_roots )
N = len ( D_roots )
if auto_scale :
    if M > 0 and N > 0 :
        size = max ( np . max ( np . abs ( N_roots ) ) , np . max ( np . abs ( D_roots ) ) ) + .1
    elif M > 0 :
        size = max ( np . max ( np . abs ( N_roots ) ) , 1.0 ) + .1
    elif N > 0 :
        size = max ( 1.0 , np . max ( np . abs ( D_roots ) ) ) + .1
    else :
        size = 1.1
plt . figure ( figsize = ( 5 , 5 ) )
plt . axis ( 'equal' )
r = np . linspace ( 0 , 2 * np . pi , 200 )
plt . plot ( np . cos ( r ) , np . sin ( r ) , 'r--' )
plt . plot ( [ - size , size ] , [ 0 , 0 ] , 'k-.' )
plt . plot ( [ 0 , 0 ] , [ - size , size ] , 'k-.' )
if M > 0 : # N _ roots = np . roots ( b )
    N_uniq , N_mult = unique_cpx_roots ( N_roots , tol = tol )
    plt . plot ( np . real ( N_uniq ) , np . imag ( N_uniq ) , 'ko' , mfc = 'None' , ms = 8 )
    idx_N_mult = np . nonzero ( np . ravel ( N_mult > 1 ) ) [ 0 ]
    for k in range ( len ( idx_N_mult ) ) :
        x_loc = np . real ( N_uniq [ idx_N_mult [ k ] ] ) + x_off * x_scale
        y_loc = np . imag ( N_uniq [ idx_N_mult [ k ] ] ) + y_off * y_scale
        plt . text ( x_loc , y_loc , str ( N_mult [ idx_N_mult [ k ] ] ) , ha = 'center' , va = 'bottom' , fontsize = 10 )
if N > 0 : # D _ roots = np . roots ( a )
    D_uniq , D_mult = unique_cpx_roots ( D_roots , tol = tol )
    plt . plot ( np . real ( D_uniq ) , np . imag ( D_uniq ) , 'kx' , ms = 8 )
    idx_D_mult = np . nonzero ( np . ravel ( D_mult > 1 ) ) [ 0 ]
    for k in range ( len ( idx_D_mult ) ) :
        x_loc = np . real ( D_uniq [ idx_D_mult [ k ] ] ) + x_off * x_scale
        y_loc = np . imag ( D_uniq [ idx_D_mult [ k ] ] ) + y_off * y_scale
        plt . text ( x_loc , y_loc , str ( D_mult [ idx_D_mult [ k ] ] ) , ha = 'center' , va = 'bottom' , fontsize = 10 )
if M - N < 0 :
    plt . plot ( 0.0 , 0.0 , 'bo' , mfc = 'None' , ms = 8 )
elif M - N > 0 :
    plt . plot ( 0.0 , 0.0 , 'kx' , ms = 8 )
if abs ( M - N ) > 1 :
    plt . text ( x_off * x_scale , y_off * y_scale , str ( abs ( M - N ) ) , ha = 'center' , va = 'bottom' , fontsize = 10 )
plt . xlabel ( 'Real Part' )
plt . ylabel ( 'Imaginary Part' )
plt . title ( 'Pole-Zero Plot' )
# plt . grid ( )
plt . axis ( [ - size , size , - size , size ] )
return M , N


signature:
def firwin_bpf ( N_taps , f1 , f2 , fs = 1.0 , pass_zero = False ) :
    """Design a windowed FIR bandpass filter in terms of passband
 critical frequencies f1 < f2 in Hz relative to sampling rate
 fs in Hz . The number of taps must be provided .
 Mark Wickert October 2016"""

implementation: return signal . firwin ( N_taps , 2 * ( f1 , f2 ) / fs , pass_zero = pass_zero )


signature:
def firwin_kaiser_lpf ( f_pass , f_stop , d_stop , fs = 1.0 , N_bump = 0 ) :
    """Design an FIR lowpass filter using the sinc ( ) kernel and
 a Kaiser window . The filter order is determined based on
 f _ pass Hz , f _ stop Hz , and the desired stopband attenuation
 d _ stop in dB , all relative to a sampling rate of fs Hz .
 Note : the passband ripple cannot be set independent of the
 stopband attenuation .
 Mark Wickert October 2016"""

implementation:
wc = 2 * np . pi * ( f_pass + f_stop ) / 2 / fs
delta_w = 2 * np . pi * ( f_stop - f_pass ) / fs
# Find the filter order
M = np . ceil ( ( d_stop - 8 ) / ( 2.285 * delta_w ) )
# Adjust filter order up or down as needed
M += N_bump
N_taps = M + 1
# Obtain the Kaiser window
beta = signal . kaiser_beta ( d_stop )
w_k = signal . kaiser ( N_taps , beta )
n = np . arange ( N_taps )
b_k = wc / np . pi * np . sinc ( wc / np . pi * ( n - M / 2 ) ) * w_k
b_k /= np . sum ( b_k )
print ( 'Kaiser Win filter taps = %d.' % N_taps )
return b_k


signature:
def firwin_kaiser_bpf ( f_stop1 , f_pass1 , f_pass2 , f_stop2 , d_stop , fs = 1.0 , N_bump = 0 ) :
    """Design an FIR bandpass filter using the sinc ( ) kernel and
 a Kaiser window . The filter order is determined based on
 f _ stop1 Hz , f _ pass1 Hz , f _ pass2 Hz , f _ stop2 Hz , and the
 desired stopband attenuation d _ stop in dB for both stopbands ,
 all relative to a sampling rate of fs Hz .
 Note : the passband ripple cannot be set independent of the
 stopband attenuation .
 Mark Wickert October 2016"""

implementation:
# Design BPF starting from simple LPF equivalent
# The upper and lower stopbands are assumed to have
# the same attenuation level . The LPF equivalent critical
# frequencies :
f_pass = ( f_pass2 - f_pass1 ) / 2
f_stop = ( f_stop2 - f_stop1 ) / 2
# Continue to design equivalent LPF
wc = 2 * np . pi * ( f_pass + f_stop ) / 2 / fs
delta_w = 2 * np . pi * ( f_stop - f_pass ) / fs
# Find the filter order
M = np . ceil ( ( d_stop - 8 ) / ( 2.285 * delta_w ) )
# Adjust filter order up or down as needed
M += N_bump
N_taps = M + 1
# Obtain the Kaiser window
beta = signal . kaiser_beta ( d_stop )
w_k = signal . kaiser ( N_taps , beta )
n = np . arange ( N_taps )
b_k = wc / np . pi * np . sinc ( wc / np . pi * ( n - M / 2 ) ) * w_k
b_k /= np . sum ( b_k )
# Transform LPF to BPF
f0 = ( f_pass2 + f_pass1 ) / 2
w0 = 2 * np . pi * f0 / fs
n = np . arange ( len ( b_k ) )
b_k_bp = 2 * b_k * np . cos ( w0 * ( n - M / 2 ) )
print ( 'Kaiser Win filter taps = %d.' % N_taps )
return b_k_bp


signature:
def lowpass_order ( f_pass , f_stop , dpass_dB , dstop_dB , fsamp = 1 ) :
    """Optimal FIR ( equal ripple ) Lowpass Order Determination
 Text reference : Ifeachor , Digital Signal Processing a Practical Approach ,
 second edition , Prentice Hall , 2002.
 Journal paper reference : Herriman et al . , Practical Design Rules for Optimum
 Finite Imulse Response Digitl Filters , Bell Syst . Tech . J . , vol 52 , pp .
 769-799 , July - Aug . , 1973 . IEEE , 1973."""


signature:
def bandpass_order ( f_stop1 , f_pass1 , f_pass2 , f_stop2 , dpass_dB , dstop_dB , fsamp = 1 ) :
    """Optimal FIR ( equal ripple ) Bandpass Order Determination
 Text reference : Ifeachor , Digital Signal Processing a Practical Approach ,
 second edition , Prentice Hall , 2002.
 Journal paper reference : F . Mintzer & B . Liu , Practical Design Rules for Optimum
 FIR Bandpass Digital Filters , IEEE Transactions on Acoustics and Speech , pp .
 204-206 , April , 1979."""

implementation:
dpass = 1 - 10 ** ( - dpass_dB / 20 )
dstop = 10 ** ( - dstop_dB / 20 )
Df1 = ( f_pass1 - f_stop1 ) / fsamp
Df2 = ( f_stop2 - f_pass2 ) / fsamp
b1 = 0.01201
b2 = 0.09664
b3 = - 0.51325
b4 = 0.00203
b5 = - 0.5705
b6 = - 0.44314
Df = min ( Df1 , Df2 )
Cinf = np . log10 ( dstop ) * ( b1 * np . log10 ( dpass ) ** 2 + b2 * np . log10 ( dpass ) + b3 ) + ( b4 * np . log10 ( dpass ) ** 2 + b5 * np . log10 ( dpass ) + b6 )
g = - 14.6 * np . log10 ( dpass / dstop ) - 16.9
N = Cinf / Df + g * Df + 1
ff = 2 * np . array ( [ 0 , f_stop1 , f_pass1 , f_pass2 , f_stop2 , fsamp / 2 ] ) / fsamp
aa = np . array ( [ 0 , 0 , 1 , 1 , 0 , 0 ] )
wts = np . array ( [ dpass / dstop , 1 , dpass / dstop ] )
return int ( N ) , ff , aa , wts


signature:
def fir_remez_lpf ( f_pass , f_stop , d_pass , d_stop , fs = 1.0 , N_bump = 5 ) :
    """Design an FIR lowpass filter using remez with order
 determination . The filter order is determined based on
 f _ pass Hz , fstop Hz , and the desired passband ripple
 d _ pass dB and stopband attenuation d _ stop dB all
 relative to a sampling rate of fs Hz .
 Mark Wickert October 2016 , updated October 2018"""

implementation:
n , ff , aa , wts = lowpass_order ( f_pass , f_stop , d_pass , d_stop , fsamp = fs )
# Bump up the order by N _ bump to bring down the final d _ pass & d _ stop
N_taps = n
N_taps += N_bump
b = signal . remez ( N_taps , ff , aa [ 0 : : 2 ] , wts , Hz = 2 )
print ( 'Remez filter taps = %d.' % N_taps )
return b


signature:
def fir_remez_hpf ( f_stop , f_pass , d_pass , d_stop , fs = 1.0 , N_bump = 5 ) :
    """Design an FIR highpass filter using remez with order
 determination . The filter order is determined based on
 f _ pass Hz , fstop Hz , and the desired passband ripple
 d _ pass dB and stopband attenuation d _ stop dB all
 relative to a sampling rate of fs Hz .
 Mark Wickert October 2016 , updated October 2018"""

implementation:
# Transform HPF critical frequencies to lowpass equivalent
f_pass_eq = fs / 2. - f_pass
f_stop_eq = fs / 2. - f_stop
# Design LPF equivalent
n , ff , aa , wts = lowpass_order ( f_pass_eq , f_stop_eq , d_pass , d_stop , fsamp = fs )
# Bump up the order by N _ bump to bring down the final d _ pass & d _ stop
N_taps = n
N_taps += N_bump
b = signal . remez ( N_taps , ff , aa [ 0 : : 2 ] , wts , Hz = 2 )
# Transform LPF equivalent to HPF
n = np . arange ( len ( b ) )
b *= ( - 1 ) ** n
print ( 'Remez filter taps = %d.' % N_taps )
return b


signature:
def fir_remez_bpf ( f_stop1 , f_pass1 , f_pass2 , f_stop2 , d_pass , d_stop , fs = 1.0 , N_bump = 5 ) :
    """Design an FIR bandpass filter using remez with order
 determination . The filter order is determined based on
 f _ stop1 Hz , f _ pass1 Hz , f _ pass2 Hz , f _ stop2 Hz , and the
 desired passband ripple d _ pass dB and stopband attenuation
 d _ stop dB all relative to a sampling rate of fs Hz .
 Mark Wickert October 2016 , updated October 2018"""

implementation:
n , ff , aa , wts = bandpass_order ( f_stop1 , f_pass1 , f_pass2 , f_stop2 , d_pass , d_stop , fsamp = fs )
# Bump up the order by N _ bump to bring down the final d _ pass & d _ stop
N_taps = n
N_taps += N_bump
b = signal . remez ( N_taps , ff , aa [ 0 : : 2 ] , wts , Hz = 2 )
print ( 'Remez filter taps = %d.' % N_taps )
return b


signature:
def fir_remez_bsf ( f_pass1 , f_stop1 , f_stop2 , f_pass2 , d_pass , d_stop , fs = 1.0 , N_bump = 5 ) :
    """Design an FIR bandstop filter using remez with order
 determination . The filter order is determined based on
 f _ pass1 Hz , f _ stop1 Hz , f _ stop2 Hz , f _ pass2 Hz , and the
 desired passband ripple d _ pass dB and stopband attenuation
 d _ stop dB all relative to a sampling rate of fs Hz .
 Mark Wickert October 2016 , updated October 2018"""

implementation:
n , ff , aa , wts = bandstop_order ( f_pass1 , f_stop1 , f_stop2 , f_pass2 , d_pass , d_stop , fsamp = fs )
# Bump up the order by N _ bump to bring down the final d _ pass & d _ stop
# Initially make sure the number of taps is even so N _ bump needs to be odd
if np . mod ( n , 2 ) != 0 :
    n += 1
N_taps = n
N_taps += N_bump
b = signal . remez ( N_taps , ff , aa [ 0 : : 2 ] , wts , Hz = 2 , maxiter = 25 , grid_density = 16 )
print ( 'N_bump must be odd to maintain odd filter length' )
print ( 'Remez filter taps = %d.' % N_taps )
return b


signature:
def CIC ( M , K ) :
    """A functional form implementation of a cascade of integrator comb ( CIC ) filters .
 Parameters
 M : Effective number of taps per section ( typically the decimation factor ) .
 K : The number of CIC sections cascaded ( larger K gives the filter a wider image rejection bandwidth ) .
 Returns
 b : FIR filter coefficients for a simple direct form implementation using the filter ( ) function .
 Notes
 Commonly used in multirate signal processing digital down - converters and digital up - converters . A true CIC filter
 requires no multiplies , only add and subtract operations . The functional form created here is a simple FIR requiring
 real coefficient multiplies via filter ( ) .
 Mark Wickert July 2013"""

implementation:
if K == 1 :
    b = np . ones ( M )
else :
    h = np . ones ( M )
    b = h
    for i in range ( 1 , K ) :
        b = signal . convolve ( b , h )
# cascade by convolving impulse responses
# Make filter have unity gain at DC
return b / np . sum ( b )


signature:
def ten_band_eq_filt ( x , GdB , Q = 3.5 ) :
    """Filter the input signal x with a ten - band equalizer having octave gain values in ndarray GdB .
 The signal x is filtered using octave - spaced peaking filters starting at 31.25 Hz and
 stopping at 16 kHz . The Q of each filter is 3.5 , but can be changed . The sampling rate
 is assumed to be 44.1 kHz .
 Parameters
 x : ndarray of the input signal samples
 GdB : ndarray containing ten octave band gain values [ G0dB , . . . , G9dB ]
 Q : Quality factor vector for each of the NB peaking filters
 Returns
 y : ndarray of output signal samples
 Examples
 > > > # Test with white noise
 > > > w = randn ( 100000)
 > > > y = ten _ band _ eq _ filt ( x , GdB )
 > > > psd ( y , 2 * * 10,44.1)"""

implementation:
fs = 44100.0
# Hz
NB = len ( GdB )
if not NB == 10 :
    raise ValueError ( "GdB length not equal to ten" )
Fc = 31.25 * 2 ** np . arange ( NB )
B = np . zeros ( ( NB , 3 ) )
A = np . zeros ( ( NB , 3 ) )
# Create matrix of cascade coefficients
for k in range ( NB ) :
    [ b , a ] = peaking ( GdB [ k ] , Fc [ k ] , Q )
    B [ k , : ] = b
    A [ k , : ] = a
# Pass signal x through the cascade of ten filters
y = np . zeros ( len ( x ) )
for k in range ( NB ) :
    if k == 0 :
        y = signal . lfilter ( B [ k , : ] , A [ k , : ] , x )
    else :
        y = signal . lfilter ( B [ k , : ] , A [ k , : ] , y )
return y


signature:
def ten_band_eq_resp ( GdB , Q = 3.5 ) :
    """Create a frequency response magnitude plot in dB of a ten band equalizer
 using a semilogplot ( semilogx ( ) ) type plot
 Parameters
 GdB : Gain vector for 10 peaking filters [ G0 , . . . , G9]
 Q : Quality factor for each peaking filter ( default 3.5)
 Returns
 Nothing : two plots are created
 Examples
 > > > import matplotlib . pyplot as plt
 > > > from sk _ dsp _ comm import sigsys as ss
 > > > ss . ten _ band _ eq _ resp ( [ 0,10.0,0,0 , - 1,0,5,0 , - 4,0 ] )
 > > > plt . show ( )"""

implementation:
fs = 44100.0
# Hz
NB = len ( GdB )
if not NB == 10 :
    raise ValueError ( "GdB length not equal to ten" )
Fc = 31.25 * 2 ** np . arange ( NB )
B = np . zeros ( ( NB , 3 ) ) ;
A = np . zeros ( ( NB , 3 ) ) ;
# Create matrix of cascade coefficients
for k in range ( NB ) :
    b , a = peaking ( GdB [ k ] , Fc [ k ] , Q , fs )
    B [ k , : ] = b
    A [ k , : ] = a
# Create the cascade frequency response
F = np . logspace ( 1 , np . log10 ( 20e3 ) , 1000 )
H = np . ones ( len ( F ) ) * np . complex ( 1.0 , 0.0 )
for k in range ( NB ) :
    w , Htemp = signal . freqz ( B [ k , : ] , A [ k , : ] , 2 * np . pi * F / fs )
    H *= Htemp
plt . figure ( figsize = ( 6 , 4 ) )
plt . subplot ( 211 )
plt . semilogx ( F , 20 * np . log10 ( abs ( H ) ) )
plt . axis ( [ 10 , fs / 2 , - 12 , 12 ] )
plt . grid ( )
plt . title ( 'Ten-Band Equalizer Frequency Response' )
plt . xlabel ( 'Frequency (Hz)' )
plt . ylabel ( 'Gain (dB)' )
plt . subplot ( 212 )
plt . stem ( np . arange ( NB ) , GdB , 'b' , 'bs' )
# plt . bar ( np . arange ( NB ) - . 1 , GdB , 0.2)
plt . axis ( [ 0 , NB - 1 , - 12 , 12 ] )
plt . xlabel ( 'Equalizer Band Number' )
plt . ylabel ( 'Gain Set (dB)' )
plt . grid ( )


signature:
def peaking ( GdB , fc , Q = 3.5 , fs = 44100. ) :
    """A second - order peaking filter having GdB gain at fc and approximately
 and 0 dB otherwise .
 The filter coefficients returns correspond to a biquadratic system function
 containing five parameters .
 Parameters
 GdB : Lowpass gain in dB
 fc : Center frequency in Hz
 Q : Filter Q which is inversely proportional to bandwidth
 fs : Sampling frquency in Hz
 Returns
 b : ndarray containing the numerator filter coefficients
 a : ndarray containing the denominator filter coefficients
 Examples
 > > > import matplotlib . pyplot as plt
 > > > import numpy as np
 > > > from sk _ dsp _ comm . sigsys import peaking
 > > > from scipy import signal
 > > > b , a = peaking ( 2.0,500)
 > > > f = np . logspace ( 1,5,400)
 > > > w , H = signal . freqz ( b , a , 2 * np . pi * f / 44100)
 > > > plt . semilogx ( f , 20 * np . log10 ( abs ( H ) ) )
 > > > plt . ylabel ( " Power Spectral Density ( dB ) " )
 > > > plt . xlabel ( " Frequency ( Hz ) " )
 > > > plt . show ( )
 > > > b , a = peaking ( - 5.0,500,4)
 > > > w , H = signal . freqz ( b , a , 2 * np . pi * f / 44100)
 > > > plt . semilogx ( f , 20 * np . log10 ( abs ( H ) ) )
 > > > plt . ylabel ( " Power Spectral Density ( dB ) " )
 > > > plt . xlabel ( " Frequency ( Hz ) " )"""


signature:
def ex6_2 ( n ) :
    """Generate a triangle pulse as described in Example 6-2
 of Chapter 6.
 You need to supply an index array n that covers at least [ - 2 , 5 ] .
 The function returns the hard - coded signal of the example .
 Parameters
 n : time index ndarray covering at least - 2 to + 5.
 Returns
 x : ndarray of signal samples in x
 Examples
 > > > import numpy as np
 > > > import matplotlib . pyplot as plt
 > > > from sk _ dsp _ comm import sigsys as ss
 > > > n = np . arange ( - 5,8)
 > > > x = ss . ex6_2 ( n )
 > > > plt . stem ( n , x ) # creates a stem plot of x vs n"""

implementation:
x = np . zeros ( len ( n ) )
for k , nn in enumerate ( n ) :
    if nn >= - 2 and nn <= 5 :
        x [ k ] = 8 - nn
return x


signature:
def position_CD ( Ka , out_type = 'fb_exact' ) :
    """CD sled position control case study of Chapter 18.
 The function returns the closed - loop and open - loop
 system function for a CD / DVD sled position control
 system . The loop amplifier gain is the only variable
 that may be changed . The returned system function can
 however be changed .
 Parameters
 Ka : loop amplifier gain , start with 50.
 out _ type : ' open _ loop ' for open loop system function
 out _ type : ' fb _ approx ' for closed - loop approximation
 out _ type : ' fb _ exact ' for closed - loop exact
 Returns
 b : numerator coefficient ndarray
 a : denominator coefficient ndarray
 Notes
 With the exception of the loop amplifier gain , all
 other parameters are hard - coded from Case Study example .
 Examples
 > > > b , a = position _ CD ( Ka , ' fb _ approx ' )
 > > > b , a = position _ CD ( Ka , ' fb _ exact ' )"""

implementation:
rs = 10 / ( 2 * np . pi )
# Load b and a ndarrays with the coefficients
if out_type . lower ( ) == 'open_loop' :
    b = np . array ( [ Ka * 4000 * rs ] )
    a = np . array ( [ 1 , 1275 , 31250 , 0 ] )
elif out_type . lower ( ) == 'fb_approx' :
    b = np . array ( [ 3.2 * Ka * rs ] )
    a = np . array ( [ 1 , 25 , 3.2 * Ka * rs ] )
elif out_type . lower ( ) == 'fb_exact' :
    b = np . array ( [ 4000 * Ka * rs ] )
    a = np . array ( [ 1 , 1250 + 25 , 25 * 1250 , 4000 * Ka * rs ] )
else :
    raise ValueError ( 'out_type must be: open_loop, fb_approx, or fc_exact' )
return b , a


signature:
def cruise_control ( wn , zeta , T , vcruise , vmax , tf_mode = 'H' ) :
    """Cruise control with PI controller and hill disturbance .
 This function returns various system function configurations
 for a the cruise control Case Study example found in
 the supplementary article . The plant model is obtained by the
 linearizing the equations of motion and the controller contains a
 proportional and integral gain term set via the closed - loop parameters
 natuarl frequency wn ( rad / s ) and damping zeta .
 Parameters
 wn : closed - loop natural frequency in rad / s , nominally 0.1
 zeta : closed - loop damping factor , nominally 1.0
 T : vehicle time constant , nominally 10 s
 vcruise : cruise velocity set point , nominally 75 mph
 vmax : maximum vehicle velocity , nominally 120 mph
 tf _ mode : ' H ' , ' HE ' , ' HVW ' , or ' HED ' controls the system function returned by the function
 ' H ' : closed - loop system function V ( s ) / R ( s )
 ' HE ' : closed - loop system function E ( s ) / R ( s )
 ' HVW ' : closed - loop system function V ( s ) / W ( s )
 ' HED ' : closed - loop system function E ( s ) / D ( s ) , where D is the hill disturbance input
 Returns
 b : numerator coefficient ndarray
 a : denominator coefficient ndarray
 Examples
 > > > # return the closed - loop system function output / input velocity
 > > > b , a = cruise _ control ( wn , zeta , T , vcruise , vmax , tf _ mode = ' H ' )
 > > > # return the closed - loop system function loop error / hill disturbance
 > > > b , a = cruise _ control ( wn , zeta , T , vcruise , vmax , tf _ mode = ' HED ' )"""

implementation:
tau = T / 2. * vmax / vcruise
g = 9.8
g *= 3 * 60 ** 2 / 5280.
# m / s to mph conversion
Kp = T * ( 2 * zeta * wn - 1 / tau ) / vmax
Ki = T * wn ** 2. / vmax
K = Kp * vmax / T
print ( 'wn = ' , np . sqrt ( K / ( Kp / Ki ) ) )
print ( 'zeta = ' , ( K + 1 / tau ) / ( 2 * wn ) )
a = np . array ( [ 1 , 2 * zeta * wn , wn ** 2 ] )
if tf_mode == 'H' :
    b = np . array ( [ K , wn ** 2 ] )
elif tf_mode == 'HE' :
    b = np . array ( [ 1 , 2 * zeta * wn - K , 0. ] )
elif tf_mode == 'HVW' :
    b = np . array ( [ 1 , wn ** 2 / K + 1 / tau , wn ** 2 / ( K * tau ) ] )
    b *= Kp
elif tf_mode == 'HED' :
    b = np . array ( [ g , 0 ] )
else :
    raise ValueError ( 'tf_mode must be: H, HE, HVU, or HED' )
return b , a


signature:
def splane ( b , a , auto_scale = True , size = [ - 1 , 1 , - 1 , 1 ] ) :
    """Create an s - plane pole - zero plot .
 As input the function uses the numerator and denominator
 s - domain system function coefficient ndarrays b and a respectively .
 Assumed to be stored in descending powers of s .
 Parameters
 b : numerator coefficient ndarray .
 a : denominator coefficient ndarray .
 auto _ scale : True
 size : [ xmin , xmax , ymin , ymax ] plot scaling when scale = False
 Returns
 ( M , N ) : tuple of zero and pole counts + plot window
 Notes
 This function tries to identify repeated poles and zeros and will
 place the multiplicity number above and to the right of the pole or zero .
 The difficulty is setting the tolerance for this detection . Currently it
 is set at 1e - 3 via the function signal . unique _ roots .
 Examples
 > > > # Here the plot is generated using auto _ scale
 > > > splane ( b , a )
 > > > # Here the plot is generated using manual scaling
 > > > splane ( b , a , False , [ - 10,1 , - 10,10 ] )"""

implementation:
M = len ( b ) - 1
N = len ( a ) - 1
plt . figure ( figsize = ( 5 , 5 ) )
# plt . axis ( ' equal ' )
N_roots = np . array ( [ 0.0 ] )
if M > 0 :
    N_roots = np . roots ( b )
D_roots = np . array ( [ 0.0 ] )
if N > 0 :
    D_roots = np . roots ( a )
if auto_scale :
    size [ 0 ] = min ( np . min ( np . real ( N_roots ) ) , np . min ( np . real ( D_roots ) ) ) - 0.5
    size [ 1 ] = max ( np . max ( np . real ( N_roots ) ) , np . max ( np . real ( D_roots ) ) ) + 0.5
    size [ 1 ] = max ( size [ 1 ] , 0.5 )
    size [ 2 ] = min ( np . min ( np . imag ( N_roots ) ) , np . min ( np . imag ( D_roots ) ) ) - 0.5
    size [ 3 ] = max ( np . max ( np . imag ( N_roots ) ) , np . max ( np . imag ( D_roots ) ) ) + 0.5
plt . plot ( [ size [ 0 ] , size [ 1 ] ] , [ 0 , 0 ] , 'k--' )
plt . plot ( [ 0 , 0 ] , [ size [ 2 ] , size [ 3 ] ] , 'r--' )
# Plot labels if multiplicity greater than 1
x_scale = size [ 1 ] - size [ 0 ]
y_scale = size [ 3 ] - size [ 2 ]
x_off = 0.03
y_off = 0.01
if M > 0 : # N _ roots = np . roots ( b )
    N_uniq , N_mult = signal . unique_roots ( N_roots , tol = 1e-3 , rtype = 'avg' )
    plt . plot ( np . real ( N_uniq ) , np . imag ( N_uniq ) , 'ko' , mfc = 'None' , ms = 8 )
    idx_N_mult = np . nonzero ( np . ravel ( N_mult > 1 ) ) [ 0 ]
    for k in range ( len ( idx_N_mult ) ) :
        x_loc = np . real ( N_uniq [ idx_N_mult [ k ] ] ) + x_off * x_scale
        y_loc = np . imag ( N_uniq [ idx_N_mult [ k ] ] ) + y_off * y_scale
        plt . text ( x_loc , y_loc , str ( N_mult [ idx_N_mult [ k ] ] ) , ha = 'center' , va = 'bottom' , fontsize = 10 )
if N > 0 : # D _ roots = np . roots ( a )
    D_uniq , D_mult = signal . unique_roots ( D_roots , tol = 1e-3 , rtype = 'avg' )
    plt . plot ( np . real ( D_uniq ) , np . imag ( D_uniq ) , 'kx' , ms = 8 )
    idx_D_mult = np . nonzero ( np . ravel ( D_mult > 1 ) ) [ 0 ]
    for k in range ( len ( idx_D_mult ) ) :
        x_loc = np . real ( D_uniq [ idx_D_mult [ k ] ] ) + x_off * x_scale
        y_loc = np . imag ( D_uniq [ idx_D_mult [ k ] ] ) + y_off * y_scale
        plt . text ( x_loc , y_loc , str ( D_mult [ idx_D_mult [ k ] ] ) , ha = 'center' , va = 'bottom' , fontsize = 10 )
plt . xlabel ( 'Real Part' )
plt . ylabel ( 'Imaginary Part' )
plt . title ( 'Pole-Zero Plot' )
# plt . grid ( )
plt . axis ( np . array ( size ) )
return M , N


signature:
def OS_filter ( x , h , N , mode = 0 ) :
    """Overlap and save transform domain FIR filtering .
 This function implements the classical overlap and save method of
 transform domain filtering using a length P FIR filter .
 Parameters
 x : input signal to be filtered as an ndarray
 h : FIR filter coefficients as an ndarray of length P
 N : FFT size > P , typically a power of two
 mode : 0 or 1 , when 1 returns a diagnostic matrix
 Returns
 y : the filtered output as an ndarray
 y _ mat : an ndarray whose rows are the individual overlap outputs .
 Notes
 y _ mat is used for diagnostics and to gain understanding of the algorithm .
 Examples
 > > > n = arange ( 0,100)
 > > > x = cos ( 2 * pi * 0.05 * n )
 > > > b = ones ( 10)
 > > > y = OS _ filter ( x , h , N )
 > > > # set mode = 1
 > > > y , y _ mat = OS _ filter ( x , h , N , 1)"""

implementation:
P = len ( h )
# zero pad start of x so first frame can recover first true samples of x
x = np . hstack ( ( np . zeros ( P - 1 ) , x ) )
L = N - P + 1
Nx = len ( x )
Nframe = int ( np . ceil ( Nx / float ( L ) ) )
# zero pad end of x to full number of frames needed
x = np . hstack ( ( x , np . zeros ( Nframe * L - Nx ) ) )
y = np . zeros ( int ( Nframe * N ) )
# create an instrumentation matrix to observe the overlap and save behavior
y_mat = np . zeros ( ( Nframe , int ( Nframe * N ) ) )
H = fft . fft ( h , N )
# begin the filtering operation
for k in range ( Nframe ) :
    xk = x [ k * L : k * L + N ]
    Xk = fft . fft ( xk , N )
    Yk = H * Xk
    yk = np . real ( fft . ifft ( Yk ) )
    # imag part should be zero
    y [ k * L + P - 1 : k * L + N ] = yk [ P - 1 : ]
    y_mat [ k , k * L : k * L + N ] = yk
if mode == 1 :
    return y [ P - 1 : Nx ] , y_mat [ : , P - 1 : Nx ]
else :
    return y [ P - 1 : Nx ]


signature:
def lp_samp ( fb , fs , fmax , N , shape = 'tri' , fsize = ( 6 , 4 ) ) :
    """Lowpass sampling theorem plotting function .
 Display the spectrum of a sampled signal after setting the bandwidth ,
 sampling frequency , maximum display frequency , and spectral shape .
 Parameters
 fb : spectrum lowpass bandwidth in Hz
 fs : sampling frequency in Hz
 fmax : plot over [ - fmax , fmax ]
 shape : ' tri ' or ' line '
 N : number of translates , N positive and N negative
 fsize : the size of the figure window , default ( 6,4)
 Returns
 Nothing : A plot window opens containing the spectrum plot
 Examples
 > > > import matplotlib . pyplot as plt
 > > > from sk _ dsp _ comm . sigsys import lp _ samp
 No aliasing as bandwidth 10 Hz < 25/2 ; fs > fb .
 > > > lp _ samp ( 10,25,50,10)
 > > > plt . show ( )
 Now aliasing as bandwidth 15 Hz > 25/2 ; fs < fb .
 > > > lp _ samp ( 15,25,50,10)"""

implementation:
plt . figure ( figsize = fsize )
# define the plot interval
f = np . arange ( - fmax , fmax + fmax / 200. , fmax / 200. )
A = 1.0
line_ampl = A / 2. * np . array ( [ 0 , 1 ] )
# plot the lowpass spectrum in black
shapes = [ 'tri' , 'line' ]
if shape . lower ( ) not in shapes :
    raise ValueError ( 'shape must be tri or line' )
if shape . lower ( ) == 'tri' :
    plt . plot ( f , lp_tri ( f , fb ) )
    # overlay positive and negative frequency translates
    for n in range ( N ) :
        plt . plot ( f , lp_tri ( f - ( n + 1 ) * fs , fb ) , '--r' )
        plt . plot ( f , lp_tri ( f + ( n + 1 ) * fs , fb ) , '--g' )
elif shape . lower ( ) == 'line' :
    plt . plot ( [ fb , fb ] , line_ampl , 'b' , linewidth = 2 )
    plt . plot ( [ - fb , - fb ] , line_ampl , 'b' , linewidth = 2 )
    # overlay positive and negative frequency translates
    for n in range ( N ) :
        plt . plot ( [ fb + ( n + 1 ) * fs , fb + ( n + 1 ) * fs ] , line_ampl , '--r' , linewidth = 2 )
        plt . plot ( [ - fb + ( n + 1 ) * fs , - fb + ( n + 1 ) * fs ] , line_ampl , '--r' , linewidth = 2 )
        plt . plot ( [ fb - ( n + 1 ) * fs , fb - ( n + 1 ) * fs ] , line_ampl , '--g' , linewidth = 2 )
        plt . plot ( [ - fb - ( n + 1 ) * fs , - fb - ( n + 1 ) * fs ] , line_ampl , '--g' , linewidth = 2 )
plt . ylabel ( 'Spectrum Magnitude' )
plt . xlabel ( 'Frequency in Hz' )
plt . axis ( [ - fmax , fmax , 0 , 1 ] )
plt . grid ( )


signature:
def lp_tri ( f , fb ) :
    """Triangle spectral shape function used by : func : ` lp _ samp ` .
 Parameters
 f : ndarray containing frequency samples
 fb : the bandwidth as a float constant
 Returns
 x : ndarray of spectrum samples for a single triangle shape
 Notes
 This is a support function for the lowpass spectrum plotting function
 : func : ` lp _ samp ` .
 Examples
 > > > x = lp _ tri ( f , fb )"""

implementation:
x = np . zeros ( len ( f ) )
for k in range ( len ( f ) ) :
    if abs ( f [ k ] ) <= fb :
        x [ k ] = 1 - abs ( f [ k ] ) / float ( fb )
return x


signature:
def sinusoidAWGN ( x , SNRdB ) :
    """Add white Gaussian noise to a single real sinusoid .
 Input a single sinusoid to this function and it returns a noisy
 sinusoid at a specific SNR value in dB . Sinusoid power is calculated
 using np . var .
 Parameters
 x : Input signal as ndarray consisting of a single sinusoid
 SNRdB : SNR in dB for output sinusoid
 Returns
 y : Noisy sinusoid return vector
 Examples
 > > > # set the SNR to 10 dB
 > > > n = arange ( 0,10000)
 > > > x = cos ( 2 * pi * 0.04 * n )
 > > > y = sinusoidAWGN ( x , 10.0)"""

implementation:
# Estimate signal power
x_pwr = np . var ( x )
# Create noise vector
noise = np . sqrt ( x_pwr / 10 ** ( SNRdB / 10. ) ) * np . random . randn ( len ( x ) ) ;
return x + noise


signature:
def discrim ( x ) :
    """function disdata = discrim ( x )
 where x is an angle modulated signal in complex baseband form .
 Mark Wickert"""

implementation:
X = np . real ( x )
# X is the real part of the received signal
Y = np . imag ( x )
# Y is the imaginary part of the received signal
b = np . array ( [ 1 , - 1 ] )
# filter coefficients for discrete derivative
a = np . array ( [ 1 , 0 ] )
# filter coefficients for discrete derivative
derY = signal . lfilter ( b , a , Y )
# derivative of Y ,
derX = signal . lfilter ( b , a , X )
disdata = ( X * derY - Y * derX ) / ( X ** 2 + Y ** 2 )
return disdata


signature:
def mono_FM ( x , fs = 2.4e6 , file_name = 'test.wav' ) :
    """Decimate complex baseband input by 10
 Design 1st decimation lowpass filter ( f _ c = 200 KHz )"""

implementation:
b = signal . firwin ( 64 , 2 * 200e3 / float ( fs ) )
# Filter and decimate ( should be polyphase )
y = signal . lfilter ( b , 1 , x )
z = ss . downsample ( y , 10 )
# Apply complex baseband discriminator
z_bb = discrim ( z )
# Design 2nd decimation lowpass filter ( fc = 12 KHz )
bb = signal . firwin ( 64 , 2 * 12e3 / ( float ( fs ) / 10 ) )
# Filter and decimate
zz_bb = signal . lfilter ( bb , 1 , z_bb )
# Decimate by 5
z_out = ss . downsample ( zz_bb , 5 )
# Save to wave file
ss . to_wav ( file_name , 48000 , z_out / 2 )
print ( 'Done!' )
return z_bb , z_out


signature:
def stereo_FM ( x , fs = 2.4e6 , file_name = 'test.wav' ) :
    """Stereo demod from complex baseband at sampling rate fs .
 Assume fs is 2400 ksps
 Mark Wickert July 2017"""

implementation:
N1 = 10
b = signal . firwin ( 64 , 2 * 200e3 / float ( fs ) )
# Filter and decimate ( should be polyphase )
y = signal . lfilter ( b , 1 , x )
z = ss . downsample ( y , N1 )
# Apply complex baseband discriminator
z_bb = discrim ( z )
# Work with the ( 3 ) stereo multiplex signals :
# Begin by designing a lowpass filter for L + R and DSP demoded ( L - R )
# ( fc = 12 KHz )
b12 = signal . firwin ( 128 , 2 * 12e3 / ( float ( fs ) / N1 ) )
# The L + R term is at baseband , we just lowpass filter to remove
# other terms above 12 kHz .
y_lpr = signal . lfilter ( b12 , 1 , z_bb )
b19 = signal . firwin ( 128 , 2 * 1e3 * np . array ( [ 19 - 5 , 19 + 5 ] ) / ( float ( fs ) / N1 ) , pass_zero = False ) ;
z_bb19 = signal . lfilter ( b19 , 1 , z_bb )
# Lock PLL to 19 kHz pilot
# A type 2 loop with bandwidth Bn = 10 Hz and damping zeta = 0.707
# The VCO quiescent frequency is set to 19000 Hz .
theta , phi_error = pilot_PLL ( z_bb19 , 19000 , fs / N1 , 2 , 10 , 0.707 )
# Coherently demodulate the L - R subcarrier at 38 kHz .
# theta is the PLL output phase at 19 kHz , so to double multiply
# by 2 and wrap with cos ( ) or sin ( ) .
# First bandpass filter
b38 = signal . firwin ( 128 , 2 * 1e3 * np . array ( [ 38 - 5 , 38 + 5 ] ) / ( float ( fs ) / N1 ) , pass_zero = False ) ;
x_lmr = signal . lfilter ( b38 , 1 , z_bb )
# Coherently demodulate using the PLL output phase
x_lmr = 2 * np . sqrt ( 2 ) * np . cos ( 2 * theta ) * x_lmr
# Lowpass at 12 kHz to recover the desired DSB demod term
y_lmr = signal . lfilter ( b12 , 1 , x_lmr )
# Matrix the y _ lmr and y _ lpr for form right and left channels :
y_left = y_lpr + y_lmr
y_right = y_lpr - y_lmr
# Decimate by N2 ( nominally 5)
N2 = 5
fs2 = float ( fs ) / ( N1 * N2 )
# ( nominally 48 ksps )
y_left_DN2 = ss . downsample ( y_left , N2 )
y_right_DN2 = ss . downsample ( y_right , N2 )
# Deemphasize with 75 us time constant to ' undo ' the preemphasis
# applied at the transmitter in broadcast FM .
# A 1 - pole digital lowpass works well here .
a_de = np . exp ( - 2.1 * 1e3 * 2 * np . pi / fs2 )
z_left = signal . lfilter ( [ 1 - a_de ] , [ 1 , - a_de ] , y_left_DN2 )
z_right = signal . lfilter ( [ 1 - a_de ] , [ 1 , - a_de ] , y_right_DN2 )
# Place left and righ channels as side - by - side columns in a 2D array
z_out = np . hstack ( ( np . array ( [ z_left ] ) . T , ( np . array ( [ z_right ] ) . T ) ) )
ss . to_wav ( file_name , 48000 , z_out / 2 )
print ( 'Done!' )
# return z _ bb , z _ out
return z_bb , theta , y_lpr , y_lmr , z_out


signature:
def pilot_PLL ( xr , fq , fs , loop_type , Bn , zeta ) :
    """theta , phi _ error = pilot _ PLL ( xr , fq , fs , loop _ type , Bn , zeta )
 Mark Wickert , April 2014"""

implementation:
T = 1 / float ( fs )
# Set the VCO gain in Hz / V
Kv = 1.0
# Design a lowpass filter to remove the double freq term
Norder = 5
b_lp , a_lp = signal . butter ( Norder , 2 * ( fq / 2. ) / float ( fs ) )
fstate = np . zeros ( Norder )
# LPF state vector
Kv = 2 * np . pi * Kv
# convert Kv in Hz / v to rad / s / v
if loop_type == 1 : # First - order loop parameters
    fn = Bn
    Kt = 2 * np . pi * fn
    # loop natural frequency in rad / s
elif loop_type == 2 : # Second - order loop parameters
    fn = 1 / ( 2 * np . pi ) * 2 * Bn / ( zeta + 1 / ( 4 * zeta ) )
    # given Bn in Hz
    Kt = 4 * np . pi * zeta * fn
    # loop natural frequency in rad / s
    a = np . pi * fn / zeta
else :
    print ( 'Loop type must be 1 or 2' )
# Initialize integration approximation filters
filt_in_last = 0
filt_out_last = 0
vco_in_last = 0
vco_out = 0
vco_out_last = 0
# Initialize working and final output vectors
n = np . arange ( 0 , len ( xr ) )
theta = np . zeros ( len ( xr ) )
ev = np . zeros ( len ( xr ) )
phi_error = np . zeros ( len ( xr ) )
# Normalize total power in an attemp to make the 19kHz sinusoid
# component have amplitude ~ 1.
# xr = xr / ( 2/3 * std ( xr ) ) ;
# Begin the simulation loop
for kk in range ( len ( n ) ) : # Sinusoidal phase detector ( simple multiplier )
    phi_error [ kk ] = 2 * xr [ kk ] * np . sin ( vco_out )
    # LPF to remove double frequency term
    phi_error [ kk ] , fstate = signal . lfilter ( b_lp , a_lp , np . array ( [ phi_error [ kk ] ] ) , zi = fstate )
    pd_out = phi_error [ kk ]
    # pd _ out = 0
    # Loop gain
    gain_out = Kt / Kv * pd_out
    # apply VCO gain at VCO
    # Loop filter
    if loop_type == 2 :
        filt_in = a * gain_out
        filt_out = filt_out_last + T / 2. * ( filt_in + filt_in_last )
        filt_in_last = filt_in
        filt_out_last = filt_out
        filt_out = filt_out + gain_out
    else :
        filt_out = gain_out
    # VCO
    vco_in = filt_out + fq / ( Kv / ( 2 * np . pi ) )
    # bias to quiescent freq .
    vco_out = vco_out_last + T / 2. * ( vco_in + vco_in_last )
    vco_in_last = vco_in
    vco_out_last = vco_out
    vco_out = Kv * vco_out
    # apply Kv
    # Measured loop signals
    ev [ kk ] = filt_out
    theta [ kk ] = np . mod ( vco_out , 2 * np . pi ) ;
    # The vco phase mod 2pi
return theta , phi_error


signature:
def sccs_bit_sync ( y , Ns ) :
    """rx _ symb _ d , clk , track = sccs _ bit _ sync ( y , Ns )
 Symbol synchronization algorithm using SCCS
 y = baseband NRZ data waveform
 Ns = nominal number of samples per symbol
 Reworked from ECE 5675 Project
 Translated from m - code version
 Mark Wickert April 2014"""

implementation:
# decimated symbol sequence for SEP
rx_symb_d = np . zeros ( int ( np . fix ( len ( y ) / Ns ) ) )
track = np . zeros ( int ( np . fix ( len ( y ) / Ns ) ) )
bit_count = - 1
y_abs = np . zeros ( len ( y ) )
clk = np . zeros ( len ( y ) )
k = Ns + 1
# initial 1 - of - Ns symbol synch clock phase
# Sample - by - sample processing required
for i in range ( len ( y ) ) : # y _ abs ( i ) = abs ( round ( real ( y ( i ) ) ) )
    if i >= Ns : # do not process first Ns samples
    # Collect timing decision unit ( TDU ) samples
        y_abs [ i ] = np . abs ( np . sum ( y [ i - Ns + 1 : i + 1 ] ) )
        # Update sampling instant and take a sample
        # For causality reason the early sample is ' i ' ,
        # the on - time or prompt sample is ' i - 1 ' , and
        # the late sample is ' i - 2 ' .
        if ( k == 0 ) : # Load the samples into the 3x1 TDU register w _ hat .
        # w _ hat [ 1 ] = late , w _ hat [ 2 ] = on - time ; w _ hat [ 3 ] = early .
            w_hat = y_abs [ i - 2 : i + 1 ]
            bit_count += 1
            if w_hat [ 1 ] != 0 :
                if w_hat [ 0 ] < w_hat [ 2 ] :
                    k = Ns - 1
                    clk [ i - 2 ] = 1
                    rx_symb_d [ bit_count ] = y [ i - 2 - int ( np . round ( Ns / 2 ) ) - 1 ]
                elif w_hat [ 0 ] > w_hat [ 2 ] :
                    k = Ns + 1
                    clk [ i ] = 1
                    rx_symb_d [ bit_count ] = y [ i - int ( np . round ( Ns / 2 ) ) - 1 ]
                else :
                    k = Ns
                    clk [ i - 1 ] = 1
                    rx_symb_d [ bit_count ] = y [ i - 1 - int ( np . round ( Ns / 2 ) ) - 1 ]
            else :
                k = Ns
                clk [ i - 1 ] = 1
                rx_symb_d [ bit_count ] = y [ i - 1 - int ( np . round ( Ns / 2 ) ) ]
            track [ bit_count ] = np . mod ( i , Ns )
    k -= 1
# Trim the final output to bit _ count
rx_symb_d = rx_symb_d [ : bit_count ]
return rx_symb_d , clk , track


signature:
def fsk_BEP ( rx_data , m , flip ) :
    """fsk _ BEP ( rx _ data , m , flip )
 Estimate the BEP of the data bits recovered
 by the RTL - SDR Based FSK Receiver .
 The reference m - sequence generated in Python
 was found to produce sequences running in the opposite
 direction relative to the m - sequences generated by the
 mbed . To allow error detection the reference m - sequence
 is flipped .
 Mark Wickert April 2014"""

implementation:
Nbits = len ( rx_data )
c = dc . m_seq ( m )
if flip == 1 : # Flip the sequence to compenstate for mbed code difference
# First make it a 1xN array
    c . shape = ( 1 , len ( c ) )
    c = np . fliplr ( c ) . flatten ( )
L = int ( np . ceil ( Nbits / float ( len ( c ) ) ) )
tx_data = np . dot ( c . reshape ( len ( c ) , 1 ) , np . ones ( ( 1 , L ) ) )
tx_data = tx_data . T . reshape ( ( 1 , len ( c ) * L ) ) . flatten ( )
tx_data = tx_data [ : Nbits ]
# Convert to + 1 / - 1 bits
tx_data = 2 * tx_data - 1
Bit_count , Bit_errors = dc . BPSK_BEP ( rx_data , tx_data )
print ( 'len rx_data = %d, len tx_data = %d' % ( len ( rx_data ) , len ( tx_data ) ) )
Pe = Bit_errors / float ( Bit_count )
print ( '/////////////////////////////////////' )
print ( 'Bit Errors: %d' % Bit_errors )
print ( 'Bits Total: %d' % Bit_count )
print ( '       BEP: %2.2e' % Pe )
print ( '/////////////////////////////////////' )


signature:
def complex2wav ( filename , rate , x ) :
    """Save a complex signal vector to a wav file for compact binary
 storage of 16 - bit signal samples . The wav left and right channels
 are used to save real ( I ) and imaginary ( Q ) values . The rate is
 just a convent way of documenting the original signal sample rate .
 complex2wav ( filename , rate , x )
 Mark Wickert April 2014"""

implementation:
x_wav = np . hstack ( ( np . array ( [ x . real ] ) . T , np . array ( [ x . imag ] ) . T ) )
ss . to_wav ( filename , rate , x_wav )
print ( 'Saved as binary wav file with (I,Q)<=>(L,R)' )


signature:
def wav2complex ( filename ) :
    """Return a complex signal vector from a wav file that was used to store
 the real ( I ) and imaginary ( Q ) values of a complex signal ndarray .
 The rate is included as means of recalling the original signal sample
 rate .
 fs , x = wav2complex ( filename )
 Mark Wickert April 2014"""

implementation:
fs , x_LR_cols = ss . from_wav ( filename )
x = x_LR_cols [ : , 0 ] + 1j * x_LR_cols [ : , 1 ]
return fs , x


signature:
def FIR_header ( fname_out , h ) :
    """Write FIR Filter Header Files
 Mark Wickert February 2015"""


signature:
def FIR_fix_header ( fname_out , h ) :
    """Write FIR Fixed - Point Filter Header Files
 Mark Wickert February 2015"""


signature:
def IIR_sos_header ( fname_out , SOS_mat ) :
    """Write IIR SOS Header Files
 File format is compatible with CMSIS - DSP IIR
 Directform II Filter Functions
 Mark Wickert March 2015 - October 2016"""


signature:
def CA_code_header ( fname_out , Nca ) :
    """Write 1023 bit CA ( Gold ) Code Header Files
 Mark Wickert February 2015"""

implementation:
dir_path = os . path . dirname ( os . path . realpath ( __file__ ) )
ca = loadtxt ( dir_path + '/ca1thru37.txt' , dtype = int16 , usecols = ( Nca - 1 , ) , unpack = True )
M = 1023
# code period
N = 23
# code bits per line
Sca = 'ca' + str ( Nca )
f = open ( fname_out , 'wt' )
f . write ( '//define a CA code\n\n' )
f . write ( '#include <stdint.h>\n\n' )
f . write ( '#ifndef N_CA\n' )
f . write ( '#define N_CA %d\n' % M )
f . write ( '#endif\n' )
f . write ( '/*******************************************************************/\n' ) ;
f . write ( '/*                    1023 Bit CA Gold Code %2d                     */\n' % Nca ) ;
f . write ( 'int8_t ca%d[N_CA] = {' % Nca )
kk = 0 ;
for k in range ( M ) : # k _ mod = k % M
    if ( kk < N - 1 ) and ( k < M - 1 ) :
        f . write ( '%d,' % ca [ k ] )
        kk += 1
    elif ( kk == N - 1 ) & ( k < M - 1 ) :
        f . write ( '%d,\n' % ca [ k ] )
        if k < M :
            if Nca < 10 :
                f . write ( '                    ' )
            else :
                f . write ( '                     ' )
            kk = 0
    else :
        f . write ( '%d' % ca [ k ] )
f . write ( '};\n' )
f . write ( '/*******************************************************************/\n' )
f . close ( )


signature:
def farrow_resample ( x , fs_old , fs_new ) :
    """Parameters
 x : Input list representing a signal vector needing resampling .
 fs _ old : Starting / old sampling frequency .
 fs _ new : New sampling frequency .
 Returns
 y : List representing the signal vector resampled at the new frequency .
 Notes
 A cubic interpolator using a Farrow structure is used resample the
 input data at a new sampling rate that may be an irrational multiple of
 the input sampling rate .
 Time alignment can be found for a integer value M , found with the following :
 . . math : : f _ { s , out } = f _ { s , in } ( M - 1 ) / M
 The filter coefficients used here and a more comprehensive listing can be
 found in H . Meyr , M . Moeneclaey , & S . Fechtel , " Digital Communication
 Receivers , " Wiley , 1998 , Chapter 9 , pp . 521-523.
 Another good paper on variable interpolators is : L . Erup , F . Gardner , &
 R . Harris , " Interpolation in Digital Modems - - Part II : Implementation
 and Performance , " IEEE Comm . Trans . , June 1993 , pp . 998-1008.
 A founding paper on the subject of interpolators is : C . W . Farrow , " A
 Continuously variable Digital Delay Element , " Proceedings of the IEEE
 Intern . Symp . on Circuits Syst . , pp . 2641-2645 , June 1988.
 Mark Wickert April 2003 , recoded to Python November 2013
 Examples
 The following example uses a QPSK signal with rc pulse shaping , and time alignment at M = 15.
 > > > import matplotlib . pyplot as plt
 > > > from numpy import arange
 > > > from sk _ dsp _ comm import digitalcom as dc
 > > > Ns = 8
 > > > Rs = 1.
 > > > fsin = Ns * Rs
 > > > Tsin = 1 / fsin
 > > > N = 200
 > > > ts = 1
 > > > x , b , data = dc . MPSK _ bb ( N + 12 , Ns , 4 , ' rc ' )
 > > > x = x [ 12 * Ns : ]
 > > > xxI = x . real
 > > > M = 15
 > > > fsout = fsin * ( M - 1 ) / M
 > > > Tsout = 1 . / fsout
 > > > xI = dc . farrow _ resample ( xxI , fsin , fsin )
 > > > tx = arange ( 0 , len ( xI ) ) / fsin
 > > > yI = dc . farrow _ resample ( xxI , fsin , fsout )
 > > > ty = arange ( 0 , len ( yI ) ) / fsout
 > > > plt . plot ( tx - Tsin , xI )
 > > > plt . plot ( tx [ ts : : Ns ] - Tsin , xI [ ts : : Ns ] , ' r . ' )
 > > > plt . plot ( ty [ ts : : Ns ] - Tsout , yI [ ts : : Ns ] , ' g . ' )
 > > > plt . title ( r ' Impact of Asynchronous Sampling ' )
 > > > plt . ylabel ( r ' Real Signal Amplitude ' )
 > > > plt . xlabel ( r ' Symbol Rate Normalized Time ' )
 > > > plt . xlim ( [ 0 , 20 ] )
 > > > plt . grid ( )
 > > > plt . show ( )"""

implementation:
# Cubic interpolator over 4 samples .
# The base point receives a two sample delay .
v3 = signal . lfilter ( [ 1 / 6. , - 1 / 2. , 1 / 2. , - 1 / 6. ] , [ 1 ] , x )
v2 = signal . lfilter ( [ 0 , 1 / 2. , - 1 , 1 / 2. ] , [ 1 ] , x )
v1 = signal . lfilter ( [ - 1 / 6. , 1 , - 1 / 2. , - 1 / 3. ] , [ 1 ] , x )
v0 = signal . lfilter ( [ 0 , 0 , 1 ] , [ 1 ] , x )
Ts_old = 1 / float ( fs_old )
Ts_new = 1 / float ( fs_new )
T_end = Ts_old * ( len ( x ) - 3 )
t_new = np . arange ( 0 , T_end + Ts_old , Ts_new )
if x . dtype == np . dtype ( 'complex128' ) or x . dtype == np . dtype ( 'complex64' ) :
    y = np . zeros ( len ( t_new ) ) + 1j * np . zeros ( len ( t_new ) )
else :
    y = np . zeros ( len ( t_new ) )
for n in range ( len ( t_new ) ) :
    n_old = int ( np . floor ( n * Ts_new / Ts_old ) )
    mu = ( n * Ts_new - n_old * Ts_old ) / Ts_old
    # Combine outputs
    y [ n ] = ( ( v3 [ n_old + 1 ] * mu + v2 [ n_old + 1 ] ) * mu + v1 [ n_old + 1 ] ) * mu + v0 [ n_old + 1 ]
return y


signature:
def eye_plot ( x , L , S = 0 ) :
    """Eye pattern plot of a baseband digital communications waveform .
 The signal must be real , but can be multivalued in terms of the underlying
 modulation scheme . Used for BPSK eye plots in the Case Study article .
 Parameters
 x : ndarray of the real input data vector / array
 L : display length in samples ( usually two symbols )
 S : start index
 Returns
 None : A plot window opens containing the eye plot
 Notes
 Increase S to eliminate filter transients .
 Examples
 1000 bits at 10 samples per bit with ' rc ' shaping .
 > > > import matplotlib . pyplot as plt
 > > > from sk _ dsp _ comm import digitalcom as dc
 > > > x , b , data = dc . NRZ _ bits ( 1000,10 , ' rc ' )
 > > > dc . eye _ plot ( x , 20,60)
 > > > plt . show ( )"""

implementation:
plt . figure ( figsize = ( 6 , 4 ) )
idx = np . arange ( 0 , L + 1 )
plt . plot ( idx , x [ S : S + L + 1 ] , 'b' )
k_max = int ( ( len ( x ) - S ) / L ) - 1
for k in range ( 1 , k_max ) :
    plt . plot ( idx , x [ S + k * L : S + L + 1 + k * L ] , 'b' )
plt . grid ( )
plt . xlabel ( 'Time Index - n' )
plt . ylabel ( 'Amplitude' )
plt . title ( 'Eye Plot' )
return 0


signature:
def scatter ( x , Ns , start ) :
    """Sample a baseband digital communications waveform at the symbol spacing .
 Parameters
 x : ndarray of the input digital comm signal
 Ns : number of samples per symbol ( bit )
 start : the array index to start the sampling
 Returns
 xI : ndarray of the real part of x following sampling
 xQ : ndarray of the imaginary part of x following sampling
 Notes
 Normally the signal is complex , so the scatter plot contains
 clusters at point in the complex plane . For a binary signal
 such as BPSK , the point centers are nominally + / - 1 on the real
 axis . Start is used to eliminate transients from the FIR
 pulse shaping filters from appearing in the scatter plot .
 Examples
 > > > import matplotlib . pyplot as plt
 > > > from sk _ dsp _ comm import digitalcom as dc
 > > > x , b , data = dc . NRZ _ bits ( 1000,10 , ' rc ' )
 Add some noise so points are now scattered about + / - 1.
 > > > y = dc . cpx _ AWGN ( x , 20,10)
 > > > yI , yQ = dc . scatter ( y , 10,60)
 > > > plt . plot ( yI , yQ , ' . ' )
 > > > plt . grid ( )
 > > > plt . xlabel ( ' In - Phase ' )
 > > > plt . ylabel ( ' Quadrature ' )
 > > > plt . axis ( ' equal ' )
 > > > plt . show ( )"""

implementation:
xI = np . real ( x [ start : : Ns ] )
xQ = np . imag ( x [ start : : Ns ] )
return xI , xQ


signature:
def strips ( x , Nx , fig_size = ( 6 , 4 ) ) :
    """Plots the contents of real ndarray x as a vertical stacking of
 strips , each of length Nx . The default figure size is ( 6,4 ) inches .
 The yaxis tick labels are the starting index of each strip . The red
 dashed lines correspond to zero amplitude in each strip .
 strips ( x , Nx , my _ figsize = ( 6,4 ) )
 Mark Wickert April 2014"""

implementation:
plt . figure ( figsize = fig_size )
# ax = fig . add _ subplot ( 111)
N = len ( x )
Mx = int ( np . ceil ( N / float ( Nx ) ) )
x_max = np . max ( np . abs ( x ) )
for kk in range ( Mx ) :
    plt . plot ( np . array ( [ 0 , Nx ] ) , - kk * Nx * np . array ( [ 1 , 1 ] ) , 'r-.' )
    plt . plot ( x [ kk * Nx : ( kk + 1 ) * Nx ] / x_max * 0.4 * Nx - kk * Nx , 'b' )
plt . axis ( [ 0 , Nx , - Nx * ( Mx - 0.5 ) , Nx * 0.5 ] )
plt . yticks ( np . arange ( 0 , - Nx * Mx , - Nx ) , np . arange ( 0 , Nx * Mx , Nx ) )
plt . xlabel ( 'Index' )
plt . ylabel ( 'Strip Amplitude and Starting Index' )
return 0


signature:
def bit_errors ( tx_data , rx_data , Ncorr = 1024 , Ntransient = 0 ) :
    """Count bit errors between a transmitted and received BPSK signal .
 Time delay between streams is detected as well as ambiquity resolution
 due to carrier phase lock offsets of : math : ` k * \\ pi ` , k = 0,1.
 The ndarray tx _ data is Tx 0/1 bits as real numbers I .
 The ndarray rx _ data is Rx 0/1 bits as real numbers I .
 Note : Ncorr needs to be even"""

implementation:
# Remove Ntransient symbols and level shift to { - 1 , + 1}
tx_data = 2 * tx_data [ Ntransient : ] - 1
rx_data = 2 * rx_data [ Ntransient : ] - 1
# Correlate the first Ncorr symbols at four possible phase rotations
R0 = np . fft . ifft ( np . fft . fft ( rx_data , Ncorr ) * np . conj ( np . fft . fft ( tx_data , Ncorr ) ) )
R1 = np . fft . ifft ( np . fft . fft ( - 1 * rx_data , Ncorr ) * np . conj ( np . fft . fft ( tx_data , Ncorr ) ) )
# Place the zero lag value in the center of the array
R0 = np . fft . fftshift ( R0 )
R1 = np . fft . fftshift ( R1 )
R0max = np . max ( R0 . real )
R1max = np . max ( R1 . real )
R = np . array ( [ R0max , R1max ] )
Rmax = np . max ( R )
kphase_max = np . where ( R == Rmax ) [ 0 ]
kmax = kphase_max [ 0 ]
# Correlation lag value is zero at the center of the array
if kmax == 0 :
    lagmax = np . where ( R0 . real == Rmax ) [ 0 ] - Ncorr / 2
elif kmax == 1 :
    lagmax = np . where ( R1 . real == Rmax ) [ 0 ] - Ncorr / 2
taumax = lagmax [ 0 ]
print ( 'kmax =  %d, taumax = %d' % ( kmax , taumax ) )
# Count bit and symbol errors over the entire input ndarrays
# Begin by making tx and rx length equal and apply phase rotation to rx
if taumax < 0 :
    tx_data = tx_data [ int ( - taumax ) : ]
    tx_data = tx_data [ : min ( len ( tx_data ) , len ( rx_data ) ) ]
    rx_data = ( - 1 ) ** kmax * rx_data [ : len ( tx_data ) ]
else :
    rx_data = ( - 1 ) ** kmax * rx_data [ int ( taumax ) : ]
    rx_data = rx_data [ : min ( len ( tx_data ) , len ( rx_data ) ) ]
    tx_data = tx_data [ : len ( rx_data ) ]
# Convert to 0 ' s and 1 ' s
Bit_count = len ( tx_data )
tx_I = np . int16 ( ( tx_data . real + 1 ) / 2 )
rx_I = np . int16 ( ( rx_data . real + 1 ) / 2 )
Bit_errors = tx_I ^ rx_I
return Bit_count , np . sum ( Bit_errors )


signature:
def QAM_bb ( N_symb , Ns , mod_type = '16qam' , pulse = 'rect' , alpha = 0.35 ) :
    """QAM _ BB _ TX : A complex baseband transmitter
 x , b , tx _ data = QAM _ bb ( K , Ns , M )
 / / / / / Inputs / / / / /
 N _ symb = the number of symbols to process
 Ns = number of samples per symbol
 mod _ type = modulation type : qpsk , 16qam , 64qam , or 256qam
 alpha = squareroot raised codine pulse shape bandwidth factor .
 For DOCSIS alpha = 0.12 to 0.18 . In general alpha can
 range over 0 < alpha < 1.
 SRC = pulse shape : 0 - > rect , 1 - > SRC
 / / / / / Outputs / / / / /
 x = complex baseband digital modulation
 b = transmitter shaping filter , rectangle or SRC
 tx _ data = xI + 1j * xQ = inphase symbol sequence +
 1j * quadrature symbol sequence
 Mark Wickert November 2014"""

implementation:
# Filter the impulse train waveform with a square root raised
# cosine pulse shape designed as follows :
# Design the filter to be of duration 12 symbols and
# fix the excess bandwidth factor at alpha = 0.35
# If SRC = 0 use a simple rectangle pulse shape
if pulse . lower ( ) == 'src' :
    b = sqrt_rc_imp ( Ns , alpha , 6 )
elif pulse . lower ( ) == 'rc' :
    b = rc_imp ( Ns , alpha , 6 )
elif pulse . lower ( ) == 'rect' :
    b = np . ones ( int ( Ns ) )
    # alt . rect . pulse shape
else :
    raise ValueError ( 'pulse shape must be src, rc, or rect' )
if mod_type . lower ( ) == 'qpsk' :
    M = 2
    # bits per symbol
elif mod_type . lower ( ) == '16qam' :
    M = 4
elif mod_type . lower ( ) == '64qam' :
    M = 8
elif mod_type . lower ( ) == '256qam' :
    M = 16
else :
    raise ValueError ( 'Unknown mod_type' )
# Create random symbols for the I & Q channels
xI = np . random . randint ( 0 , M , N_symb )
xI = 2 * xI - ( M - 1 )
xQ = np . random . randint ( 0 , M , N_symb )
xQ = 2 * xQ - ( M - 1 )
# Employ differential encoding to counter phase ambiquities
# Create a zero padded ( interpolated by Ns ) symbol sequence .
# This prepares the symbol sequence for arbitrary pulse shaping .
symbI = np . hstack ( ( xI . reshape ( N_symb , 1 ) , np . zeros ( ( N_symb , int ( Ns ) - 1 ) ) ) )
symbI = symbI . flatten ( )
symbQ = np . hstack ( ( xQ . reshape ( N_symb , 1 ) , np . zeros ( ( N_symb , int ( Ns ) - 1 ) ) ) )
symbQ = symbQ . flatten ( )
symb = symbI + 1j * symbQ
if M > 2 :
    symb /= ( M - 1 )
# The impulse train waveform contains one pulse per Ns ( or Ts ) samples
# imp _ train = [ ones ( K , 1 ) zeros ( K , Ns - 1 ) ] ' ;
# imp _ train = reshape ( imp _ train , Ns * K , 1 ) ;
# Filter the impulse train signal
x = signal . lfilter ( b , 1 , symb )
x = x . flatten ( )
# out is a 1D vector
# Scale shaping filter to have unity DC gain
b = b / sum ( b )
return x , b , xI + 1j * xQ


signature:
def QAM_SEP ( tx_data , rx_data , mod_type , Ncorr = 1024 , Ntransient = 0 , SEP_disp = True ) :
    """Nsymb , Nerr , SEP _ hat =
 QAM _ symb _ errors ( tx _ data , rx _ data , mod _ type , Ncorr = 1024 , Ntransient = 0)
 Count symbol errors between a transmitted and received QAM signal .
 The received symbols are assumed to be soft values on a unit square .
 Time delay between streams is detected .
 The ndarray tx _ data is Tx complex symbols .
 The ndarray rx _ data is Rx complex symbols .
 Note : Ncorr needs to be even"""

implementation:
# Remove Ntransient symbols and makes lengths equal
tx_data = tx_data [ Ntransient : ]
rx_data = rx_data [ Ntransient : ]
Nmin = min ( [ len ( tx_data ) , len ( rx_data ) ] )
tx_data = tx_data [ : Nmin ]
rx_data = rx_data [ : Nmin ]
# Perform level translation and quantize the soft symbol values
if mod_type . lower ( ) == 'qpsk' :
    M = 2
    # bits per symbol
elif mod_type . lower ( ) == '16qam' :
    M = 4
elif mod_type . lower ( ) == '64qam' :
    M = 8
elif mod_type . lower ( ) == '256qam' :
    M = 16
else :
    raise ValueError ( 'Unknown mod_type' )
rx_data = np . rint ( ( M - 1 ) * ( rx_data + ( 1 + 1j ) ) / 2. )
# Fix - up edge points real part
s1r = np . nonzero ( np . ravel ( rx_data . real > M - 1 ) ) [ 0 ]
s2r = np . nonzero ( np . ravel ( rx_data . real < 0 ) ) [ 0 ]
rx_data . real [ s1r ] = ( M - 1 ) * np . ones ( len ( s1r ) )
rx_data . real [ s2r ] = np . zeros ( len ( s2r ) )
# Fix - up edge points imag part
s1i = np . nonzero ( np . ravel ( rx_data . imag > M - 1 ) ) [ 0 ]
s2i = np . nonzero ( np . ravel ( rx_data . imag < 0 ) ) [ 0 ]
rx_data . imag [ s1i ] = ( M - 1 ) * np . ones ( len ( s1i ) )
rx_data . imag [ s2i ] = np . zeros ( len ( s2i ) )
rx_data = 2 * rx_data - ( M - 1 ) * ( 1 + 1j )
# Correlate the first Ncorr symbols at four possible phase rotations
R0 , lags = xcorr ( rx_data , tx_data , Ncorr )
R1 , lags = xcorr ( rx_data * ( 1j ) ** 1 , tx_data , Ncorr )
R2 , lags = xcorr ( rx_data * ( 1j ) ** 2 , tx_data , Ncorr )
R3 , lags = xcorr ( rx_data * ( 1j ) ** 3 , tx_data , Ncorr )
# Place the zero lag value in the center of the array
R0max = np . max ( R0 . real )
R1max = np . max ( R1 . real )
R2max = np . max ( R2 . real )
R3max = np . max ( R3 . real )
R = np . array ( [ R0max , R1max , R2max , R3max ] )
Rmax = np . max ( R )
kphase_max = np . where ( R == Rmax ) [ 0 ]
kmax = kphase_max [ 0 ]
# Find correlation lag value is zero at the center of the array
if kmax == 0 :
    lagmax = lags [ np . where ( R0 . real == Rmax ) [ 0 ] ]
elif kmax == 1 :
    lagmax = lags [ np . where ( R1 . real == Rmax ) [ 0 ] ]
elif kmax == 2 :
    lagmax = lags [ np . where ( R2 . real == Rmax ) [ 0 ] ]
elif kmax == 3 :
    lagmax = lags [ np . where ( R3 . real == Rmax ) [ 0 ] ]
taumax = lagmax [ 0 ]
if SEP_disp :
    print ( 'Phase ambiquity = (1j)**%d, taumax = %d' % ( kmax , taumax ) )
# Count symbol errors over the entire input ndarrays
# Begin by making tx and rx length equal and apply
# phase rotation to rx _ data
if taumax < 0 :
    tx_data = tx_data [ - taumax : ]
    tx_data = tx_data [ : min ( len ( tx_data ) , len ( rx_data ) ) ]
    rx_data = ( 1j ) ** kmax * rx_data [ : len ( tx_data ) ]
else :
    rx_data = ( 1j ) ** kmax * rx_data [ taumax : ]
    rx_data = rx_data [ : min ( len ( tx_data ) , len ( rx_data ) ) ]
    tx_data = tx_data [ : len ( rx_data ) ]
# Convert QAM symbol difference to symbol errors
errors = np . int16 ( abs ( rx_data - tx_data ) )
# Detect symbols errors
# Could decode bit errors from symbol index difference
idx = np . nonzero ( np . ravel ( errors != 0 ) ) [ 0 ]
if SEP_disp :
    print ( 'Symbols = %d, Errors %d, SEP = %1.2e' % ( len ( errors ) , len ( idx ) , len ( idx ) / float ( len ( errors ) ) ) )
return len ( errors ) , len ( idx ) , len ( idx ) / float ( len ( errors ) )


signature:
def GMSK_bb ( N_bits , Ns , MSK = 0 , BT = 0.35 ) :
    """MSK / GMSK Complex Baseband Modulation
 x , data = gmsk ( N _ bits , Ns , BT = 0.35 , MSK = 0)
 Parameters
 N _ bits : number of symbols processed
 Ns : the number of samples per bit
 MSK : 0 for no shaping which is standard MSK , MSK < > 0 - - > GMSK is generated .
 BT : premodulation Bb * T product which sets the bandwidth of the Gaussian lowpass filter
 Mark Wickert Python version November 2014"""

implementation:
x , b , data = NRZ_bits ( N_bits , Ns )
# pulse length 2 * M * Ns
M = 4
n = np . arange ( - M * Ns , M * Ns + 1 )
p = np . exp ( - 2 * np . pi ** 2 * BT ** 2 / np . log ( 2 ) * ( n / float ( Ns ) ) ** 2 ) ;
p = p / np . sum ( p ) ;
# Gaussian pulse shape if MSK not zero
if MSK != 0 :
    x = signal . lfilter ( p , 1 , x )
y = np . exp ( 1j * np . pi / 2 * np . cumsum ( x ) / Ns )
return y , data


signature:
def MPSK_bb ( N_symb , Ns , M , pulse = 'rect' , alpha = 0.25 , MM = 6 ) :
    """Generate a complex baseband MPSK signal with pulse shaping .
 Parameters
 N _ symb : number of MPSK symbols to produce
 Ns : the number of samples per bit ,
 M : MPSK modulation order , e . g . , 4 , 8 , 16 , . . .
 pulse _ type : ' rect ' , ' rc ' , ' src ' ( default ' rect ' )
 alpha : excess bandwidth factor ( default 0.25)
 MM : single sided pulse duration ( default = 6)
 Returns
 x : ndarray of the MPSK signal values
 b : ndarray of the pulse shape
 data : ndarray of the underlying data bits
 Notes
 Pulse shapes include ' rect ' ( rectangular ) , ' rc ' ( raised cosine ) ,
 ' src ' ( root raised cosine ) . The actual pulse length is 2 * M + 1 samples .
 This function is used by BPSK _ tx in the Case Study article .
 Examples
 > > > from sk _ dsp _ comm import digitalcom as dc
 > > > import scipy . signal as signal
 > > > import matplotlib . pyplot as plt
 > > > x , b , data = dc . MPSK _ bb ( 500,10,8 , ' src ' , 0.35)
 > > > # Matched filter received signal x
 > > > y = signal . lfilter ( b , 1 , x )
 > > > plt . plot ( y . real [ 12*10 : ] , y . imag [ 12*10 : ] )
 > > > plt . xlabel ( ' In - Phase ' )
 > > > plt . ylabel ( ' Quadrature ' )
 > > > plt . axis ( ' equal ' )
 > > > # Sample once per symbol
 > > > plt . plot ( y . real [ 12*10 : : 10 ] , y . imag [ 12*10 : : 10 ] , ' r . ' )
 > > > plt . show ( )"""

implementation:
data = np . random . randint ( 0 , M , N_symb )
xs = np . exp ( 1j * 2 * np . pi / M * data )
x = np . hstack ( ( xs . reshape ( N_symb , 1 ) , np . zeros ( ( N_symb , int ( Ns ) - 1 ) ) ) )
x = x . flatten ( )
if pulse . lower ( ) == 'rect' :
    b = np . ones ( int ( Ns ) )
elif pulse . lower ( ) == 'rc' :
    b = rc_imp ( Ns , alpha , MM )
elif pulse . lower ( ) == 'src' :
    b = sqrt_rc_imp ( Ns , alpha , MM )
else :
    raise ValueError ( 'pulse type must be rec, rc, or src' )
x = signal . lfilter ( b , 1 , x )
if M == 4 :
    x = x * np . exp ( 1j * np . pi / 4 ) ;
    # For QPSK points in quadrants
return x , b / float ( Ns ) , data


signature:
def QPSK_rx ( fc , N_symb , Rs , EsN0 = 100 , fs = 125 , lfsr_len = 10 , phase = 0 , pulse = 'src' ) :
    """This function generates"""

implementation:
Ns = int ( np . round ( fs / Rs ) )
print ( 'Ns = ' , Ns )
print ( 'Rs = ' , fs / float ( Ns ) )
print ( 'EsN0 = ' , EsN0 , 'dB' )
print ( 'phase = ' , phase , 'degrees' )
print ( 'pulse = ' , pulse )
x , b , data = QPSK_bb ( N_symb , Ns , lfsr_len , pulse )
# Add AWGN to x
x = cpx_AWGN ( x , EsN0 , Ns )
n = np . arange ( len ( x ) )
xc = x * np . exp ( 1j * 2 * np . pi * fc / float ( fs ) * n ) * np . exp ( 1j * phase )
return xc , b , data


signature:
def QPSK_BEP ( tx_data , rx_data , Ncorr = 1024 , Ntransient = 0 ) :
    """Count bit errors between a transmitted and received QPSK signal .
 Time delay between streams is detected as well as ambiquity resolution
 due to carrier phase lock offsets of : math : ` k * \\ frac { \\ pi } { 4 } ` , k = 0,1,2,3.
 The ndarray sdata is Tx + / - 1 symbols as complex numbers I + j * Q .
 The ndarray data is Rx + / - 1 symbols as complex numbers I + j * Q .
 Note : Ncorr needs to be even"""

implementation:
# Remove Ntransient symbols
tx_data = tx_data [ Ntransient : ]
rx_data = rx_data [ Ntransient : ]
# Correlate the first Ncorr symbols at four possible phase rotations
R0 = np . fft . ifft ( np . fft . fft ( rx_data , Ncorr ) * np . conj ( np . fft . fft ( tx_data , Ncorr ) ) )
R1 = np . fft . ifft ( np . fft . fft ( 1j * rx_data , Ncorr ) * np . conj ( np . fft . fft ( tx_data , Ncorr ) ) )
R2 = np . fft . ifft ( np . fft . fft ( - 1 * rx_data , Ncorr ) * np . conj ( np . fft . fft ( tx_data , Ncorr ) ) )
R3 = np . fft . ifft ( np . fft . fft ( - 1j * rx_data , Ncorr ) * np . conj ( np . fft . fft ( tx_data , Ncorr ) ) )
# Place the zero lag value in the center of the array
R0 = np . fft . fftshift ( R0 )
R1 = np . fft . fftshift ( R1 )
R2 = np . fft . fftshift ( R2 )
R3 = np . fft . fftshift ( R3 )
R0max = np . max ( R0 . real )
R1max = np . max ( R1 . real )
R2max = np . max ( R2 . real )
R3max = np . max ( R3 . real )
R = np . array ( [ R0max , R1max , R2max , R3max ] )
Rmax = np . max ( R )
kphase_max = np . where ( R == Rmax ) [ 0 ]
kmax = kphase_max [ 0 ]
# Correlation lag value is zero at the center of the array
if kmax == 0 :
    lagmax = np . where ( R0 . real == Rmax ) [ 0 ] - Ncorr / 2
elif kmax == 1 :
    lagmax = np . where ( R1 . real == Rmax ) [ 0 ] - Ncorr / 2
elif kmax == 2 :
    lagmax = np . where ( R2 . real == Rmax ) [ 0 ] - Ncorr / 2
elif kmax == 3 :
    lagmax = np . where ( R3 . real == Rmax ) [ 0 ] - Ncorr / 2
taumax = lagmax [ 0 ]
print ( 'kmax =  %d, taumax = %d' % ( kmax , taumax ) )
# Count bit and symbol errors over the entire input ndarrays
# Begin by making tx and rx length equal and apply phase rotation to rx
if taumax < 0 :
    tx_data = tx_data [ - taumax : ]
    tx_data = tx_data [ : min ( len ( tx_data ) , len ( rx_data ) ) ]
    rx_data = 1j ** kmax * rx_data [ : len ( tx_data ) ]
else :
    rx_data = 1j ** kmax * rx_data [ taumax : ]
    rx_data = rx_data [ : min ( len ( tx_data ) , len ( rx_data ) ) ]
    tx_data = tx_data [ : len ( rx_data ) ]
# Convert to 0 ' s and 1 ' s
S_count = len ( tx_data )
tx_I = np . int16 ( ( tx_data . real + 1 ) / 2 )
tx_Q = np . int16 ( ( tx_data . imag + 1 ) / 2 )
rx_I = np . int16 ( ( rx_data . real + 1 ) / 2 )
rx_Q = np . int16 ( ( rx_data . imag + 1 ) / 2 )
I_errors = tx_I ^ rx_I
Q_errors = tx_Q ^ rx_Q
# A symbol errors occurs when I or Q or both are in error
S_errors = I_errors | Q_errors
# return 0
return S_count , np . sum ( I_errors ) , np . sum ( Q_errors ) , np . sum ( S_errors )


signature:
def BPSK_tx ( N_bits , Ns , ach_fc = 2.0 , ach_lvl_dB = - 100 , pulse = 'rect' , alpha = 0.25 , M = 6 ) :
    """Generates biphase shift keyed ( BPSK ) transmitter with adjacent channel interference .
 Generates three BPSK signals with rectangular or square root raised cosine ( SRC )
 pulse shaping of duration N _ bits and Ns samples per bit . The desired signal is
 centered on f = 0 , which the adjacent channel signals to the left and right
 are also generated at dB level relative to the desired signal . Used in the
 digital communications Case Study supplement .
 Parameters
 N _ bits : the number of bits to simulate
 Ns : the number of samples per bit
 ach _ fc : the frequency offset of the adjacent channel signals ( default 2.0)
 ach _ lvl _ dB : the level of the adjacent channel signals in dB ( default - 100)
 pulse : the pulse shape ' rect ' or ' src '
 alpha : square root raised cosine pulse shape factor ( default = 0.25)
 M : square root raised cosine pulse truncation factor ( default = 6)
 Returns
 x : ndarray of the composite signal x0 + ach _ lvl * ( x1p + x1m )
 b : the transmit pulse shape
 data0 : the data bits used to form the desired signal ; used for error checking
 Notes
 Examples
 > > > x , b , data0 = BPSK _ tx ( 1000,10 , ' src ' )"""

implementation:
x0 , b , data0 = NRZ_bits ( N_bits , Ns , pulse , alpha , M )
x1p , b , data1p = NRZ_bits ( N_bits , Ns , pulse , alpha , M )
x1m , b , data1m = NRZ_bits ( N_bits , Ns , pulse , alpha , M )
n = np . arange ( len ( x0 ) )
x1p = x1p * np . exp ( 1j * 2 * np . pi * ach_fc / float ( Ns ) * n )
x1m = x1m * np . exp ( - 1j * 2 * np . pi * ach_fc / float ( Ns ) * n )
ach_lvl = 10 ** ( ach_lvl_dB / 20. )
return x0 + ach_lvl * ( x1p + x1m ) , b , data0


signature:
def rc_imp ( Ns , alpha , M = 6 ) :
    """A truncated raised cosine pulse used in digital communications .
 The pulse shaping factor : math : ` 0 < \\ alpha < 1 ` is required as well as the
 truncation factor M which sets the pulse duration to be : math : ` 2 * M * T _ { symbol } ` .
 Parameters
 Ns : number of samples per symbol
 alpha : excess bandwidth factor on ( 0 , 1 ) , e . g . , 0.35
 M : equals RC one - sided symbol truncation factor
 Returns
 b : ndarray containing the pulse shape
 See Also
 sqrt _ rc _ imp
 Notes
 The pulse shape b is typically used as the FIR filter coefficients
 when forming a pulse shaped digital communications waveform .
 Examples
 Ten samples per symbol and : math : ` \\ alpha = 0.35 ` .
 > > > import matplotlib . pyplot as plt
 > > > from sk _ dsp _ comm . digitalcom import rc _ imp
 > > > from numpy import arange
 > > > b = rc _ imp ( 10,0.35)
 > > > n = arange ( - 10*6,10*6 + 1)
 > > > plt . stem ( n , b )
 > > > plt . show ( )"""

implementation:
# Design the filter
n = np . arange ( - M * Ns , M * Ns + 1 )
b = np . zeros ( len ( n ) )
a = alpha
Ns *= 1.0
for i in range ( len ( n ) ) :
    if ( 1 - 4 * ( a * n [ i ] / Ns ) ** 2 ) == 0 :
        b [ i ] = np . pi / 4 * np . sinc ( 1 / ( 2. * a ) )
    else :
        b [ i ] = np . sinc ( n [ i ] / Ns ) * np . cos ( np . pi * a * n [ i ] / Ns ) / ( 1 - 4 * ( a * n [ i ] / Ns ) ** 2 )
return b


signature:
def sqrt_rc_imp ( Ns , alpha , M = 6 ) :
    """A truncated square root raised cosine pulse used in digital communications .
 The pulse shaping factor : math : ` 0 < \\ alpha < 1 ` is required as well as the
 truncation factor M which sets the pulse duration to be : math : ` 2 * M * T _ { symbol } ` .
 Parameters
 Ns : number of samples per symbol
 alpha : excess bandwidth factor on ( 0 , 1 ) , e . g . , 0.35
 M : equals RC one - sided symbol truncation factor
 Returns
 b : ndarray containing the pulse shape
 Notes
 The pulse shape b is typically used as the FIR filter coefficients
 when forming a pulse shaped digital communications waveform . When
 square root raised cosine ( SRC ) pulse is used to generate Tx signals and
 at the receiver used as a matched filter ( receiver FIR filter ) , the
 received signal is now raised cosine shaped , thus having zero
 intersymbol interference and the optimum removal of additive white
 noise if present at the receiver input .
 Examples
 Ten samples per symbol and : math : ` \\ alpha = 0.35 ` .
 > > > import matplotlib . pyplot as plt
 > > > from numpy import arange
 > > > from sk _ dsp _ comm . digitalcom import sqrt _ rc _ imp
 > > > b = sqrt _ rc _ imp ( 10,0.35)
 > > > n = arange ( - 10*6,10*6 + 1)
 > > > plt . stem ( n , b )
 > > > plt . show ( )"""