Problem 1¶

In [3]:

fn = '/home/bhardin2/ece458/VIPIR_10min_UT_2015014143004.npz.npz'
npzfile = load(fn)
print npzfile.keys()
npzfile.close()

['UT_times', 'freqs_kHz', 'volts', 'DST', 'time_zone', 'channels', 'hts_km']

For the JRO instrument we have the following variables:

In [4]:

fn = '/home/bhardin2/ece458/JRO_10min_Valley2015014_11.npz.npz'
npzfile = load(fn)
print npzfile.keys()
npzfile.close()

['epochtime', 'DST', 'time_zone', 'heights', 'data_ch0', 'data_ch1']

Problem 2¶

In [6]:

fn = '/home/bhardin2/ece458/VIPIR_10min_UT_2015014143004.npz.npz'
npzfile = load(fn)
V = npzfile['volts']
print shape(V) # (runs, pulses, ranges, antennas)

(10, 1632, 512, 8)

In [8]:

freqs = npzfile['freqs_kHz']
hts   = npzfile['hts_km']

In [33]:

fg = list(set(freqs))
fg.sort()
hg = hts.copy()
F,H = meshgrid(fg,hg)
P = zeros(shape(F))
Nh,Nf = shape(I)
V0 = V[0,:,:,0] # 1st time, 1st antenna

In [34]:

for j in range(Nf):
    idx = freqs == fg[j]
    for i in range(Nh):
        V0i = V0[idx, i]
        P[i,j] = mean(abs(V0i)**2)

In [45]:

pcolormesh(F/1e3, H, P, norm=mpl.colors.LogNorm())
colorbar()#shrink=0.7, aspect=5)
axis('tight');
xlabel('Freq [MHz]')
ylabel('Height [km]')

Out[45]:

<matplotlib.text.Text at 0x7f5d9de11450>

Problem 3¶

Part a¶

In [63]:

fidx = array([979,983,987,991])
h = 292
hidx = argmin(abs(hts-h))

ivec = [1, 3, 1, 7, 7, 7]
jvec = [3, 5, 5, 1, 3, 5]
visib_vec = zeros(len(ivec))
N = len(ivec)
Nt = 10 # num times

visibs = zeros((N,Nt),dtype=complex)
P = zeros((N,Nt))

for ii in range(N):
    for k in range(Nt):

        i = ivec[ii] # ant 1
        j = jvec[ii] # ant 2

        Vi = V[k,fidx,hidx,i]
        Vj = V[k,fidx,hidx,j]

        P_ij = sqrt( mean(abs(Vi)**2) * mean(abs(Vj)**2))
        visib_ij = mean(Vi.conjugate() * Vj) / P_ij
        
        visibs[ii,k] = visib_ij
        P[ii,k] = P_ij

In [72]:

for ii in range(N):
    C = abs(visibs[ii,:]) # coherence
    plot(C,label='%i,%i'%(ivec[ii],jvec[ii]))
legend(loc='best')
xlabel('Time [min]')
ylabel('Coherence')
ylim((0.0,1.1))

Out[72]:

(0.0, 1.1)

In [80]:

for ii in range(N):
    ph = 180/pi*angle(visibs[ii,:]) # phase
    plot(ph,label='%i,%i'%(ivec[ii],jvec[ii]))
legend(loc='best')
xlabel('Time [min]')
ylabel('Phase [deg]')

Out[80]:

<matplotlib.text.Text at 0x7f5d93f6e210>

In [75]:

for ii in range(N):
    plot(P[ii],label='%i,%i'%(ivec[ii],jvec[ii]))
legend(loc='best')
xlabel('Time [min]')
ylabel('Power')

Out[75]:

<matplotlib.text.Text at 0x7f5d93fe7a50>

Part b¶

The coherence is very close to 1 when the signal power is strong. Since this is true for all antenna separations, this suggests that we are dealing with a single plane wave source.

Part c¶

