%pylab inline

Populating the interactive namespace from numpy and matplotlib

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)

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

fidx = array([979,983,987,991])-1 # 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]
ivec = array(ivec)-1
jvec = array(jvec)-1

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 = array([1,3,5,7])-1
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

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+1,N))
y = zeros(2*Nb+1)
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])
H[-1,:] = 1.
y[-1] = 1.

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

<matplotlib.text.Text at 0x7faa1f13ae90>

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

U,s,VT = np.linalg.svd(H)
Vs = VT.T
r = len(s)
Ns = Vs[r:,:]

Direct Fourier Inversion

fstar = H.T.dot(y)
fstar = fstar/sum(fstar)
plot(180/pi*th,fstar,'k')
xlabel('Angle [deg]')
ylabel('Normalized Power')
title('Direct Fourier "Solution"')

<matplotlib.text.Text at 0x7faa1f1f7c50>

Min 2-norm

$$~~~~\textrm{min} ~~~~~~~~~||c||_2 $$$$ \textrm{such that}~~~\textbf{g} = M \textbf{c} $$

lam0 = 1e-10
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('Min 2-norm')

<matplotlib.text.Text at 0x7faa1f031ad0>

Min 1-norm

$$~~~~\textrm{min} ~~~~~~~~~||c||_1 $$$$ \textrm{such that}~~~\textbf{g} = M \textbf{c} $$

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

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

<matplotlib.text.Text at 0x7faa1ef68490>

def get_slice(tidx, hidx, N=200, fidx = array([979,983,987,991])-1):
    ''' 
    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]
    ivec = array(ivec)-1
    jvec = array(jvec)-1
    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 = array([1,3,5,7])-1
    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])
        

    fstar = H.T.dot(y)
    

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

    return f,th

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
    

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

TH,H = meshgrid(th,hts)
X = H*TH
Y = H*sqrt(1-TH**2)

figure(figsize=(16,10))
Ip = 10*log10(Imask)
pcolormesh(X,Y,Ip,cmap=cmap)
clim((-10,Ip.max()))
cb = colorbar()
cb.set_label('Power [dB]')
axis('equal')
xlabel('Horizontal Distance [km]')
ylabel('Vertical Distance [km]')
#axis('tight')
title('Direct Fourier "Solution"')

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

<matplotlib.text.Text at 0x7faa1ede4d50>

def get_slice(tidx, hidx, N=200, fidx = array([979,983,987,991])-1):
    ''' 
    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]
    ivec = array(ivec)-1
    jvec = array(jvec)-1
    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 = array([1,3,5,7])-1
    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-5 # 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

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
    

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

TH,H = meshgrid(th,hts)
X = H*TH
Y = H*sqrt(1-TH**2)

figure(figsize=(16,10))
Ip = 10*log10(Imask)
pcolormesh(X,Y,Ip,cmap=cmap)
clim((-10,Ip.max()))
cb = colorbar()
cb.set_label('Power [dB]')
axis('equal')
xlabel('Horizontal Distance [km]')
ylabel('Vertical Distance [km]')
#axis('tight')
title('Min 1-norm solution')

/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

<matplotlib.text.Text at 0x7faa1ec45750>

import subprocess
import shutil
nb_name = 'ECE458_Lecture_graphics.ipynb' # There's probably a way to get this programmatically
cmd = ['ipython', 'nbconvert', '--to', 'slides', nb_name, \
       '--reveal-prefix', '"http://cdn.jsdelivr.net/reveal.js/2.6.2"',
      '--template','output_toggle']
# The --template command points to a custom template which opens and closes the input code cells with a click
subprocess.call(cmd)
html_name = '%s.slides.html' % (nb_name.split('.')[0])
shutil.copy(html_name, '/home/bhardin2/public_html/slideshows/')