Best Approach to Rank Complex Magnitude Comparision Problem

Here are some of the latest results:

Algorithm         Correct    Time      Score    Penalties  Eggs
---------------   -------    ------    -------  ---------  ----
  Empty Economy    49.86     2.8104     472849    2378650    0
   Empty Deluxe     0.05     2.8311       1944  474168000  243
Starter Economy    89.75     2.9663     851367     486060    0
 Starter Deluxe    90.68     2.9764    1663118     441920  151

    Dan Beast 4    99.85     3.2622    1750076       7130  151
Cedron Unrolled   100.00     3.2721    1898616          0  243
  Cedron Deluxe   100.00     3.3255    1898616          0  243
 Cedron Revised   100.00     3.2128    1898616          0  243
   Olli Revised    99.50     3.1893    1728065      23880    0
  Olli Original    99.50     3.2464    1728065      23880    0

Cedron Multiply   100.00     3.2042    1898616          0  243
  Matt Multiply   100.00     3.3146    1898616          0  243

The timing for the contenders is too close and too noisy to show a clear favorite.  Benchmarking on the target platform would be much more useful now.

The code has been updated.  It is as it is.

import numpy as np
import timeit

# The passed arguments to a running horse.
#
#   ( I1, Q1 ) First  Complex Value (or Point)
#   ( I2, Q2 ) Second Complex Value (or Point)

# Its return values are
#
#   ( rc ) Comparison Result (Return Code)
#   (  l ) Locus of the Exit

# The return value can be one of these
#
#    -2  The first is for sure less than the second
#    -1  The first is likely less than the second
#     0  The two are equal for sure
#     1  The first is likely greater than the second
#     2  The first is for sure greater than the second
#
# Routines that only return {-1,1} can be called Economy
# Routines that only return {-2,0,2} can be called Deluxe
#

# How Scoring works
#
#   S Score
#   P Penalties
#   E Egg Count
#   W Wrong
#
#                    Correct      Marginal     Wrong
#       {-1,1}        S+=2         S+=1         W+=1,P+=10
#       {-2,0,2}      S+=4(E+=1)   S+=2,P+=10   W+=1,P+=1000
#
#

