phase-amplitude-extraction-test.py

import numpy as np
import scipy
from scipy import signal
import pylab as pl

if __name__ == "__main__":
    N  = 1024*8
    A  = 2**14
    An = 400
    phase = np.pi/5
    fs = 250E6
    t = [(i/fs) for i in range(N)]

    window = signal.flattop(N) #np.hamming(N)
#    f_current = 300E6
    for f_current in np.arange(50E6,500E6,10E6):
        if f_current % 125E6 == 0:
            continue
        print("")
        print("%g Hz"%(f_current))
        noise1 = np.random.normal(0,1,N)
        noise2 = np.random.normal(0,1,N)
        
        ref1 = [A*np.sin(2*np.pi*f_current*i) for i in t]
        ref1 = [(An*noise1[i] + ref1[i])*window[i] for i in range(N)]
        comp = [.8*A*np.sin(2*np.pi*f_current*i + phase) for i in t]
        comp = [(An*noise2[i] + comp[i])*window[i] for i in range(N)]


        ref1_f = np.fft.rfft(ref1)
        comp_f = np.fft.rfft(comp)
      #  pl.plot(ref1)
        #pl.show()
        #pl.subplot(1,2,1)
        #pl.plot(abs(ref1_f))
        #pl.plot(abs(comp_f))
        #pl.subplot(1,2,2)
        #pl.plot(np.angle(ref1_f))
        #pl.plot(np.angle(comp_f))
        
        f_div = fs/(N)
        
        nyquist_region = int(f_current/(.5*fs))
        f_mod = f_current % fs
        if (nyquist_region % 2) == 1:
            idx = N/2 - int((f_mod/f_div))
        else:
            idx = int((f_mod/f_div))
            #print idx
            
        print(abs(comp_f[idx])/abs(ref1_f[idx]))
        print((np.angle(comp_f[idx])-np.angle(ref1_f[idx])))

#        pl.show()