We'll use the longest baseline (5,7) to determine the direction of arrival, ignoring phase calibration errors. We'll calculate it for the 5th minute, when there is large signal power.

In [82]:

freqs[fidx]

Out[82]:

array([ 5485.50195312,  5485.50195312,  5485.50195312,  5485.50195312], dtype=float32)

In [84]:

ii = 5
kk = 5
phase = angle(visibs[ii,kk])
lam = 3e8/(1e3*freqs[fidx][0])
k = 2*pi/lam
dx = 60+45
theta = 180/pi * arcsin(phase/(k*dx))
print 'Angle of arrival is %.1f degrees' % theta

Angle of arrival is -1.3 degrees

Part d: Try Imaging¶

In [525]:

fidx = array([979,983,987,991]) # frequency
hidx = 190 # height
tidx = 5 # time
ivec = [1, 3, 1, 7, 7, 7]
jvec = [3, 5, 5, 1, 3, 5]
dxvec = [15, 30, 45, 60, 75, 105]

f = 1e3*freqs[fidx][0] # Hz
lam = 2.9979e8/f
    
visib_vec = zeros(len(ivec))
Nb = len(ivec)

visibs = zeros(Nb,dtype=complex)
Pa = zeros(4)
avec = [1,3,5,7]
for ai in range(len(avec)):
    Vai = V[tidx,fidx,hidx,avec[ai]]
    Pa[ai] = mean(abs(Vai)**2)

for ii in range(Nb):

    i = ivec[ii] # ant 1
    j = jvec[ii] # ant 2

    Vi = V[tidx,fidx,hidx,i]
    Vj = V[tidx,fidx,hidx,j]

    P_ij = sqrt( mean(abs(Vi)**2) * mean(abs(Vj)**2))
    visib_ij = mean(Vi.conjugate() * Vj) / P_ij

    visibs[ii] = visib_ij

In [526]:

N = 200 # number of pixels
th_max = 45*pi/180
th = linspace(-th_max, th_max, N)
# Construct matrix and RHS
H = zeros((2*Nb,N))
y = zeros(2*Nb)
for ii in range(Nb):
    phase = 2*pi*th*dxvec[ii]/lam
    H[2*ii,:] = cos(phase)
    y[2*ii] = real(visibs[ii])
    H[2*ii+1,:] = sin(phase)
    y[2*ii+1] = imag(visibs[ii])

In [527]:

TH,NB = meshgrid(180/pi*th,range(2*Nb))
pcolormesh(TH,NB,H,cmap='RdBu_r')
colorbar()
axis('tight')
xlabel('Degrees')
ylabel('Measurement index')
title('Measurement matrix')

Out[527]:

<matplotlib.text.Text at 0x7f5d57a84690>

So our (severely underdetermined) matrix equation is:

$$ \textbf{y} = \textbf{H} \textbf{f} $$

How to solve it? First, try minimum 2-norm solution.

Minimum 2-norm solution¶

In [528]:

lam0 = 1e3
L = eye(N) # penalize 2-norm of vector
M = H.T.dot(H) + lam0*L.T.dot(L)
fstar = np.linalg.solve(M, H.T.dot(y))

plot(180/pi*th,fstar,'k')
xlabel('Angle [deg]')
ylabel('Normalized Power')
title('Appx min 2-norm solution')

Out[528]:

<matplotlib.text.Text at 0x7f5d57a02e90>

Minimum 1-norm solution¶

Next try minimum 1-norm solution (to promote sparsity). Although it's possible to solve this better by incorporating knowledge of the noise, for now I'll just solve it approximately using a Basis-Pursuit-style inversion similar to CLEAN.

In [529]:

from sklearn import linear_model
alpha = 1e-10 # make it small so as to fit data ~exactly
clf = linear_model.Lasso(alpha=alpha,fit_intercept=False, positive=True)
clf.fit(H,y)
fstar = clf.coef_

plot(180/pi*th,fstar,'k')
xlabel('Angle [deg]')
ylabel('Normalized Power')
title('Appx min 1-norm solution')