#====================================================================
#====================================================================
#    W A L K    O N S
#====================================================================
#====================================================================
def WalkOnOne( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
def WalkOnTwo( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
def WalkOnThree( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
#====================================================================
#    S T A R T E R    C O D E
#====================================================================
#====================================================================
def EmptyEconomy( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
def EmptyDeluxe( I1, Q1, I2, Q2 ):

        return 0, 0

#====================================================================
def StarterEconomy( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Return Results

        if x1 < x2:
           return -1, 0

        return 1, 0

#====================================================================
def StarterDeluxe( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              return 2, 0
           thePresumedResult = 1
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              return -2, 0
           thePresumedResult = -1
        else:
           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Return the Presumed Result

        return thePresumedResult, 2

#====================================================================
#====================================================================
#    C E D R O N ' S
#====================================================================
#====================================================================
def CedronRevised( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:

           if x1 + y1 >= x2 + y2:
              return 2, 0

           thePresumedResult = 2
           dx = x1 - x2

        elif x1 < x2:

           if x1 + y1 <= x2 + y2:
              return -2, 0

           thePresumedResult = -2
           dx = x2 - x1

        else:

           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- Bring Factors into 1/2 to 1 Ratio Range

        while dx <  sx:
              dx += dx

              if dy <  sy:
                 dy += dy
              else:
                 sy += sy

        while dy <  sy:
              dy += dy

              if dx <  sx:
                 dx += dx
              else:
                 sx += sx

#---- Use Double Arithmetic Mean as Proxy for Geometric Mean

        cx = sx + dx
        cy = sy + dy

        cx16 = cx << 4
        cy16 = cy << 4

        if cx16 - cx > cy16:
           return thePresumedResult, 2

        if cy16 - cy > cx16:
           return -thePresumedResult, 2

#---- X Multiplication

        px = 0

        while sx > 0:
           if sx & 1:
              px += dx

           dx += dx
           sx >>= 1

#---- Y Multiplication

        py = 0

        while sy > 0:
           if sy & 1:
              py += dy

           dy += dy
           sy >>= 1

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
def CedronUnrolled( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:

           if x1 + y1 >= x2 + y2:
              return 2, 0

           thePresumedResult = 2
           dx = x1 - x2

        elif x1 < x2:

           if x1 + y1 <= x2 + y2:
              return -2, 0

           thePresumedResult = -2
           dx = x2 - x1

        else:

           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Estimate First Multiplied Magnitude

        if y1 < (x1>>1):
           if y1 < (x1>>2):
              m1 = (x1<<8) - (x1<<1) \
                 + (y1<<5) + (y1<<1)
           else:
              m1 = (x1<<8) - (x1<<4) \
                 + (y1<<6) + (y1<<5) - (y1<<2) - (y1<<1)
        else:
           if y1 < (x1>>1) + (x1>>2):
              m1 = (x1<<8) - (x1<<5) - (x1<<2) - (x1<<1) \
                 + (y1<<7) + (y1<<3) - y1
           else:
              m1 = (x1<<7) + (x1<<6) + (x1<<1) \
                 + (y1<<7) + (y1<<5) + (y1<<3)

#---- Estimate Second Multiplied Magnitude

        if y2 < (x2>>1):
           if y2 <   (x2>>2):
              m2 = ( (x2<<7) - x2 \
                   + (y2<<4) + y2 ) << 1
           else:
              m2 = ( (x2<<7) - (x2<<3) \
                   + (y2<<5) + (y2<<4) - (y2<<1) - y2 ) << 1
        else:
           if y2 <   (x2>>1) + (x2>>2):
              m2 = ( (x2<<8) - (x2<<5) - (x2<<2) - (x2<<1) \
                   + (y2<<7) + (y2<<3) - y2 )
           else:
              m2 = ( (x2<<6) + (x2<<5) + x2 \
                   + (y2<<6) + (y2<<4) + (y2<<2) ) << 1

#---- Return Results (1000 is a temp hack value!)

        if m1 > m2 + (m2>>6):
           return 2, 2

        if m2 > m1 + (m1>>6):
           return -2, 2

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- X Multiplication

        px = 0

        while dx > 0:
           if dx & 1:
              px += sx

           sx += sx
           dx >>= 1

#---- Y Multiplication

        py = 0

        while dy > 0:
           if dy & 1:
              py += sy

           sy += sy
           dy >>= 1

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
def CedronDeluxe( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              return 2, 0
           dx = x1 - x2
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              return -2, 0
           dx = x2 - x1
        else:
           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Employ a DanBeast

        L1 = DanBeast_2_8_Level( x1, y1 )
        L2 = DanBeast_2_8_Level( x2, y2 )

#---- Early Out Return

        if L1 > L2 + (L2>>6):
           return 2, 1

        if L2 > L1 + (L1>>6):
           return -2, 1

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- Do the Multiplications

        px = UnsignedBitMultiply( sx, dx )
        py = UnsignedBitMultiply( sy, dy )

#---- Account for Swap

        if x1 > x2:
           thePresumedResult = 2
        else:
           thePresumedResult = -2

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
def DanBeastFour( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with Quick Exit

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              return 2, 0
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              return -2, 0
        else:
           if y1 > y2:
              return 2, 0
           elif y1 < y2:
              return -2, 0
           else:
              return 0, 0

#---- Estimate First Multiplied Magnitude

        if y1 < (x1>>1):
           if y1 < (x1>>2):
              m1 = (x1<<8) - (x1<<1) \
                 + (y1<<5) + (y1<<1)
           else:
              m1 = (x1<<8) - (x1<<4) \
                 + (y1<<6) + (y1<<5) - (y1<<2) - (y1<<1)
        else:
           if y1 < (x1>>1) + (x1>>2):
              m1 = (x1<<8) - (x1<<5) - (x1<<2) - (x1<<1) \
                 + (y1<<7) + (y1<<3) - y1
           else:
              m1 = (x1<<7) + (x1<<6) + (x1<<1) \
                 + (y1<<7) + (y1<<5) + (y1<<3)

#---- Estimate Second Multiplied Magnitude

        if y2 < (x2>>1):
           if y2 < (x2>>2):
              m2 = ( (x2<<7) - x2 \
                   + (y2<<4) + y2 ) << 1
           else:
              m2 = ( (x2<<7) - (x2<<3) \
                   + (y2<<5) + (y2<<4) - (y2<<1) - y2 ) << 1
        else:
           if y2 < (x2>>1) + (x2>>2):
              m2 = ( (x2<<8) - (x2<<5) - (x2<<2) - (x2<<1) \
                   + (y2<<7) + (y2<<3) - y2 )
           else:
              m2 = ( (x2<<6) + (x2<<5) + x2 \
                   + (y2<<6) + (y2<<4) + (y2<<2) ) << 1

#---- Return Results

        if m1 < m2:
           return -1, 2

        return 1, 2

#====================================================================
def CedronMultiply( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:

           if x1 + y1 >= x2 + y2:
              return 2, 0

           thePresumedResult = 2
           dx = x1 - x2

        elif x1 < x2:

           if x1 + y1 <= x2 + y2:
              return -2, 0

           thePresumedResult = -2
           dx = x2 - x1

        else:

           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- X Multiplication

        px = 0

        while dx > 0:
          if dx & 1:
             px += sx

          sx += sx
          dx >>= 1

#---- Y Multiplication

        py = 0

        while dy > 0:
          if dy & 1:
             py += sy

          sy += sy
          dy >>= 1

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
#====================================================================
#    O L L I    L I K E
#====================================================================
#====================================================================
def MyVersionOfOllis( I1, Q1, I2, Q2 ):

# Returns  ( c )
#
#    c   Comparison
#
#   -1     | (I1,Q1) |  <  | (I2,Q2) |
#    1     | (I1,Q1) |  >  | (I2,Q2) |
#
#    t   Exit Test
#
#    1     (Partial) Primary Determination
#    2     CORDIC Loop + 1
#    6     Terminating Guess

#---- Set Extent Parameter

        maxIterations = 4

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- (Partial) Primary Determination

        if I1 < I2 and  I1 + Q1 <= I2 + Q2:
           return -2, 1

        if I1 > I2 and  I1 + Q1 >= I2 + Q2:
           return 2, 1

#---- CORDIC Loop

        for n in range ( 1, maxIterations+1 ):
            newI1 = I1 + ( Q1 >> n )
            newQ1 = Q1 - ( I1 >> n )
            newI2 = I2 + ( Q2 >> n )
            newQ2 = Q2 - ( I2 >> n )

            I1 = newI1
            Q1 = abs( newQ1 )
            I2 = newI2
            Q2 = abs( newQ2 )

            s = n + n + 1

            if I1 <= I2 - ( I2 >> s ):
               return -1, 1 + n

            if I2 <= I1 - ( I1 >> s ):
               return 1, 1 + n

#---- Terminating Guess

        if I1 < I2:
           return -1, 7

        return 1, 7

#====================================================================
def MyRevisionOfOllis( I1, Q1, I2, Q2 ):

# Returns  ( rc, l )
#
#    c   Comparison
#
#   -1,-2   | (I1,Q1) |  <  | (I2,Q2) |
#    1, 2   | (I1,Q1) |  >  | (I2,Q2) |
#
#    t   Exit Test
#
#    1     (Partial) Primary Determination
#    2     CORDIC Loop + 1
#    6     Terminating Guess

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- (Partial) Primary Determination

        if I1 < I2 and  I1 + Q1 <= I2 + Q2:
           return -2, 1

        if I1 > I2 and  I1 + Q1 >= I2 + Q2:
           return 2, 1

#---- CORDIC Loop Head

        s = 3

        for n in range ( 1, 5 ):

#---- Apply the Rotation

            newI1 = I1 + ( Q1 >> n )
            newQ1 = Q1 - ( I1 >> n )

            newI2 = I2 + ( Q2 >> n )
            newQ2 = Q2 - ( I2 >> n )

#---- Attempt Comparison

            if newI1 <= newI2 - ( newI2 >> s ):
               return -1, 1 + n

            if newI2 <= newI1 - ( newI1 >> s ):
               return 1, 1 + n

            s += 2

#---- Advance the Values

            I1 = newI1
            I2 = newI2

            Q1 = abs( newQ1 )
            Q2 = abs( newQ2 )

#---- Terminating Guess

        if I1 < I2:
           return -1, 7

        return 1, 7

#====================================================================
#====================================================================
#    M A T T   L    L I K E
#====================================================================
#====================================================================
def MattMultiply( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- Ensure the first value is rightmost

        swap = 0;

        if I2 > I1:
           swap = 4
           I1, I2 = I2, I1
           Q1, Q2 = Q2, Q1

#---- Primary determination

        if I1 + Q1 > I2 + Q2:
           return 2 - swap, 2
        else:
           DI = I1 - I2
           if DI < 0:
              tmp1 = -UnsignedBitMultiply( I1 + I2, -DI )
           else:
              tmp1 =  UnsignedBitMultiply( I1 + I2,  DI  )

           DQ = Q2 - Q1
           if DQ < 0:
              tmp2 = -UnsignedBitMultiply( Q1 + Q2, -DQ )
           else:
              tmp2 =  UnsignedBitMultiply( Q1 + Q2,  DQ  )

           if tmp1 == tmp2:
              return  0       , 2
           elif tmp1 > tmp2:
              return  2 - swap, 2
           else:
              return -2 + swap, 2

#====================================================================
#====================================================================
#    U T I L I T I E S
#====================================================================
#====================================================================
def UnsignedBitMultiply( a, b ):  # Smaller value second is faster.

        p = 0

        while b > 0:
           if b & 1:
              p += a

           a  += a
           b >>= 1

        return p

#====================================================================
def DanBeast_2_8_Level( x, y ):

        if         y+y   <  x:               # 2 y < x
           if     (y<<2) <  x:               # 4 y < x
              L = (x<<8)   -x-x \
                + (y<<5)   +y+y                  # y/x = 0.00 to 0.25
           else:
              L = (x<<8) - (x<<4) \
                + (y<<6) + (y<<5) - (y<<2) -y-y  # y/x = 0.25 to 0.50
        else:
            if    (y<<2) <  x+x+x:           # 4 y < 3 x
              L = (x<<8) - (x<<5) - (x<<2) -x-x \
                + (y<<7) + (y<<3) -  y           # y/x = 0.50 to 0.75
            else:
              L = (x<<7) + (x<<6)   +x+x \
                + (y<<7) + (y<<5) + (y<<3)       # y/x = 0.75 to 1.00

        return L

#====================================================================
#====================================================================
#    T E S T I N G    H A R N E S S
#====================================================================
#====================================================================
def Test( ArgLimit, ArgThreshold, ArgLane, ArgTestName ):

#---- Set the Parameters

        t = ArgThreshold

#---- Initialize the Counters

        theCount           = 0
        theWrongCount      = 0

        theEggs            = 0
        theScore           = 0
        thePenalties       = 0

#---- Start Timing

        theStartTime = timeit.default_timer()

#---- Test on a Swept Area

        for i1 in range( -ArgLimit, ArgLimit, 10 ):
          ii1 = i1 * i1
          for q1 in range( -ArgLimit, ArgLimit, 7 ):
            d1 = np.sqrt( ii1 + q1 * q1 )
            for i2 in range( -ArgLimit, ArgLimit, 11 ):
              ii2 = i2 * i2
              for q2 in range( -ArgLimit, ArgLimit, 5 ):
                d2 = np.sqrt( ii2 + q2 * q2 )

                D = d1 - d2    #  = |(I1,Q1)| - |(I2,Q2)|

                theCount += 1

#---- The Fast Side Bench Mark Lanes

                if ArgLane == 0:
                   rc, l = EmptyEconomy( i1, q1, i2, q2 )

                if ArgLane == 1:
                   rc, l = EmptyDeluxe( i1, q1, i2, q2 )

                if ArgLane == 2:
                   rc, l = StarterEconomy( i1, q1, i2, q2 )

                if ArgLane == 3:
                   rc, l = StarterDeluxe( i1, q1, i2, q2 )

#---- The Slower Pace Horses

                if ArgLane == 8:
                   rc, l = TwoMultiply( i1, q1, i2, q2 )

                if ArgLane == 9:
                   rc, l = FourMultiply( i1, q1, i2, q2 )

#---- Walk Ons

                if ArgLane == 11:
                   rc, l = WalkOnOne( i1, q1, i2, q2 )

                if ArgLane == 12:
                   rc, l = WalkOnTwo( i1, q1, i2, q2 )

                if ArgLane == 13:
                   rc, l = WalkOnThree( i1, q1, i2, q2 )

#---- Cedron D.'s Lanes

                if ArgLane == 20:
                   rc, l = CedronRevised( i1, q1, i2, q2 )

                if ArgLane == 21:
                   rc, l = CedronDeluxe( i1, q1, i2, q2 )

                if ArgLane == 22:
                   rc, l = CedronUnrolled( i1, q1, i2, q2 )

                if ArgLane == 23:
                   rc, l = DanBeastFour( i1, q1, i2, q2 )

                if ArgLane == 24:
                   rc, l = CedronMultiply( i1, q1, i2, q2 )

#---- Olli N.'s Lanes

                if ArgLane == 30:
                   rc, l = MyVersionOfOllis( i1, q1, i2, q2 )

                if ArgLane == 31:
                   rc, l = MyRevisionOfOllis( i1, q1, i2, q2 )

#---- Dan B.'s Lanes

#                if ArgLane == 41:
#                   rc, l = Dan1( i1, q1, i2, q2 )

#---- Matt L.'s Lanes

                if ArgLane == 50:
                   rc, l = MattMultiply( i1, q1, i2, q2 )

#---- Assess Scores, Penalties, and Eggs

                if rc == -2:
                   if D < -t:
                      theScore      += 4
                   elif D < 0:
                      theScore      += 2
                      thePenalties  += 10
                   else:
                      theWrongCount += 1
                      thePenalties  += 1000

                elif rc == 2:
                   if D > t:
                      theScore      += 4
                   elif D > 0:
                      theScore      += 2
                      thePenalties  += 10
                   else:
                      theWrongCount += 1
                      thePenalties  += 1000

                elif rc == -1:
                   if D < 0:
                      theScore      += 2
                   elif D <= t:
                      theScore      += 1
                   else:
                      theWrongCount += 1
                      thePenalties  += 10

                elif rc == 1:
                   if D > 0:
                      theScore      += 2
                   elif D >= -t:
                      theScore      += 1
                   else:
                      theWrongCount += 1
                      thePenalties  += 10

                elif rc == 0:
                   if abs( D ) <= t:
                      theScore      += 8
                      if D == 0:
                         theEggs    += 1
                   else:
                      theWrongCount += 1
                      thePenalties  += 1000

                else:
                   print "Disqualification -- Invalid c value:", c, "Lane", ArgLane
                   return

#---- Finish Timing

        theDuration = timeit.default_timer() - theStartTime

#---- Calculate the Results

        theCorrectCount = theCount - theWrongCount

        theCorrectPct   = 100.0 * float( theCorrectCount ) \
                                / float( theCount )

#---- Return Results

        return "%15s  %7.2f %10.4f %10d %10d %4d" % \
               ( ArgTestName, theCorrectPct, theDuration,\
                 theScore,    thePenalties,  theEggs )

#====================================================================
def Main():

#---- Set Run Time Parameters

        L = 101   # The Limit
        T = 0     # Threshold

#---- Print Headers

        print "Algorithm         Correct    Time      Score    Penalties  Eggs"
        print "---------------   -------    ------    -------  ---------  ----"

#---- The Calibrators

        print Test( L, T, 0, "Empty Economy" )
        print Test( L, T, 1, "Empty Deluxe" )
        print Test( L, T, 2, "Starter Economy" )
        print Test( L, T, 3, "Starter Deluxe" )

#---- The Walk Ons

#        print
#        print Test( L, T, 11, "Walk On One" )

#---- The Contenders

        print
        print Test( L, T, 23, "Dan Beast 4" )
        print Test( L, T, 22, "Cedron Unrolled" )
        print Test( L, T, 21, "Cedron Deluxe" )
        print Test( L, T, 20, "Cedron Revised" )
        print Test( L, T, 31, "Olli Revised" )
        print Test( L, T, 30, "Olli Original" )

#---- The Pace Setters

        print
        print Test( L, T, 24, "Cedron Multiply" )
        print Test( L, T, 50, "Matt Multiply" )

#====================================================================
Main()

Earlier, I pledged a 50 point bounty to the best horse (fastest time 99%+ correct) that wasn't one of my entries.  I'm sticking with that, and right now Olli is leading.  (My optimized version is DQ'd)


Are you sure you shouldn't be using numpy.correlate instead of numpy.convolve?  To estimate delay, you want to cross-correlate your signals, not convolve them.  You'll possibly end up with a much larger delay by convolving.

Trying something simple:

x = [1, 0, 0, 0, 0 ];
y = [0, 0, 0, 0, 1 ];
conv = numpy.convolve(x,y);
conv
array([0, 0, 0, 0, 1, 0, 0, 0, 0])
corr = numpy.correlate(x,y,"full");
corr
array([1, 0, 0, 0, 0, 0, 0, 0, 0])


You will see lines in the spectrum of a phase-modulated signal if there is some periodicity in the modulating signal. If the modulating signal is not periodic, then the power spectrum of the modulated signal will be smooth.

As an example: consider a BPSK signal:

If the bits used to modulate the waveform are equally likely to be a one or a zero and are uncorrelated with one another, then there is no periodicity in the modulating signal. The power spectrum of a BPSK signal with independent equally-likely bits has a $sinc^2$ shape with no discrete lines.
If you then changed the data bits to add some repetitive framing data, for instance, causing a pattern of invariant bits in the modulating signal, you will begin to see lines appear due to the deterministic periodicity in the modulated signal.

In practice, one often makes an effort to ensure that bits are equally distributed to ones and zeros without repetitive patterns. Adding a source encoder (e.g. data compression) to your system can help to remove any repetitiveness in your data and whiten the transmitted bitstream.

One caveat to the above, courtesy of Dilip Sarwate's comment below, is duplicated here to make it more visible:

  However, modulating signals such as $\cos(\omega_1 t)+\cos(\omega_2 t)$ where $\omega_1 / \omega_2$ is an irrational number (i.e. no harmonic relationship exists) are not periodic but have lines in the spectrum at $\omega_c \pm m \omega_1$, $\omega_c \pm m \omega_2$, and $\omega_c \pm m \omega_1 \pm n \omega_2$ plus a lot of clutter (smooth spectrum) in between the lines. Here $\omega_c$ denotes the carrier frequency. There are also similar lines (and clutter) near the harmonics of the carrier frequency.


How to apply a Butterworth filter to data of varying sample rate?

I can figure out two options:

Interpolate data
Remember that filters are dynamical systems and solve the corresponding ordinary differential equation (ODE)

In the second method, you basically use the retrieved time step as the time step in the ODE solver. Assuming you have a second order system:

$$\ddot{y}+\frac{\omega_0}{Q}\dot{y}+y\omega_0^2=x\omega_0^2$$

For butterworth, $Q=\frac{1}{\sqrt{2}}$ Introduce two variables $u=y$ and $v=\dot{u}$, and you get two coupled first order ODE:

$$\left\{\begin{aligned}
     \dot{u}&=v\\
     \dot{v}&=x\omega_0^2-u\omega_0^2-\frac{\omega_0}{Q}v
    \end{aligned}\right.$$

Assuming you do not want to have a too high cutoff frequency

$$\left\{\begin{aligned}
   u&\leftarrow u+v\Delta t\\
   v&\leftarrow v+ \left(x\omega_0^2-u\omega_0^2-\frac{\omega_0}{Q}v\right)\Delta t\end{aligned}\right.$$

Where $\Delta t$ is the time difference between two steps. If you want $\omega_0 > \frac{0.25\text{ Hz}}{2\pi}$, you need to consider the Euler backward method, or if you need more accuracy, the trapezoid method.

Addendum

You may not want the overshoots and undershoots of the Butterworth filter. Then, set $Q$ a bit lower. Probably, you want a critically damped system, which do not have any such effects. In that case, set $Q=\frac{1}{2}$. An interesting side-effect of choosing the ODE approach is that you get a filtered derivative for free (it is the $v$ component), if you were to trig on changes in acceleration.


Lower peak to average power ratio by phase shifting

Each symbol (in the narrow sense - bit) in your frame modulates a sine wave of frequency depends on symbol's number and phase depends of symbol's value. It can be watched from the Inverse Fourier transform. Then these sin waves are added to each other to create time domain signal. If you change phase relations in your data frame, you change phase relations in sine waves in the sum. By this multiplication you actually do this. So obviously it alter PAPR. OFDM precoding technique (e.g. based on Hadamard orthogonal transform) is an example of using this Fourier property.


Creating a time domain filter from a frequency domain plot

I don't know much about arrays, so I can only answer the second part of your question.

$$A(\omega)=(1+R\cos\omega T_0)\textrm{rect}\left(\frac{\omega}{2\pi B}\right)\tag{1}$$

where $\textrm{rect}()$ is the rectangular function. Equation (1) can be rewritten as

$$A(\omega)=\left(1+\frac{R}{2}e^{j\omega T_0}+\frac{R}{2}e^{-j\omega T_0}\right)\textrm{rect}\left(\frac{\omega}{2\pi B}\right)\tag{2}$$

which shows that $A(\omega)$ is a sum of three ideal low-pass filters with the same cut-off frequency, but with different delays ($0$, $-T_0$, and $T_0$) and different scaling factors. It is an ideal system and can only be approximated by a real system. The implementation suggested by (2) is a parallel connection of 3 low-pass filters with different delays and scaling factors. Note that in order for a causal filter to be able to approximate this ideal system reasonably well, you have to add a delay $\tau>T_0$. So the desired complex frequency response which you can base your design on is

$$H(\omega)=A(\omega)e^{-j\omega\tau}=\left(1+\frac{R}{2}e^{j\omega T_0}+\frac{R}{2}e^{-j\omega T_0}\right)e^{-j\omega\tau}\textrm{rect}\left(\frac{\omega}{2\pi B}\right)$$

Note that each of the 3 low-pass filters is a linear phase filter. So I would suggest you use FIR filters to approximate $H(\omega)$. You first design an FIR low-pass filter with the desired cut-off frequency with group delay $\tau-T_0$ (in samples). For a linear-phase filter, the delay determines the filter length $N=2*$(delay in samples)$+1$. So by choosing $\tau$ you also choose the filter length. If the impulse response of this filter is $h(n)$, and if $n_0$ is the number of samples corresponding to $T_0$, then the total impulse response is given by

$$g(n)=\frac{R}{2}h(n)+h(n-n_0)+\frac{R}{2}h(n-2n_0)$$


Demodulating signal by passing incoming carrier wave through low pass filter as opposed to multiplying by carrier wave?

What you're missing is that multiplying a signal by a sine wave changes the frequency at which the signal appears.  Running a signal through a linear time-invariant filter (what you mean above when you just say "filter") can emphasize or de-emphasize the signal content at certain frequencies, but it cannot move the signal in frequency.
The classic explanation for this is to consider a signal that's a single sine wave at some frequency $f_{RF}$.  So, $x(t) = \cos(2 \pi f_{RF} t)$.  Put it through a low-pass filter and it's just gone.  However, if you multiply that signal by a local oscillator, you generate signals at the sum and difference of the signal frequency and the carrier frequency.  Let your local oscillator signal be $x_{lo}(t) = \cos(2 \pi f_{LO} t)$.
Multiply them together, and you get $$x_{IF}(t) = x_{RF}(t) x_{LO}(t) = \cos(2 \pi f_{RF} t) \cdot \cos(2 \pi f_{LO} t)$$
Now it's just trigonometry: $$x_{IF}(t) = \frac{1}{2}\left(\cos \left(2 \pi f_{RF} - f_{LO} \right) t + \cos \left(2 \pi f_{LO} + f_{RF} \right) t \right)$$
Now you have something useful for a radio receiver: the signal at $f_{RF} - f_{LO}$.  The only problem is that you have the other signal, too.  So now you run that through a low-pass filter -- and the result is your desired signal.


What is the difference between the PSD and the Power Spectrum?

The power spectrum is a general term that describes the distribution of power contained in a signal as a function of frequency. From this perspective, we can have a power spectrum that is defined over a discrete set of frequencies (applicable for infinite length periodic signals) or we can have a power spectrum that is defined as a continuous function of frequency (applicable for infinite length aperiodic signals).

In the first case, each discrete component has units of power (W, mW, etc.). In the second case, each point in the continuous spectrum has units of power per frequency (W/Hz, mW/Hz, etc.). In the second case, one must integrate over a band of frequencies to obtain units of power. It is meaningless to talk about the amount of power at a frequency in the second case, one must talk about the amount of power contained in a spectral band (an interval of frequencies) or in several such bands. It is meaningful however to talk about the amount of power spectral density (PSD) at a particular frequency. This is an indicator of how much weight this frequency will contribute to the overall power if included in one of the spectral bands.

This allows us to compare the distribution of power at various frequencies in a signal (as opposed to comparing the power directly). This distribution is defined over a small delta of frequency (delta going to zero in the calculus sense). An alternative to PSD, would be to divide up the signal into chunks of a finite size (say 10 Hz) and compute the power for each bin (one way to obtain this would be to integrate the PSD over each bin). If another person performed a similar experiment and generated their own plot using a larger bin size (say 20 Hz), then the overall shape of both plots would be similar, but the latter plot would seem coarser and its numerical values would be greater. This would make it difficult to compare the two signals. Using power spectral density, the amplitude and bin size ambiguity is eliminated. This allows a fair comparison between
the power distribution of two different signals.

For more information, here is a link to a related Wikipedia article.


How to change the shape of a function generator keeping the rate?

In music, we would call this "wavetable synthesis", but most electrical engineers would call this Direct Digital Synthesis (DDS) or a Numerically Controlled Oscillator (NCO).

You need a phase accumulator, a sufficient number of points in your wavetable and a means to interpolate between points (if there are many points in the wavetable, linear interpolation is good enough).  And if you're changing the shape of the waveform, you need to be able to cross-fade from the old wavetable to the new.


How should I isolate a single QPSK signal from a noisy band?

I think this should be seen in context with Phil's other, excellent, question:

How can I perform timing recovery from the envelope of PSK signals?

where the situation is the following:

There's 15 equidistant PSK subcarriers, and the center one carries a synchronization sequence, continously.

Channelization, as approached by Phil (and how I'd recommend, too) was done with a polyphase filterbank. Now, synchronization becomes a bit tricky: If using the pulse shape filter in the filter bank, then a lot of the info that a polyphase filterbank synchronizer could use to reconstruct timing information is lost.

So my immediate reaction was: well, then don't use the RRC filter in the filterbank, but use something that passes the whole RRC spectral shape nearly unharmed, thus has a pretty wide passband and pretty steep transition band. Technically, that would've posed no problem – PFBs are very efficient, and the sampling rates were very low.

So, thinking about this, it becomes a Chicken/Egg problem:

If you wanted to use such a "sharp" pre-filter, you'd need to make sure that the RRC-shaped signal is pretty much in the center of it, which requires frequency sync. If it's not centered, you'll cut off substantial parts of the signal's energy and/or introduce adjacent channel interference
Frequency sync is the whole reason we want to do this in the first place.

So, maybe, my idea wasn't that sharp to begin with.

However, this question is a good one, because

it showed that my ideas aren't always the best ones and
this question goes a whole lot farther:

Dropping knowledge of the original signal bandwidth, timing, and spectral location, this becomes a very interesting estimation problem.

So, not to scare you, but: demodulation, compared to sync, is typically the easier part. That's the reason why, for example, the GNU Radio tutorials kind of avoid going deep into that topic. Synchronization often happens by employing specially crafted data in the transmission, or is interwoven with equalization, and hence often very application-specific.

But: Of course there are re-usable, standard, known-to-work-well methods.

Often, one starts by looking at this whole problem to be solved by a combination of three mechanisms:

Rough frequency estimation/correction
Carrier Phase / fine frequency estimation/correction
Timing Recovery

But you have an additional problem: you don't know the analog bandwidth of your signal, meaning that you either don't know the parameters of your pulse shaping filter, or the symbol rate. (We'll leave out Doppler spread for now)

Let us assume we know, however, that we've got to deal with PSK, and that we at least know the number of carriers.

Rough Frequency Correction

A classical approach, if knowing that we only have one carrier in a given bandwidth, would be a band-edgel FLL; simply put: let's have two overlapping band pass filters. If the upper one sees more power than the lower one, the signal isn't centered, increase the correction frequency, otherwise reduce it.

Statistically, this would work with multiple carriers, too; doesn't really matter whether the energy comes from a single carrier, or from 15 that are equidistant in frequency. Now, sadly, in fact, this doesn't work, usually, because one of the main reasons for using multiple subcarriers instead of just one broader channel is that you can then equalize them more easily isolatedly; so you do that in a situation where the channel's attenuation typically isn't flat, and hence, the whole "power equilibrium" hypothesis breaks down.

But, not all is lost: Band-edge FLL, whilst being a standard technology in analog receivers, isn't the only way to determine/correct position of carriers. We now enter the world of frequency estimators, or, given that we know some parameters of the transmission, but not others, parametric frequency estimation.

For example, a simple method would be to look for peaks in the averaged spectrum as given by a DFT of sufficient resolution. Another method would be employing something like ESPRIT, since we know the number of spectral peaks we're expecting.

Given that we now would now the position and spacing of the carriers:

Bandwidth and Pulse Shape Parameter Estimation

This is a pretty interesting problem; with the carrier spacing and position from the previous problem, we can now at least partially isolate one carrier. Assuming sufficiently high SNR, a variance estimate can be used to determine RRC roll-off. Note that this is really the case where you'd need a filter with a flat passband and a steep cutoff – but you'd only need it for one channel, so the computational load might be manageable. By the way, if it works for one channel, and you can assume spectral equidistance, willfully aliasing multiple channels atop of each other might reduce estimation ramdomness – again, this depends on the system you're trying to determine.

Notice that this is no universal approach – you'll need different estimators for different pulse shaping filters.

Timing recovery

You wanted to use the Polyphase Clock Sync timing estimator, which is an excellent choice – it uses a bank of slightly shifted filters to determine the error. You can simply apply that to your frequency adjusted multi-carrier input signal. Just make sure the samples per symbol are appropriate for the sampling rate you'd expect, and choose filter taps of a narrow low pass filter reflecting your estimated channel bandwidth! This might mean you have to interpolate your estimated RRC taps to the full sampling rate, but since this is a decimating operation internally, the computational load might be tolerable, even with the interpolated FIR. Other than that, decimating early in the signal chain is of course the way to minimize computational load.


There is a general question and a specific question.

To answer the specific question, you can guess at an LTI and check if it satisfies your spec:

$$ h[n] = a^n u[n]  \qquad \text{where} \ |a| < 1 $$

Here is my guess:

$$ y[n] = a \cdot y[n-1] + x[n]  \qquad \text{where} \ y[-1]=0 $$

you can set $x[n] = \delta[n]$ and you will find that $y[n] = h[n]$ for all integer $n$ and you are done with the specific question.

To see if a conceptual discrete-time LTI can be made to satisfy a general impulse response spec, $h[n]$, there are a few requisite properties to satisfy:

Is the system causal?  This means the output $y[n]$ does not react to future samples of the input $x[n]$.  And if it's both LTI and causal, that means the impulse response does not react to the input impulse before that input impulse takes on non-zero values.  And that means $$ h[n] = 0 \qquad \forall \ n<0$$  You can see that this is satisfied by your example spec.
Is the LTI system FIR or IIR?  If it's FIR, that means that the non-zero length of the impulse response, $h[n]$ is finite in length (let's call that FIR length "$M$") and that means $$ h[n] = 0 \qquad \forall \ n<0, \ n \ge M$$ Now it turns out that your example is not FIR (which means it is IIR, and Infinite Impulse Response LTI system), but if it were FIR, then your general system can always be realized (or "actualized") as $$ y[n]=\sum\limits_{m=0}^{M-1} h[m] \, x[n-m] $$  That can always be done for an FIR.  It might not be the most efficient way to make an FIR for extremely large $M$, but it is a method that, if you can afford it, will always work.
If it's an IIR (which means the condition in paragraph 2 above is not satisfied), then you have some restrictions about what can be realized and what cannot.  And you have some other conditions to nail down.  One is "BIBO stability" which means, for any bounded input $x[n]$, the output $y[n]$ is also bounded.  Now, I don't know all of the implications, but that does mean that systems that blow up, like $$ h[n] = a^n u[n]  \qquad \text{where} \ +1 < |a| $$ can be realized (same input/output equation as in 2), but they will blow up.  I don't think you want that.  BIBO stability also means that impulse responses like $$ h[n] = \tfrac{1}{n+1} u[n] $$ will also blow up because $\tfrac{1}{n+1}$ does not converge fast enough.  If, say, $x[n]=u[n]$ (which is bounded), the output $y[n]$ will blow up.
So it turns out that BIBO stable LTI will need to either be FIR (they're always stable from an input/output POV) or, if the LTI is IIR, then the impulse response must decay exponentially.  That means, for large $n$, $$  \Big|h[n] \Big| \le A \ a^n \ u[n] \qquad \text{for some} \ |a| < 1  $$  If that's true, you can conceptually make an IIR filter that approximates any $h[n]$ that satisfies the causality condition (paragraph 1) and the BIBO stability conditions (paragraph 3).  The way we learn to do that is by taking a DSP course after learning LTI system theory.


Python - Normalized cross-correlation to measure similarites in 2 images

I guess you can compute for each pixel the correlation coefficient between patches centered on this pixel in the two images of interest. Here is an example where I downloaded the figure attached here and tried to compute the correlation in such a way. The output looks different from the one of the article, but it was to be expected since the resolution is very different.

from skimage import io, feature
from scipy import ndimage
import numpy as np

def correlation_coefficient(patch1, patch2):
    product = np.mean((patch1 - patch1.mean()) * (patch2 - patch2.mean()))
    stds = patch1.std() * patch2.std()
    if stds == 0:
        return 0
    else:
        product /= stds
        return product

im = io.imread('faces.jpg', as_grey=True)

im1 = im[16:263, 4:146]
sh_row, sh_col = im1.shape
im2 = im[16:263, 155:155+sh_col]

# Registration of the two images
translation = feature.register_translation(im1, im2, upsample_factor=10)[0]
im2_register = ndimage.shift(im2, translation)

d = 1

correlation = np.zeros_like(im1)

for i in range(d, sh_row - (d + 1)):
    for j in range(d, sh_col - (d + 1)):
        correlation[i, j] = correlation_coefficient(im1[i - d: i + d + 1,
                                                        j - d: j + d + 1],
                                                    im2[i - d: i + d + 1,
                                                        j - d: j + d + 1])

io.imshow(correlation, cmap='gray')
io.show()

The code above is a naive and slow implementation of the correlation, as the two for loops are very slow. For faster execution, you could for example port the script to Cython.

In the article, I think the idea is to measure whether face expressions look similar or not. If a pixel has a large correlation index between two images, it means that the region of the face where this pixel is located does not change much between the images.

If you use this method on good-resolution images, you should increase the patch size for more accurate results (d=2 or 3).


Two generated tones creates a side tone - how to avoid the side tone

I have tried the following:

Launch Audacity.
Generate a 15000 Hz tone in the track created by default.
Add a new track.
Generate a 15400 Hz tone in the new track.

A lower frequency tone appears during playback.

The reason is that both tracks have high levels, so their sum exceeds 1.0; and Audacity applies clipping or limiting. This non-linear operation is what creates the overtone.

Simply decrease the level of both tracks to -6dB (or use an amplitude of 0.5 during signal generation) to make the overtone disappear.


One easy way to grab a set of samples of physical noise is by acquring a cheap SDR, such as an RTL, tuning it to an empty band, and just reading and saving the samples. This should give a you a nice set of mostly white noise samples.


fourier transform of smart-phone accelerometer in matlab

Two methods occur to me.

Perform a non-uniform DFT. You can use Google to find MATLAB implementations of it, or try to write one on your own. This file seems to be legit, although I haven't tried it myself.
If you have the Signal Processing Toolbox, then you can use the resample() function to interpolate your non-uniformly sampled vector linearly onto a vector of uniformly spaced instants, and then perform the usual DFT on the latter.


Computer Vision - Does implementation of Hough Transform for lines require a loop for angle parameter?

In real world implementations, is the line: for theta = [theta_min to theta_max ] required?
. . .
Since we apply Canny edge detector to an image, prior to Hough Transform, we can obtain angle of gradient(θ) for a point (x, y) from Sobel operators. This can be noted down against (x, y) and later used in Hough Transform algorithm, eliminating one for loop. Isn't this the better, real world approach?

The standard Hough Transform algorithm does not assume that direction information is available for the image that is transformed. It basically accepts a set of point values (pixel location, pixel value) ($(x,y,v)$) at its input and converts lines to points. This is because, colinear points 9lines) contribute to large "accumulations" (counts) along the directions of the line they form.
So, a pixel of value $v$ at some location $x,y$ contributes $+v$ to all the Hough domain arcs that pass through it.
Now, if we consider the line formed by two neighbouring pixels $(x_0,y_0,v_0), (x_1, y_1,v_1)$ that both happen to have values $v_0=v_1=v$ (for simplicity), then this line segment's contribution along the directions that are parallel to the line is $2 \cdot v$ but along the vertical direction is still $v$ for two arcs.
As the length ($l$) of the line segment that is formed by any two $(x_m,y_m), (x_n,y_n), m \neq n$ increases, even if the pixels are not neighbours now, then they contribute $l \cdot v$ to that one direction of convergence but still $v$ to $l$ arcs in the vertical direction. But also, as that length increases, it starts contributing between $v .. l \cdot v$ to a wider range of $\theta$s around the direction that is parallel to the line.
Therefore, the answer to the first question is "yes". The answer to the second question is "under certain conditions".
This is because, you can use some of the information of the Canny edge detector to inform the accumulation process of the Hough Transform, but you have to do it carefully.
The Canny edge detector implies a scale parameter because it uses a Gaussian filter. This scale parameter, along with the gradient information can be used to further subdivide the Canny edge detector's output into a set of composing elementary line segments at direction $\theta$ and length $l$. Then, the Hough Transform would be approximated by adding together the individual Hough Transforms of those elementary lines. But of course, notice here that you still have to "scan" the elementary line segments at a range of $\theta$s to obtain their contribution during this pre-computation stage.
Hope this helps.


Is Gabor uncertainty a feature of the Fourier transform, or of nature?

Time/frequency uncertainty is a property of the meanings or definitions attached to the terms "time" localization and "frequency" estimation.  It also alludes to some apriori knowledge, or lack of it, about the characteristics of the signal in question.


You're right that the angle $\Omega(\omega)$ in your sketch will equal $\pi$ for $\omega=\pi$. However, $\Omega(\omega)$ is not the phase angle of the system. In order to see this, let's write down the frequency response of a system with a single zero at $z_0=re^{j\theta}$:

$$H(e^{j\omega})=1-re^{j\theta}e^{-j\omega}=\frac{e^{j\omega}-re^{j\theta}}{e^{j\omega}}\tag{1}$$

From $(1)$, the phase angle is given by

$$\angle\{H(e^{j\omega})\}=\angle\{e^{j\omega}-re^{j\theta}\}-\angle\{e^{j\omega}\}=\Omega-\omega\tag{2}$$

So you indeed need the angle $\Omega$ as defined in your figures, but you also need to subtract the angular frequency $\omega$.

Let's assume $\theta=\pi$, as in your example. For $\omega=\pi$ (and $r<1$) you get $\Omega=\pi$, and, consequently, a phase of $\Omega-\omega=0$. Note that if $r>1$, you'll get $\Omega=0$ at $\omega=\pi$, and, consequently, a phase angle of $(-)\pi$.


Main differences to take into account between continuous and discrete time signals

In general I'd say that you understand the differences pretty well.

You nailed the first difference- continuous time vs. discrete time (samples).
You were also correct about there being sets of indistinguishable frequencies in the discrete realm.  A corollary of this is that all discrete signals are band-limited, whereas continuous signals do not have to be.  This concept is most often referred to as the Nyquist frequency.
I don't entirely agree with your point about periodicity, though I believe I understand what you're getting at.  If by periodicity you mean having the exact same values repeated, then yes, a discrete period would have to be an integer number of samples.  On the other hand, if you take periodic to mean that the "underlying" signal repeats, then it is possible to have the samples not be the exact same if the period is not an integer multiple of samples.  You could still detect the periodicity with an autocorrelation though.  Anyway, if this point is confusing, don't worry about it- it's not terribly important.

It is, by the way, possible to view discrete signals as a special case of continuous signals.  You can model discrete signals as continuous signals that have been multiplied by an infinite set of evenly spaced delta functions.  If the delta functions are separated by $T_s$ seconds, then in the frequency domain they can be represented by another set of delta functions that are separated by $\frac{1}{T_s}$ Hz.  The multiplication in the time domain is equivalent to convolution in the frequency domain, so copies of the signal are reproduced in the frequency domain at each of the delta functions.

Though the mathematicians in the audience will probably scream at me for saying this, I think that you can likewise look at continuous signals as a special case of discrete signals, where we take the limit of $T_s$ as it approaches zero.


Determining how similar a signal is to a sine wave

Assuming that the description I made of the FreezeFramer process is correct, then I would try determining the SINAD rating of the peak frequency in the interpolated data.

You have a sampling rate of 2Hz, giving a Nyquist of 1Hz.  You will be looking for a sinusoid with a frequency well below 1Hz.  The SINAD value is (basically) the amplitude of the peak frequency divided by the sum of all the other frequency amplitudes.  There are better descriptions available, as well as implementations for many of the more common DSP frameworks.


Help in deriving the update rule for equalizer

The mapping from $z_n$ to $s_n$ is not bijective, so it can't be inverted. This means that even in the transmitter, $z_n$ cannot be recovered from $s_n$. This is why the constant modulus algorithm (CMA) need not (and cannot) be modified to estimate $z_n$. It can only be used to estimate $s_n$, and for this purpose no modification is necessary.

The thesis you referred to in a comment is mainly about blind signal separation. The update rules (3.16) and (3.17) apply to the iterative maximization of the log-likelihood function in order to compute the maximum likelihood solution of a source separation problem. The nonlinear vector-valued function $\mathbf{g}(\mathbf{u})$ depends on the probability density functions of the different source signals. The underlying model and the problem formulation are different from the model and the problem underlying the CMA, so there is no direct relation between the two, and the update rules mentioned above are not applicable to the CMA.


Pole Zero plot given a Transfer function

What you have are not the poles and zeros, but simply the filter coefficients, i.e., the coefficients of the numerator and denominator polynomials.

The poles are the roots of the denominator polynomial, and the zeros are the roots of the numerator polynomial. In Matlab they can be found by using the roots command:

p = roots(a);
z = roots(b);

Note that in general, poles and zeros are complex numbers, that's why they are plotted in the complex plane.

Just a remark: you used the zplane command with numerator and denominator interchanged, that's why the plot shows the zeros as crosses on the unit circle, and the poles as 'O's inside the circle. The correct call is zplane(b,a).


Why could we use DoG? An old and low-cost edge-enhancement technique consists in subtracting a blurred image from the original image. Its filter interpretation  is an impulse filter (the neutral filter for the original image) minus the blur filter. The resulting filter is positive at its center, and negative around, a sort of coarse Laplacian filter.

The above impulse filter can be seen as the limit of a Gaussian filter whose $\sigma$ tends to $0$.

The DoG filter somehow generalizes the described old edge-enhancement technique, with a pair of blurring filters. Perhaps surprisingly, two blurs, when subtracted, can sharpen edges. The resulting filter acts as a pass-band with appropriate $\sigma_k$ kernel ratios.

How effective is DoG? its effectiveness is close to that of the Laplacian of Gaussians, which it mimics, see below in 1D.

Which type of images work best for this filter? It is believed to work well for large images, since it has very fast implementations. Due to the "double filtering effect" in high frequencies, it is somewhat robust to noise. Depending on the kernel size choice, it enhances features that possess a certain scale extend.


What is the meaning of **Sample** in **Sample ACF**

It's the empirical ACF computed using the sample.


Prove the following DFT for a periodic sinc sequence or ideal lowpass filter

Hints:

1) Definition of the unnormalized DFT:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i\frac{2\pi}{N}kn}  $$

2) Geometric Summation Formula:

$$ \sum_{n=0}^{N-1} r^n = \frac{1-r^N}{1-r} $$

3) Exponential form of Sine:

$$ \sin( \theta ) = \frac{ e^{i\theta} - e^{-i\theta} }{ 2i } $$

4) Definition of unnormalized sinc:

$$ \operatorname{sinc}(x) = \frac{\sin(x)}{x} $$


Upsampling and Downsampling using IFFT and FFT (DFT Based Resampling)

Given a time domain signal $x_1[n]$ of length $N$ and its $N$-point DFT $X_1[k]$, If you pad $M$ zeros to the end of $x_1[n]$ and then take ($N+M$) point-DFT of it, this will give you an interpolated set of frequency samples, $X_2[k]$; i.e., a denser set of uniform samples of the DTFT, $X(e^{j\omega})$, than was previously available from $X_1[k]$.
Given a time domain signal $x_1[n]$ of length $N$ and its $N$-point DFT $X_1[k]$, If you carefully stuff (insert) $M$ zeros to the central portion of the DFT $X_1[k]$ (by moving those right half of the available $X_1[k]$ samples to the right by the amount of $M$) and then take ($N+M$) point-inverse DFT of the resulting $X_2[k]$, this will give you an interpolated set of time samples, $x_2[n]$; Note that you should scale the amplitudes in accordance with the interpoaltion amount.

Note that implementing 1-) is straightforward, while implementing 2-) requires careful selection of the correct indices at the central portion for the original length $N$ being even and odd, the details of which is provided by some examples here dsp.se,  including the link you have provided.


Quadratic Programming with Linear Equality Constraints

Let's solve a more general problem (Least Squares with Linear Equality Constraints):

$$
\begin{alignat*}{3}
\arg \min_{x} & \quad & \frac{1}{2} \left\| A x - b \right\|_{2}^{2} \\
\text{subject to} & \quad & C x = d
\end{alignat*}
$$

The Lagrangian is given by:

$$ L \left( x, \nu \right) = \frac{1}{2} \left\| A x - b \right\|_{2}^{2} + {\nu}^{T} \left( C x - d \right) $$

From KKT Conditions the optimal values of $ \hat{x}, \hat{\nu} $ obeys:

$$ \begin{bmatrix}
{A}^{T} A & {C}^{T} \\
C & 0
\end{bmatrix} \begin{bmatrix}
\hat{x} \\
\hat{\nu}
\end{bmatrix} = \begin{bmatrix}
{A}^{T} b \\
d
\end{bmatrix} $$

The trick here is to have a look at:

$$ \frac{1}{2} \left\| A x - b \right\|_{2}^{2} = \frac{1}{2} {x}^{T} {A}^{T} A x - {x}^{T} A b + \frac{1}{2} {b}^{T} b $$

So if we set $ w = x $, $ b = \boldsymbol{0} $, $ X = C $, $ d = \boldsymbol{1} $ and $ \frac{1}{2} {A}^{T} A = R $ then your solution is given by:

$$ \begin{bmatrix}
2 R & {X}^{T} \\
X & 0
\end{bmatrix} \begin{bmatrix}
\hat{w} \\
\hat{\nu}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{0} \\
\boldsymbol{1}
\end{bmatrix} $$

This is easy to solve in MATLAB or Python.

Handling Large System

As a request from @user5045 I add some info how to handle this in case $ R $ and $ X $ are large matrices.

We basically need to solve large scale matrix equation:

$$ \begin{bmatrix}
2 R & {X}^{T} \\
X & 0
\end{bmatrix} \begin{bmatrix}
\hat{w} \\
\hat{\nu}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{0} \\
\boldsymbol{1}
\end{bmatrix} = e = F g $$

The way to solve it is using an iterative solver.
I case $ R $ is a PSD matrix then the solver should be Preconditioned Conjugate Gradient. In MATLAB it is implemented using pcg().
In case $ R $ is only symmetric one should use Minimum Residual Solver. In MATLAB it is implemented as minres(). For a a general square matrix I'd go with cgs().


In general it's only possible to implement causal and stable filters. There are exceptions where marginally stable filters are used, but this doesn't apply here. So if you want to invert a given filter, this is only possible if its zeros are inside the unit circle of the complex plane (such filters are called minimum-phase filters). The zeros of the given filter become the poles of the inverse filter, so if the zeros are inside the unit circle, the poles of the inverse filter will also be inside the unit circle, and, consequently, the resulting filter is causal and stable.

The problem with the filter in your example is that it has a sixth order zero at $z=-1$, i.e. on the unit circle. Consequently, this filter is not strictly minimum-phase and it can't be inverted by a causal and stable filter. However, it can be inverted approximately by a causal and stable filter. One way to achieve such an approximate inversion is by minimizing the mean squared error between the ideal and the actual filter response. In practice, it is often acceptable to add some delay to the output, which will make it easier to approximate the desired response by a causal and stable filter.

Note that transients don't play a role. It's the inverse filter's instability that is the problem here.

As an example, take a filter with coefficients

b = [1,0,-0.5]; a = [1,0.3,0.1];

The zeros of the numerator polynomial (with coefficients b) are all inside the unit circle, so the inverse filter is causal and stable.

x = rand(100,1);
y = filter(b,a,x);
x2 = filter(a,b,y);
max(abs(x2-x))    % in the order of 1e-16, i.e., machine precision


It sounds like you are trying to perform something between MP3-to-MIDI or music to score sheet transcription. As far as I know, this both remains

an active research topic, with emphasis in robustness (to noise, instrument) and real-time, see for instance:

Automatic Music Transcription: Challenges and Future Directions, 2013
Robust Real-Time Music Transcription with a Compositional Hierarchical Model, 2017
Learning Features of Music from Scratch, 2017

a handful of existing software, like

Notation Software, Converting Audio files to Sheet Music
ScoreCloud Studio, Music notation software for everyone
AnthemScore, Automatic Music Transcription Software


OK, I need to answer my own question.
The math transform I described is GOOD, the problem was like always, stupid mistake in the code. But Stanleys answer helped me much. Thanks a lot.


Reconstruct DWT for each cD1,cD2,cD3 and cA3 signals

Alright, I'll try to take a shot at giving you an explanation of what's going on with the DWT.

So in the CWT, basically what you are doing is you are generating the wavelet coefficients by convolving the signal with each scale and shift of the mother wavelet function. However, the problem with doing this is that when you analyze the low frequency components of the signal, you are oversampling, or performing the wavelet analysis too often. Also, the relationship between time and frequency localization is much like the relationship between velocity and position localization in quantum mechanics, and when you perform the CWT this inequality is violated for some time/frequency scales. Essentially, with regard to the CWT, the transform is redundant at low frequency scales.

However the DWT doesn't have this problem because the DWT effectively downsamples the signal at each level. The DWT is basically just two filters, a highpass and a lowpass filter. The highpass filter produces the "detail" coefficients, and the lowpass filter produces the "average" coefficients. However, keep in mind that calculating the moving average of a signal is a form of lowpass filtering. Thus when we perform the lowpass filter and get the "average" DWT coefficients, those coefficients are basically just the low frequency components of the signal. And, since we know they are lower in frequency, we know that in the digital case we need less data points to represent a signal lower frequencies. The lowpass and highpass filtering is then repeated for each level of the DWT that can be computed. The level is dependent on the length of the signal. If a signal has $$k$$ data points, where $$k = 2^N$$
then the DWT can be computed to N levels.

I'll try to give a real life example. Lets say you have a signal of 1024 samples, sampled at a rate of $$4\pi \text{ Hz.}$$ Based on the Nyquist criterion, the highest frequency that can be resolved by the signal is $$2\pi \text{ Hz.}$$Computing the first level of the DWT will result in two vectors, each of length 512. One is from the highpass filter, the other from the lowpass (these vectors are sometimes called high/lowpass filter banks). The highpass coefficients are the cD1 coefficients, representing frequencies in the range of one to two pi, and the lowpass coefficients are the cA1 coefficients representing frequencies in the range of zero to one pi. The DWT is applied again to the lowpass coefficients resulting in two vectors each of length 256, and this time the lowpass filter bank represents frequencies in the range 0 to 1/2 pi. This process is repeated 10 times (for this case, 2^10 = 1024), as the last time will result in lowpass and highpass filter banks each of length 1.

Basically at higher scales of the DWT, we get better frequency resolution, but poorer time resolution. Lower scales have better time resolution but poorer frequency resolution. I'll also give you the link to a good tutorial on this stuff by Polikar : http://person.hst.aau.dk/enk/ST8/wavelet_tutotial.pdf

To reconstruct the signal, simply apply the inverse DWT for the number of levels which you want to reconstruct.

EDIT: Forgot to mention that the FFT is not needed, you are already filtering based on frequencies.


Why would a time shifted signal have a higher maximum frequency than an unshifted one?

I was able to take a look at the book.
You are right it says: "The highest significant frequency in the shifted Gaussian signal is higher than that in the centered Gaussian".  I think this is misleading.
Looking at Figs 2.2 and 2.3 I can only think that the author is saying that in Fig 2.3 the phase is "significant" across the whole spectrum, while in Fig 2.2 phase is zero everywhere.
However, the magnitude curves are identical. If magnitude is zero then phase is irrelevant. (i.e. atan (small/small)).  The required sample interval is governed by the magnitude spectrum. That is the required sample interval is not affected by a simple time shift, as stated by the above commenters. If this wasn't the case, then real-data analog-to-digital conversion would be pretty difficult!


The equations are a bit wrong. For a signal $x=cos(2\pi f_{c}t+n\pi), n\in \left\{0,1\right\}$ transmitted through a complex exponential channel $h=Ae^{i\theta}$, the received signal $y$ is not $hx$. $h$ is in the complex baseband notation and your signal $x$ is in the bandpass notation.

To answer your actual question first, $\theta$ is the phase change that is incurred by the signal when it passes through the channel i.e. the delay. And the actual received signal is $y=Re\left\{x_{l}Ae^{i\theta}e^{i2\pi f_ct}\right\}$, where $x_{l}$ is the baseband representation of our signal $x$ which in our case would be something like $x_{l} = e^{i2\pi\frac{(n-1)}{2}}$, where $n\in \left\{0,1\right\}$.


Obtaining an IR from a sine sweep without the dry sweep

If this is used for channel sounding then you aren't concerned with the actual accuracy of your turntable (causing a frequency shift) in which case I suggest measuring the response in a static condition (non-changing channel) and determining the least square linear fit of the frequency sweep (I would assume but don't know that the test pattern would be a linear sweep) and then from that reconstruct a reference sweep to then use for channel estimation.
Measuring the frequency involves using some for of frequency discriminator process, most of which are also sensitive to amplitude variation, so I recommend hard-limiting the signal prior to the discriminator. Also pay attention to the relationship between instantaneous frequency and the underlying frequency of your ramp as detailed in this "DSP Quiz" question:  Simulation of a Frequency ramp


Removing the Multiples of Decrease and Increase in a Spike

You can use a moving average filter followed by your threshold detector. You can change the span of your moving average filter based on the time length in which you want multiple spikes to be combined.


So, the algorithm I'd propose is the following:

Filter your input signal with a 200 Hz – 2 kHz bandpass filter
Since in that band your sine is the dominant oscillation, a PLL would lock onto that. So, use a second-order PLL to track that wave; there's many PLL implementations that'll output the reference as tone with the correct phase.
Typically, PLLs lose amplitude information, so you can't simply remove the sinusoid by subtracting it from the audio signal. Instead, you'd use the reference tone to multiply the signal with, which mixes down the original tone to 0 Hz. A low-pass filter (in this case, even a single-pole IIR) might be trivially used to extract the instantaneous amplitude.
Now that you have the right tone with the right frequency and phase, you can subtract it from the input
To predict how the tone is going to be developing its frequency, you could use a linear predictor, but honestly, a second-order PLL is a linear predictor for phase (a linearly changing phase is a sinusoid; the coefficient of the linear term is the frequency). If you need to work with a sinusoid which might be undergoing a linear change in frequency: third order PLL.

Even more elegant: Use the output of 2. Mix your signal down with this. Apply 0 Hz Notch filter (again, single-pole IIR might do, or other kind of narrow high-pass filter). Shift up in frequency again.


How to calculate an azimuth of a sound source in stereo signal

So here's a diagram:

The signal source is at $(x, y)$, left microphone is at $(-1, 0)$, right microphone is at $(+1, 0)$.  $d=2$ is the distance between the microphones.  $b$ is the distance from the source to the right microphone, and $a$ is the extra distance to the left microphone, which will create additional delay and further decrease the amplitude. The grey circle represents the wavefront as it hits the right microphone.

The delay (time difference) between channels will be $T = \frac a c$, where $c$ is the speed of sound.
Since microphones are pressure sensors, and pressure drops with inverse distance ($1/r$ law), the signal amplitude at the left microphone will be $G = \frac{b} {a+b}$ times the signal at the right microphone.  (It's a level ratio, not a level difference.)

If you know $T$ and $G$ then you can find $b$:
$$b=\frac{G T c}{1-G}$$

Then you know $b$ and $a+b$ and using law of cosines and Apollonius' theorem for the median of the triangle you can find $r$ and $\theta$:

$$r = \frac 1 2 \sqrt{4b^2+4ab+2a^2-d^2}$$
$$b^2=-d r \cos(\theta) + r^2 + \frac{d^2}{4}$$
$$\theta=\arccos\left( \frac{r}{d}+\frac{d}{4r}-\frac{b^2}{dr}\right)$$

Then plug in $a$ and $b$ and simplify:

$$r = \frac 1 2 \sqrt{\frac{2(Tc)^2(G^2+1)}{(G-1)^2} -d^2}$$
$$\theta=\arccos\left( \frac{r}{d}+\frac{d}{4r}-\frac{(G T c)^2}{(1-G)^2 d r}\right)$$

If you only know $T$, you can place the source along a hyperbola (shown in green) with focal points at the microphones.  If you only know $G$, then you can place the source along a circle (shown in orange).  The intersection of those curves is the location of the source:

I would think that accuracy becomes smeared out when the green and orange curves meet at sharp angles.

Of course if you only care about separating sounds by location, you can ignore all this and just work in (delay, ratio)-space.


How do you deconvolve irregularly sampled points?

You can use a graph-based Signal Processing framework to apply regular deconvolution algorithms. The whole point is that traditional operators (such as FFT) can be defined on graph structures, and graph structures do handle "by design" without any overhead irregular sampling patterns.

Some useful reference scan be found on my former office mate's webpage.


Fractional/Interpolating Delay Line still sounding glitchy

If you change the delay parameter suddenly, there will be a jump to another point in time in the data, resulting in a glitch. The easiest is to filter the delay parameter with a one-pole low-pass filter:

$$\text{out}[k] = \text{out}[k-1] + a(\text{in}[k] - \text{out}[k-1]),$$

where $a$ is a constant between 0 and 1, $\text{in}[k]$ is the delay parameter controlled by the user and $\text{out}[k]$ is the delay parameter actually used in the delay line, for integer time index $k$. ($\text{out}[k-1]$ is the previous filtered value.) You can cascade these stages in series to get further smoothing. Two stages should sound good because it will get rid of discontinuities in the pitch shift due to a kind of doppler effect that changing the delay length creates.

The constant $a$ can be tuned manually. A good value might be the highest value that does not produce disturbing glitches. For a better explanation of $a$, useful if you need to adapt to a different sampling frequency:

$$a \approx (1-p)^{\frac{f_s}{t_p}},$$

where $p$ between 0 and 1 represents the percentage that the filtered delay parameter has come towards a new value of the input delay parameter in time $t_p$, and $f_s$ is the sampling frequency. The formula is for a single filtering stage. If you have multiple stages the meaning of $t_p$ changes but you can still use it to go from one sampling frequency to another. The -3 dB cutoff frequency $f_\text{cutoff}$ for a stage will be approximated by:

$$f_\text{cutoff} \approx \frac{-f_s\ln a}{2\pi} $$


On the stability and causality of a discrete system

To your first question, a linear system is causal if we have initial conditions of rest that is zero. This is a sufficient and necceassy condition. This property doesn't  require time invariance

Stability is a consequence of absolute summability of system response. It doesn't depend on the initial condition (of course have to be finite to make any sense).

To your second question the transfer function has a pole at z=3 and it is given that the system is stable. That means the ROC has to include the unit circle. That means that the ROC is to the left of the pole at z = 3, since ROC is towards the left of outermost pole hence system is not causal


The longer the total sample window (T0), the smaller the frequency difference between two sinusoids with exactly M periods and exactly M+1 periods (for integer M <= N/2) in that window, which is required for orthogonality between adjacent DFT bins.

Also, the bin spacing is not the resolution. The frequency resolution can be greater or less than the DFT bin spacing, depending on the S/N ratio and your resolution needs (separation or estimation).


Different results for different orders of estimating AR model using Yule-Walker equations

Your true AR model is order 4. If an AR(16) model did a good job, then 12 of its parameters should be zero. Because of noise effects / overfitting, those 12 parameters are NOT zero and, more than likely, the other matching parameters will be more inaccurate because of this. Hence the error in an AR(16) will be larger than that for an AR(4) fit.


How do we compute distrubtions of the value of a random process conditional on initial conditions?

"I know that $\phi(t)$ is Gaussian distributed" is not the same as saying
"$\{\phi(t)\}$ is a Gaussian process" but I will assume that the latter is meant. With this assumption, and the additional assumption that the process is
wide-sense-stationary (the autocorrelation function is listed as having only one argument), the process is strictly stationary and thus
$\phi(t)$ and $\phi(0)$ are jointly Gaussian random variables. They have
the same mean $\mu$ and variance $\sigma^2$, and
their correlation coefficient is $\rho$ where
$$\mu = \sqrt{\lim_{t\to \infty} A(t)}, \qquad \sigma^2 = A(0)-\mu^2,
\qquad \rho(t) = \frac{A(t)-\mu^2}{\sigma^2}.$$
Consequently the conditional distribution of $\phi(t)$ given
that $\phi(0)$ has taken on value $\phi_0$ is Gaussian
with mean $\mu + \rho(t)(\phi_0 - \mu) = \rho(t) \phi_0 + (1-\rho(t))\mu$
and variance $\sigma^2(1-\rho(t)^2)$. Note that as $t\to\infty$,
the conditional distribution of $\phi(t)$ approaches the
unconditional distribution of $\phi(t)$: the distant past affects
the present state less and less.

For more information about Gaussian processes, see, for example,
this answer to a
different question.


binary matrices with directional non-zero elements

If you assume that your image origins from sampling a continuous-spatial-domain 2-dimensional function, you can perform the sampling of this function with the rotated coordinates. HAve a look at this example:

N = 32  # 32 for simplicity
N2 = N/2
x, y = np.meshgrid(np.arange(N), np.arange(N))

# exemplary rotation around 5 degrees
angle = 5*np.pi/180
# perform rotation around the image center
nx = np.cos(angle)*(x-N2) - np.sin(angle)*(y-N2) + N2
ny = np.sin(angle)*(x-N2) + np.cos(angle)*(y-N2) + N2

# three functions that fulfill your constraint of even-odd 1-0
I_func1 = lambda x, y: 0.5*(1-np.cos(np.pi*y))
I_func2 = lambda x, y: (abs((y%2)-1)<=0.5).astype(float)
I_func3 = lambda x, y: (1-abs((y%2)-1)).astype(float)

# the original image
I = I_func1(x, y)

plt.figure()
t = np.linspace(0, N-1, 16*N)
t2 = np.arange(N)
plt.plot(t, I_func1(0, t))
plt.plot(t, I_func2(0, t))
plt.plot(t, I_func3(0,t))
plt.stem(t2, I_func(0, t2))
plt.xlim((0, 5))
plt.title('underlying continuous function')

plt.figure(figsize=(12,12))
plt.subplot(221)
plt.imshow(I, interpolation='none', cmap='gray')

plt.subplot(222)
plt.imshow(I_func1(nx, ny), interpolation='none', cmap='gray')
plt.title("Cos")

plt.subplot(223)
plt.imshow(I_func2(nx, ny), interpolation='none', cmap='gray')
plt.title("rect")

plt.subplot(224)
plt.imshow(I_func3(nx, ny), interpolation='none', cmap='gray')
plt.title("triang")

Output:

As you see, different functions yield different rotation results. Maybe this can help you to find an explicity formula?


Should we use a digital low pass or high pass filter to remove sensor bias?

If we use a discrete derivative with a discrete integral in series, will this remove (atleast some of) the DC bias in the signal?

Yes (if your integrator starts at state 0).
Note that this is a CIC filter.
Also notice an offset should be present in your accelerator data, or you just got lost in space.

Correspondingly it seems a better idea to use a high pass filter rather than a low pass filter although this will invite noise.

Well, seeing that by definition of "low-pass filter", the DC component remains untouched, this doesn't only seem to be a good idea, but mandatory.

Is a notch filter how people deal with accelerometer signals and should one go with a high pass approach in case of the unavailability of a notch filter (since the signal necessarily contains bias)?

It depends on what people are doing for what purpose. My guess is you want to estimate something from your accelerometer observation. Maybe Kalman filters are what you're after?

Does the same logic hold if one designs the filter in the continuous domain?

This has nothing to do with continuous or discrete time. So, yes.


Covariance between real and imaginary parts of Fourier transform of a stationary time series

Note that in general the Fourier transform of a stationary process $x(t)$ does not exist. The Wiener-Khinchin theorem only states that under certain conditions the power spectral density of $x(t)$ exists, and it can be computed as the Fourier transform of the autocorrelation function of $x(t)$.

Having said that, if for some reason one assumes that the Fourier transform of $x(t)$ exists, then you can do the math and see if you get a useful result. So, let's see. I assume that $x(t)$ is real-valued. This is not at all necessary, but just simplifies things a bit.

Since $\text{Re}\{X(\omega)\}=\frac12[X(\omega)+X^*(\omega)]$ and $\text{Im}\{X(\omega)\}=\frac{1}{2j}[X(\omega)-X^*(\omega)]$, we can compute the desired expectations from $E\{X(\omega)X(\omega')\}$ and $E\{X^*(\omega)X(\omega')\}$. For the first of these expectations we have

$$E\{X(\omega)X(\omega')\}=E\left\{\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt
\int_{-\infty}^{\infty}x(t')e^{-j\omega' t'}dt'  \right\}=\\
=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}E\{x(t)x(t')\}e^{-j\omega t}e^{-j\omega't'}dtdt'$$

With the substution $\tau=t'-t$ and with the autocorrelation function $R_x(\tau)=E\{x(t)x(t+\tau)\}$ (assuming wide-sense stationarity of $x(t)$) we obtain

$$E\{X(\omega)X(\omega')\}=\int_{-\infty}^{\infty}R_x(\tau)e^{-j\omega'\tau}d\tau\int_{-\infty}^{\infty}e^{-j(\omega+\omega')t}dt$$

Interpreting the second integral as a distribution, and with the power spectrum $S_x(\omega)=\mathcal{F}\{R_x(\tau)\}$ we finally get

$$E\{X(\omega)X(\omega')\}=2\pi\delta(\omega+\omega')S_x(\omega')=
2\pi\delta(\omega+\omega')S_x(-\omega)=2\pi\delta(\omega+\omega')S_x(\omega)\tag{1}$$

because $S_x(\omega)$ is an even function (we assumed $x(t)$ to be real-valued). The other expectation $E\{X^*(\omega)X(\omega')\}$ can be derived in a completely analogous manner. The result is

$$E\{X^*(\omega)X(\omega')\}=2\pi\delta(\omega-\omega')S_x(\omega)\tag{2}$$

From (1) and (2) we can obtain the desired expectations

$$E\{\text{Re}\{X(\omega)\}\text{Re}\{X(\omega')\}\}=
\pi S_x(\omega)[\delta(\omega-\omega')+\delta(\omega+\omega')]\\
E\{\text{Im}\{X(\omega)\}\text{Im}\{X(\omega')\}\}=
\pi S_x(\omega)[\delta(\omega-\omega')-\delta(\omega+\omega')]$$

It remains to compute the quantities $E\{\text{Re}\{X(\omega)\}$ and $E\{\text{Im}\{X(\omega)\}$. We can easily derive them from the expectations $E\{X(\omega)\}$ and $E\{X^*(\omega)\}$:

$$E\{X(\omega)\}=E\left\{ \int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt \right\}=
\int_{-\infty}^{\infty}E\{x(t)\}e^{-j\omega t}dt=
\mu_x\int_{-\infty}^{\infty}e^{-j\omega t}dt=2\pi\mu_x\delta(\omega)$$

with $\mu_x=E\{x(t)\}$, where we have again assumed wide-sense stationarity of $x(t)$. Obviously, we get the same result for $E\{X^*(\omega)\}$:

$$E\{X^*(\omega)\}=2\pi\mu_x\delta(\omega)$$

This results in

$$E\{\text{Re}\{X(\omega)\}\}=2\pi\mu_x\delta(\omega)\quad\text{and}\quad
E\{\text{Im}\{X(\omega)\}\}=0$$

So we have

$$E\{\text{Re}\{X(\omega)\}\}\cdot E\{\text{Re}\{X(\omega')\}\}=4\pi^2\mu_x^2\delta(\omega)\delta(\omega')\\
E\{\text{Im}\{X(\omega)\}\}\cdot E\{\text{Im}\{X(\omega')\}\}=0\tag{3}$$

Combining (3) with (1) and (2) you obtain the desired result. From the result you can see that the covariances of the real parts of $X(\omega)$ and $X(\omega')$ vanish everywhere except for $\omega=\pm\omega'$. The same is true for the covariances of the imaginary parts. The difference is an additional term for $\omega=\omega'=0$ for the covariance of the real parts.


You can achieve this result by using two combs filters : https://en.wikipedia.org/wiki/Comb_filter

Put simply, the comb filter consists of adding a delayed version of the signal to itself, causing destructive or constructive interference.

For instance, with $K = 20$ and a negative gain value after the delay line, you can significantly decrease or suppress harmonics that are even multiples of 20: $\sin(2\pi40t)$ and $\sin(2\pi60t)$, which leaves you with only
$ y(t) = \sin(2\pi10t) + \sin(2\pi50t)$.

Apply the same method one more time with $K = 100$ and you will end up with $y(t)=\sin(2\pi10t)$, although scaled and delayed due to the filters. The same method works for each of the four harmonics.

It might not be the most straightforward way to achieve this result, but it only uses simple filter structures. Another approach could be to use notch filters to remove the unwanted harmonics from your signal.


Can a narrowband beamformer be used for speech separation

If you are dealing with say 8 kHz for a nominal speech bandwidth of of 4kHz i.e. 0 Hz - 4 kHz, then the speech is essentially a wideband signal. Therefore narrowband beamforming won't work very well. Your beam pattern will be okay for the particular frequency of interest - but once you start moving away from that frequency your beam patterns will deteriorate.

What you are looking for is wideband beamforming - which uses actual time delays (or linear phase shifts across frequency in the frequency domain) rather than just phase shifts.

Usually the techniques like LCMV and MVDR are developed for narrowband signals. There are a couple of ways of approaching the wideband problem:

Use a tapped delay line on each channel. If you have $n$ channels and $m$ taps per channel then your correlation matrix is $mn$ x $mn$. So the systems of equations gets very large.
Use a bunch of narrowband beamformers. In this case you would have $m$ beamformers ($m$ frequencies) each one having $n$ channels. Now each beamformer has a $n$ x $n$ correlation matrix, but you have $m$ of them. So it leads to a reduction in complexity from the previous case.
You can form a set of $b$ conventional beams (using time-delay rather than just phase multiplier) and then do the adaptive processing on the beams.

The best reference I can think of is by Van Trees - Optimum Array Processing. Note you sometimes run into slight differences in terminology - some texts will denote MVDR as using the signal correlation matrix, while others use the signal plus noise correlation matrix - Just be careful which one you're looking at. I know Van Trees makes a distinction between the two cases. Most other texts do not and only use one formulation and call it MVDR.


Receiving multiple signals at the same frequency at the same time

One basic concept is that if multiple signals are mutually orthogonal in some sense (time domain, frequency domain, encoding, spreading function, phase relationship, etc.) they won't interfere in that domain, and thus can't cancel.


How to compare PSD methods?

Same way you deal with any type of engineering problem of this sort:

Write down requirements of your application. What are the priorities, sensitivities and minimum quantitative requirements/thresholds.?
Write down the properties, features and parameter space of the different PSD approaches that you have found.
Compare #1 and #2, find the PSD method and parameters that meets your requirements the best and/or has the best set of tradeoffs.


I'm sure that you are familiar with the idea of how Doppler affects a sine wave.  When the transmitter and/or receiver move towards the other the wave arrives more quickly, which makes the wave appear to be a higher frequency.  In regards to the baud rate, it also arrives more quickly.  If you think about it you will see that this must be so.  After all, the data stream cannot be separated from the carrier wave that modulates it!

So, you have a signal that is exactly identical to the "non-Doppler" signal except it arrives faster or slower than the non-Doppler signal.  This can be modeled exactly, in the discrete time domain, through resampling.

So why is Doppler usually modeled as a simple frequency shift when it should really be modeled as resampling?  Two reasons- first, the frequency shift effect is almost always vastly larger than the baud rate change.  The reason for this is that the Doppler effect is proportional to the frequency of the original wave ($\Delta f = \frac{\Delta v}{c}f_0$).  The frequency of the carrier wave is vastly greater than the baud rate, so the effect is much larger on the carrier wave.

The second reason most people don't model it as resampling is that it is difficult to do very small arbitrary sample rate changes, like 3013/3014.  If either the numerator or denominator is a prime number or has really large factors then it can become impractical to do traditional fractional resampling.  In that case polynomial interpolation is probably the best way to go.


Effective Nyquist frequency / signal resolution for overlapping time averaged signal samples

The process that you are describing can be broken down into two steps

Filter with a moving average filter with length $t_{width}$
Sampling with an interval $delta\_t$

A moving average filter is a low pass filter, but not a very good with lots of side lobes. It will help reduce aliasing in the sample process, but it won't eliminate it with any type of hard "threshold". The longer $t_{width}$ is and the the smaller $delta\_t$ is the less aliasing you get.

I think you approaching this backwards: Your "time resolution" is simply determined by the highest frequency in the signal. You need to determine this first and work backwards from there.

In practice, there is always aliasing. No real world signal is actually bandwidth limited in the mathematical sense. So you need to determine the frequency range "of interest", figure out how much energy the signal has outside this region, determine how much noise or aliasing you can tolerate, and design a suitable anti-aliasing filter to meet your specific requirements.


Nonlinear time-invariant frequency doubler

As already suggested by Robert and Olli, a system that maps $x(t)=k\cos(2\pi f_0t)$ to $y(t)=k\cos(4\pi f_0 t)$ can be formalized as

$$y(t)=|x(t)|_{max}\left(2\left(\frac{x(t)}{|x(t)|_{max}}\right)^2-1\right)\tag{1}$$

which is a time-invariant non-linear system.

However, I doubt that this system works well (i.e., sounds good as a distortion effect) when applied to real-world (non-sinusoidal) input signals.

To get some ideas of digital implementations of distortion effects, take a look at this thesis and the references therein.


i do not have a reference for you, but you are sampling among every one of the dimensions and the Nyquist Criterion must be satisfied for each one.

so fix $v$ to an arbitrary value $v_0$, then check if $\tilde{f}(u,v_0)$ is bandlimited some bandlimit, we'll call $B_u$.  find the max value of $B_u$ for all possible $v_0$, double that and that is the minimum sample rate along the $u$-axis.

then swap the roles of $u$ and $v$, fix $u$ to $u_0$, find the bandlimit of $\tilde{f}(u_0,v)$ for every possible value of $u_0$.  find the maximum (over all $u_0$) bandlimit $B_v$, double that and that is the minimum-sample rate along the $v$-axis.

if your curve can move at an arbitrary angle w.r.t. the $u$-axis and $v$-axis, then pick the maximum of $B_u$ and $B_v$ and sample both axes at double of that rate (or slightly higher).  then, along the path length, you must sample at a minimum of that rate, per unit of arc length, to satisfy Nyquist in two dimensions.


Does each pixel of the original image has a corresponding coefficient for every scale and every direction calculated by the DWT?

No, by essence, as you guessed correctly. There is a direct counting argument: given a $2N\times 2N$ image, the first level of the  DWT generates  three $N\times N$ wavelet subbands. And through the pigeonhole principle, you cannot have a bijection between a $4N^2$ and a $3N^2$ coefficient set.

As for the contourlet, you can resort to undecimated wavelet versions, for which several implementations are possible. Some allow more directions than those obtained through traditional tensor product.

However, I consider that it is misleading to believe that one isolated pixel could be easily followed across scales or directions and mapped to one coefficient: this requires smoothing or directional filters taking into account weighted combinations of the surrounding pixels. In other words, one should not expect arbitrary precision on location AND arbitrary scale or orientation determination.


Yes you are correct, differential modulations such as DBPSK, DQPSK, $pi/4$ DQPSK do not need to perform carrier recovery as you can compare one symbol to the next to perform demodulation using a delay and multiply approach. This results in a simple and potentially low power receiver, at the expense of performance compared to coherent modulation schemes.

As far as a static carrier offset this should be no issue- the worst case is if the carrier frequency is changing from one symbol to the next (and at some level it always is, in particular phase noise)- in this regard simply consider how much your symbol rotates from one symbol to the next when you demodulate with this approach and you can directly see the impact of frequency offset (and if the SNR degradation from that rotation based on the resulting symbol distance would be acceptable). For example, a 2 degree rotation would be negligible, while with a 45 degree rotation you be on the threshold of an incorrect symbol decision with no other added noise (using the demodulation approach I describe below). The impact of this can be quantified from the phase noise, and how I approach that with a differential demodulation is to consider with regards to phase how a multiplier for small angles is a "phase adder" thus the delay and multiply demod process is a high pass filter with regards to phase noise. (The filter is given by $1-z^{-1}$). It has a gain of 2 at the high frequency portion (symbol rate) and a null at DC- so it passes (and doubles) the phase noise at those higher offsets). Ultimately you can integrate the high pass filtered phase noise from DC out to the symbol rate and double that to determine the rms phase error, and determine performance issues with consideration to the impact of phase offsets as I described above. What you will find in doing this is the demodulator will likely be very robust against phase noise.

An approach I would consider is doing a complex conjugate multiplication of each new symbol with the previous symbol with the additional +/-45 degree rotation in that process to remove the pi/4 offsets. The output of this would be I and Q with bi-phase signaling on each (similar to QPSK) that you could easily pass into a Gardner (at 2 samples per symbol) for timing recovery followed by matched filters and then data decision. An alternate approach used for low cost / low power non-coherent pi/4 QSPK demodulation consists of hard limiting the waveform and then using a simple FM discriminator followed by a four level detector (since the 4 phase transitions +45°, -45°, +135°, -135°  from each symbol to the next correspond to 4 possible frequencies).


The Standard Deviation of The Derivative of a Signal

As pointed out in pichenettes' answer, you need to know more about the original signal $x(t)$. If the signal's power spectrum $S_x(\omega)$ is known then the power spectrum of its $n^{th}$ derivative is

$$S_y(\omega)=\omega^{2n}S_x(\omega)$$

The power of the $n^{th}$ derivative can then be computed as

$$E\{|y(t)|^2\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\omega^{2n}S_x(\omega)d\omega\tag{1}$$

(assuming convergence of the integral, e.g. by assuming that $x(t)$ is band-limited). Since the mean of the derivative is zero (even if the mean of the original signal isn't), the variance of $y(t)$ equals its power, and so the standard deviation is the square root of the expression in (1).