Radar HF

#!/usr/bin/env python

import numpy as n

import matplotlib.pyplot as plt

import stuffr

import scipy.signal as s

import glob

#Input gdfdata

fl=glob.glob("/home/jm/Documents/PruebasSincronizacion/datarxcode02A/d2018289_1029/0/data*.gdf")

#fl=glob.glob("/home/jm/Documents/PruebasSincronizacion/datarxcode2B/d2018289_1034/0/data*.gdf")

#fl=glob.glob("/home/jm/Documents/PruebasSincronizacion/datarxcode2oroya/d2018289_1109/0/data*.gdf")

fl.sort()

fl=fl[-11:-1] # Just test with the last 10 files 

z=n.zeros(100000*11,dtype=n.complex64)

for i in range(11):

    print 'i:',i

    print i*100000+n.arange(100000) # rotating index from vector

    z[i*100000+n.arange(100000)]=n.fromfile(fl[i-1],dtype=n.complex64)

clen=10000 #Bauds per frame.

n_range=1000 #Number of heights

# Generating PR Code

code2=stuffr.create_pseudo_random_code(seed=2)

code2=n.roll(code2,8) #Rollingback 8 laggs

print 'len:code2:', len(code2)

zl=len(z) #Datavector lenght

##Previus check to know if the usrps aren't synchronized

#for i in range(10000):

#    print("i %d %1.2f"%(i,n.abs(n.sum(n.conj(code2)*z[n.arange(10000)+i]))))

dt=10000

#Correlation with the periodic Correlation matrix.

A=stuffr.periodic_convolution_matrix(code2,rmin=0,rmax=n_range)["A"]

# maximum likelihood: Juha's algorithm that estimate three sets of ranges

B=n.dot(n.linalg.inv(n.dot(n.transpose(n.conj(A)),A)),n.conj(n.transpose(A)))

#print(B.shape)

n_rep=100 # number of frames or traces per second.

idx=n.arange(clen)

ridx=n.arange(n_range)

S2=n.zeros([n_rep,n_range],dtype=n.complex64)

for i in range(n_rep):

    # there is a range delay. usrps are not locked.

    xhat=n.dot(B,z[i*clen+idx+dt])

    x0=xhat[ridx]

    S2[i,:]=x0    

SP2=n.zeros([n_rep,n_range])

w=s.hann(n_rep)

for i in range(n_range):

    SP2[:,i]=n.fft.fftshift(n.abs(n.fft.fft(w*S2[:,i]))**2.0)

# return SNR

def whiten_spec(S):

    n_rep=S.shape[0]

    # whiten spectrum

    for i in range(n_rep):

        m0=n.median(S[i,:])

        s0=n.median(n.abs(S[i,:]-m0))

        S[i,:]=(S[i,:]-m0)/s0

    return(S)

#Just plotting Code2 decodification

plt.pcolormesh(10.0*n.log10(SP2).transpose(),cmap='jet')

plt.colorbar()

plt.show()

#whitening spectrum =)

SP2=whiten_spec(SP2)

plt.pcolormesh(10.0*n.log10(SP2).transpose(),vmin=0,vmax=100,cmap='jet')

plt.colorbar()

plt.show()