Out[529]:

<matplotlib.text.Text at 0x7f5d579afad0>

MaxEnt solution¶

Another option is to use the Maximum Entropy method. This chooses the image with the maximum entropy (least information) that fits the measured data. This method is most useful when we have a good characterization of the errors, which I didn't have time to think about for this homework. Thus, I set the errors to be very small, which gives us rather low entropy images (such as spikes). Also, this method does not always seem to converge, which might be because of the large errors, or because of the imperfect phase calibration.

In [530]:

def modelf(lam, H, F):
    exponent = H.T.dot(lam)    
    E = max(exponent)
    fnum = exp(exponent-E)
    #print max(exponent)
    f = F*fnum/sum(fnum)
    return f

In [531]:

M = 2*Nb
F = 1. # values must sum to 1


def myfunc(lam_vals):

    sigma = 0.000001*ones(M)
    G = M

    f = modelf(lam_vals, H, F)

    # solve for Lambda and e
    Lambda = sqrt(sum(lam_vals**2 * sigma**2)/(4*G))
    e = -lam_vals*sigma**2/(2*Lambda)
    
    return y + e - H.dot(f)

In [532]:

M = 2*Nb
F = 1. # values must sum to 1


def myfunc(lam_vals):

    sigma = 0.000001*ones(M)
    G = M

    f = modelf(lam_vals, H, F)

    # solve for Lambda and e
    #Lambda = sqrt(sum(lam_vals**2 * sigma**2)/(4*G))
    #e = -lam_vals*sigma**2/(2*Lambda)
    
    return y - H.dot(f)

In [533]:

from scipy.optimize import root
lam_init = rand(M)*1e-5*ones(M)
result = root(myfunc, lam_init, method='lm', options={'maxiter':10000})
print result.message
lamstar = result.x

Both actual and predicted relative reductions in the sum of squares
  are at most 0.000000

In [534]:

fstar = modelf(lamstar, H, F)
plot(180/pi*th,fstar,'k')
xlabel('Angle [deg]')
ylabel('Normalized Power')
title('MaxEnt solution')

Out[534]:

<matplotlib.text.Text at 0x7f5d578d4e90>

Loop over all heights to create image (using L1-regularized inversion)¶

In [487]:

def get_slice(tidx, hidx, N=200, fidx = array([979,983,987,991])):
    ''' 
    Return a 1D vector representing the
    angular distribution of scattered power
    at a given time, height, and frequency
    '''

    #fidx = array([979,983,987,991]) # frequency
    ivec = [1, 3, 1, 7, 7, 7]
    jvec = [3, 5, 5, 1, 3, 5]
    dxvec = [15, 30, 45, 60, 75, 105]

    f = 1e3*freqs[fidx][0] # Hz
    lam = 2.9979e8/f

    visib_vec = zeros(len(ivec))
    Nb = len(ivec)

    visibs = zeros(Nb,dtype=complex)
    Pa = zeros(4)
    avec = [1,3,5,7]
    for ai in range(len(avec)):
        Vai = V[tidx,fidx,hidx,avec[ai]]
        Pa[ai] = mean(abs(Vai)**2)

    for ii in range(Nb):

        i = ivec[ii] # ant 1
        j = jvec[ii] # ant 2

        Vi = V[tidx,fidx,hidx,i]
        Vj = V[tidx,fidx,hidx,j]

        P_ij = sqrt( mean(abs(Vi)**2) * mean(abs(Vj)**2))
        visib_ij = mean(Vi.conjugate() * Vj) / P_ij

        visibs[ii] = visib_ij

    # Formulate problem

    th_max = 45*pi/180
    th = linspace(-th_max, th_max, N)
    # Construct matrix and RHS
    H = zeros((2*Nb,N))
    y = zeros(2*Nb)
    for ii in range(Nb):
        phase = 2*pi*th*dxvec[ii]/lam
        H[2*ii,:] = cos(phase)
        y[2*ii] = real(visibs[ii])
        H[2*ii+1,:] = sin(phase)
        y[2*ii+1] = imag(visibs[ii])

    # Solve problem
    alpha = 1e-10 # make it small so as to fit data ~exactly
    clf = linear_model.Lasso(alpha=alpha,fit_intercept=False, positive=True)
    clf.fit(H,y)
    fstar = clf.coef_

    # Un-normalize by total power
    if not all(fstar==0):
        P = mean(Pa)
        f = P*fstar/sum(fstar)

    return f,th

