<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Spectrum Estimation
by John Swoboda
This notebook will use portions of the SimISR and ISRSpectrum software modules to create IS spectra and show different ways of estimating them. There will also be examples using CWGN from the random number generator in scipy.
Import Everything
Step1:  Set up
Set the number of points desired for the averaging and length of spectra. Also set up seaborn formats.
Step2:  IS Spectra
This will create an example ISR spectra that will be used.
Step3:  White Noise
A periodogram is applied to complex white Gaussian Noise. This is here in order to show that the output of the scipy random number generator outputs uncorrelated random variables.
Step4:  Shaped Noise
A set of shaped noise is created using the IS spectrum formed earlier. Using linear pridictive coding to apply the spectrum to the noise through the MakePulseDataRepLPC function. This is similar to the method used by vocoders to encode human speech. To show the effect of the LPC coloring a periodogram esitmator is applied to noise the noise.
Step5:  Window Function
When a long pulse is used in ISR the ACF is esitmated instead of the spectrum directly through the periodogram estimator. The estimation is a two step process, first estimate the lags and then do a summation rule. This leads to a windowing of the ACF shown in this. The window is also shown in the frequency domain which will be applied as a convolution to the original spectra in frequency.
Step6:  Full ISR Data Creation and Estimator
The basics data creation and processing behind SimISR for this case for only one beam. The data is created along a set of samples by adding together a set of uncorrelated data sets together. These sets of pulses that are added together are uncorrelated because any spatial correlation of the electron density fluctuations are much smaller than a range gate. After that the ACFs are estimated they are plotted with the input ACF and spectra with the window applied.
Step7:  Plotting and Normalization of Input Spectra
Need to apply window function to spectrum.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import scipy as sp
import scipy.fftpack as scfft
from SimISR.utilFunctions import makesumrule,MakePulseDataRepLPC,spect2acf,acf2spect,CenteredLagProduct
from SimISR.IonoContainer import IonoContainer,MakeTestIonoclass
from ISRSpectrum.ISRSpectrum import ISRSpectrum
import seaborn as sns
Explanation: Spectrum Estimation
by John Swoboda
This notebook will use portions of the SimISR and ISRSpectrum software modules to create IS spectra and show different ways of estimating them. There will also be examples using CWGN from the random number generator in scipy.
Import Everything
End of explanation
# Processings parameters
spfreq=50e3 # Bandwidth
nspec=256 # length of spectrum
rep1=10000 # number of pulses
L=24. # Length of pulse in standard processings
pulse=sp.ones(int(L)) # Pulse for standard processing
pulse_pergram=sp.ones(nspec) # For periodogram
Nrg=128 # Number of Range gates for data
# Parameters for spectrum
species=['O+','e-']
databloc = sp.array([[1.66e10,1e3],[1.66e10,2.5e3]])
f_c = 440e6
# setup seaborne
sns.set_style("whitegrid")
sns.set_context("notebook")
Explanation: Set up
Set the number of points desired for the averaging and length of spectra. Also set up seaborn formats.
End of explanation
# Make spectrum
ISpec_ion = ISRSpectrum(centerFrequency=f_c, nspec=nspec, sampfreq=spfreq, dFlag=False)
f,cur_spec,rcs = ISpec_ion.getspecsep(databloc,species,rcsflag=True)
specsum = sp.absolute(cur_spec).sum()
cur_spec = len(cur_spec)*cur_spec*rcs/specsum
tau,acf=spect2acf(f,cur_spec)
fig,ax = plt.subplots(1,2,sharey=False, figsize=(8,4),facecolor='w')
rp,imp=ax[0].plot(tau*1e3,acf.real,tau*1e3,acf.imag)
ax[0].legend([rp,imp],['Real','Imag'])
ax[0].set_ylabel('Amplitude')
ax[0].set_title('ACF')
ax[0].set_xlabel(r'$\tau$ in ms')
ax[1].plot(f*1e-3,cur_spec.real)
ax[1].set_ylabel('Amplitude')
ax[1].set_title('Spectrum')
ax[1].set_xlabel(r'f in kHz')
fig.tight_layout()
Explanation: IS Spectra
This will create an example ISR spectra that will be used.
End of explanation
xin =sp.power(2,-.5)*(sp.random.randn(rep1,nspec)+1j*sp.random.randn(rep1,nspec))
Xfft=sp.power(nspec,-.5)*scfft.fftshift(scfft.fft(xin,axis=-1),axes=-1)
Xperiod=sp.power(Xfft.real,2).mean(0) +sp.power(Xfft.imag,2).mean(0)
tau2,acfperiod=spect2acf(f,Xperiod*nspec)
fig2,ax2 = plt.subplots(1,2,sharey=False, figsize=(8,4),facecolor='w')
rp,imp=ax2[0].plot(tau2*1e6,acfperiod.real,tau2*1e6,acfperiod.imag)
ax2[0].legend([rp,imp],['Real','Imag'])
ax2[0].set_ylabel('Amplitude')
ax2[0].set_title('ACF')
ax2[0].set_xlabel(r'$\tau$ in $\mu$s')
ax2[1].plot(f*1e-3,Xperiod.real)
ax2[1].set_ylabel('Amplitude')
ax2[1].set_title('Spectrum')
ax2[1].set_xlabel(r'f in kHz')
ax2[1].set_ylim([0.,1.5])
fig2.tight_layout()
Explanation: White Noise
A periodogram is applied to complex white Gaussian Noise. This is here in order to show that the output of the scipy random number generator outputs uncorrelated random variables.
End of explanation
Xdata = MakePulseDataRepLPC(pulse_pergram,cur_spec,30,rep1,numtype = sp.complex128)
Xfftd=sp.power(nspec,-.5)*scfft.fftshift(scfft.fft(Xdata,axis=-1),axes=-1)
Xperiodd=sp.power(Xfftd.real,2).mean(0) +sp.power(Xfftd.imag,2).mean(0)
tau3,acfperiodd=spect2acf(f,Xperiodd*nspec)
fig3,ax3 = plt.subplots(1,2,sharey=False, figsize=(8,4),facecolor='w')
rp,imp=ax3[0].plot(tau3*1e6,acfperiodd.real,tau3*1e6,acfperiodd.imag)
ax3[0].legend([rp,imp],['Real','Imag'])
ax3[0].set_ylabel('Amplitude')
ax3[0].set_title('ACF')
ax3[0].set_xlabel(r'$\tau$ in $\mu$s')
ax3[1].plot(f*1e-3,Xperiodd.real)
ax3[1].set_ylabel('Amplitude')
ax3[1].set_title('Spectrum')
ax3[1].set_xlabel(r'f in kHz')
fig2.tight_layout()
Explanation: Shaped Noise
A set of shaped noise is created using the IS spectrum formed earlier. Using linear pridictive coding to apply the spectrum to the noise through the MakePulseDataRepLPC function. This is similar to the method used by vocoders to encode human speech. To show the effect of the LPC coloring a periodogram esitmator is applied to noise the noise.
End of explanation
v=1
l=sp.arange(L)
W=-l**2/(L*v) + (L-v)*l/L/v+1
Wp=sp.pad(W,(int(sp.ceil(float(nspec-L)/2)),int(sp.floor(float(nspec-L)/2))),'constant',constant_values=0)
wfft=scfft.fftshift(scfft.fft(W,n=nspec))
fig4,ax4 = plt.subplots(1,2,sharey=False, figsize=(8,4),facecolor='w')
ax4[0].plot(l,W)  
ax4[0].set_ylabel('Weighting')
ax4[0].set_title('Weighting')
ax4[0].set_xlabel(r'$l$')
rp,imp,abp=ax4[1].plot(f*1e-3,wfft.real,f*1e-3,wfft.imag,f*1e-3,sp.absolute(wfft))
ax4[1].legend([rp,imp,abp],['Real','Imag','Abs'])
ax4[1].set_ylabel('Amplitude')
ax4[1].set_title('Spectrum')
ax4[1].set_xlabel(r'f in kHz')
fig4.tight_layout()
Explanation: Window Function
When a long pulse is used in ISR the ACF is esitmated instead of the spectrum directly through the periodogram estimator. The estimation is a two step process, first estimate the lags and then do a summation rule. This leads to a windowing of the ACF shown in this. The window is also shown in the frequency domain which will be applied as a convolution to the original spectra in frequency.
End of explanation
Xdata=sp.zeros((rep1,Nrg),dtype=sp.complex128)
Lint = int(L)
for i in range(int(Nrg-Lint)):
    Xdata[:,i:i+Lint] = MakePulseDataRepLPC(pulse,cur_spec,40,rep1,numtype = sp.complex128)+Xdata[:,i:i+Lint]
lagsData=CenteredLagProduct(Xdata,numtype=sp.complex128,pulse =pulse,lagtype='centered')/rep1
ptype='long'
ts = 1.
sumrule=makesumrule(ptype,L,ts,lagtype='centered')
minrg = -sumrule[0].min()
maxrg = Nrg-sumrule[1].max()
Nrng2 = maxrg-minrg
lagsDatasum = sp.zeros((Nrng2,Lint),dtype=sp.complex128)
for irngnew,irng in enumerate(sp.arange(minrg,maxrg)):
    for ilag in range(Lint):
        lagsDatasum[irngnew,ilag] = lagsData[irng+sumrule[0,ilag]:irng+sumrule[1,ilag]+1,ilag].sum(axis=0)
# divide off the gain from the pulse stacking        
lagsDatasum = lagsDatasum/L                       
Explanation: Full ISR Data Creation and Estimator
The basics data creation and processing behind SimISR for this case for only one beam. The data is created along a set of samples by adding together a set of uncorrelated data sets together. These sets of pulses that are added together are uncorrelated because any spatial correlation of the electron density fluctuations are much smaller than a range gate. After that the ACFs are estimated they are plotted with the input ACF and spectra with the window applied.
End of explanation
dt=tau[1]-tau[0]
f1,spec_all=acf2spect(l*dt,lagsDatasum,n_s=nspec)
acf_single = lagsDatasum[50]
spec_single = spec_all[50]
# Apply weighting and integrations from gain from pulse stacking 
acf_act=scfft.ifftshift(acf)[:Lint]*W
feh,spec_act=acf2spect(l*dt,acf_act,n_s=nspec)
fig5,ax5 = plt.subplots(1,2,sharey=False, figsize=(8,4),facecolor='w')
rp,imp,act_acf=ax5[0].plot(l*dt*1e6,acf_single.real,l*dt*1e6,acf_single.imag,l*dt*1e6,acf_act.real)
ax5[0].legend([rp,imp,act_acf],['Real','Imag','Actual'])
ax5[0].set_ylabel('Amplitude')
ax5[0].set_title('ACF')
ax5[0].set_xlabel(r'$\tau$ in $\mu$s')
est1,act_spec=ax5[1].plot(f*1e-3,spec_single.real,f*1e-3,spec_act.real)
ax5[1].legend([est1,act_spec],['Estimated','Actual'])
ax5[1].set_ylabel('Amplitude')
ax5[1].set_title('Spectrum')
ax5[1].set_xlabel(r'f in kHz')
fig5.tight_layout()
Explanation: Plotting and Normalization of Input Spectra
Need to apply window function to spectrum.
End of explanation
<END_TASK>