In [514]:

tidx = 5 # time
N = 200 # Ntheta
Nh = len(hts)

I = zeros((Nh,N))

for hidx in range(Nh):
    f,th = get_slice(tidx,hidx,N=N)
    I[hidx,:] = f

In [515]:

Imask = ma.masked_array(I, mask = I==0)
cmap = cm.hot
cmap.set_bad('k',1)

TH,H = meshgrid(th,hts)
X = H*sin(TH)
Y = H*cos(TH)

figure(figsize=(8,5))
Ip = 10*log10(Imask)
pcolormesh(X,Y,Ip,cmap=cmap)
clim((-20,Ip.max()))
cb = colorbar()
cb.set_label('Power [dB]')
axis('equal')
xlabel('Horizontal Distance [km]')
ylabel('Vertical Distance [km]')

/usr/local/lib/python2.7/dist-packages/ipython-3.0.0_b1-py2.7.egg/IPython/kernel/__main__.py:11: RuntimeWarning: divide by zero encountered in log10

Out[515]:

<matplotlib.text.Text at 0x7f5d788ab8d0>

Loop over time¶

In [518]:

figure(figsize=(6*3,4*3))
for tidx in range(9):
    
    subplot(3,3,tidx)

    N = 200 # Ntheta
    Nh = len(hts)

    I = zeros((Nh,N))

    for hidx in range(Nh):
        f,th = get_slice(tidx,hidx,N=N)
        I[hidx,:] = f


    Imask = ma.masked_array(I, mask = I==0)
    cmap = cm.hot
    cmap.set_bad('k',1)

    TH,H = meshgrid(th,hts)
    X = H*sin(TH)
    Y = H*cos(TH)

    Ip = 10*log10(Imask)
    pcolormesh(X,Y,Ip,cmap=cmap)
    clim((-20,Ip.max()))
    cb = colorbar()
    cb.set_label('Power [dB]')
    axis('equal')
    xlabel('Horizontal Distance [km]')
    ylabel('Vertical Distance [km]')

/usr/local/lib/python2.7/dist-packages/ipython-3.0.0_b1-py2.7.egg/IPython/kernel/__main__.py:24: RuntimeWarning: divide by zero encountered in log10

Zoom in to see reflection point(s)¶

In [522]:

figure(figsize=(6*3,4*3))
for tidx in range(9):
    
    subplot(3,3,tidx)

    N = 200 # Ntheta
    Nh = len(hts)

    I = zeros((Nh,N))

    for hidx in range(Nh):
        f,th = get_slice(tidx,hidx,N=N)
        I[hidx,:] = f


    Imask = ma.masked_array(I, mask = I==0)
    cmap = cm.hot
    cmap.set_bad('k',1)

    TH,H = meshgrid(th,hts)
    X = H*sin(TH)
    Y = H*cos(TH)

    Ip = 10*log10(Imask)
    pcolormesh(X,Y,Ip,cmap=cmap)
    clim((-20,Ip.max()))
    cb = colorbar()
    cb.set_label('Power [dB]')
    xlim((-75,75))
    ylim((225,375))
    #axis('equal')
    xlabel('Horizontal Distance [km]')
    ylabel('Vertical Distance [km]')

/usr/local/lib/python2.7/dist-packages/ipython-3.0.0_b1-py2.7.egg/IPython/kernel/__main__.py:24: RuntimeWarning: divide by zero encountered in log10

In [ ]: