

In [1]:

    
%matplotlib inline



In [2]:

    
import sys
sys.path.append('../')
from matplotlib import pyplot as plt
import numpy as np



In [3]:

    
from pyphysim.reference_signals.zadoffchu import calcBaseZC, get_shifted_root_seq
from pyphysim.channels.fading import TdlChannel, TdlChannelProfile, COST259_TUx
from pyphysim.channels.fading_generators import JakesSampleGenerator
from pyphysim.channels.multiuser import MuChannel
from pyphysim.util.conversion import linear2dB
from pyphysim.reference_signals.srs import SrsUeSequence
from pyphysim.reference_signals.channel_estimation import CazacBasedChannelEstimator
from pyphysim.reference_signals.root_sequence import RootSequence
#SrsRootSequence,









    



/home/darlan/git_files/pyphysim/pyphysim/util/misc.py:1217: UserWarning: util.misc.count_bits will be slow, since cythonized version was not used
  "util.misc.count_bits will be slow, since cythonized version was not used"



In [4]:

    
num_prbs = 25;        # Number of PRBs to simulate
Nsc = 12 * num_prbs;  # Number of subcarriers
Nzc = 139;            # Size of the sequence
u1 = 25;              # Root sequence index
u2 = u1#12;           # Root sequence index
u3 = u1#7;            # Root sequence index

Note that the sequence size Nzc is lower then the number of subcarriers that will have elements of the Zadoff-Chu sequence. That is $Nzc \leq 300/2 = 150$. Therefore, we will append new elements (creating a cyclic sequence).



In [5]:

    
# Create root sequence objects
a_u1 = RootSequence(u1, size=Nsc//2, Nzc=Nzc)
a_u2 = RootSequence(u1, size=Nsc//2, Nzc=Nzc)
a_u3 = RootSequence(u1, size=Nsc//2, Nzc=Nzc)

Create shifted sequences for 3 users

First we arbitrarely choose some cyclic shift indexes and then we call zadoffchu.getShiftedZF to get the shifted sequence.



In [6]:

    
m_u1 = 1 # Cyclic shift index
m_u2 = 4
m_u3 = 7

r1 = SrsUeSequence(a_u1, m_u1)
r2 = SrsUeSequence(a_u2, m_u2)
r3 = SrsUeSequence(a_u3, m_u3)

Generate channels from users to the BS

Now it's time to transmit the shifted sequences. We need to create the fading channels from two users to some BS.



In [7]:

    
speedTerminal = 3/3.6        # Speed in m/s
fcDbl = 2.6e9                # Central carrier frequency (in Hz)
timeTTIDbl = 1e-3            # Time of a single TTI
subcarrierBandDbl = 15e3     # Subcarrier bandwidth (in Hz)
numOfSubcarriersPRBInt = 12  # Number of subcarriers in each PRB          

# xxxxxxxxxx Dependent parametersxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
lambdaDbl = 3e8/fcDbl        # Carrier wave length
Fd = speedTerminal / lambdaDbl
Ts = 1./(Nsc * subcarrierBandDbl)

# xxxxxxxxxx Channel parameters xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
L = 16      # The number of rays for the Jakes model.



In [8]:

    
# Create the MuSisoChannel
jakes = JakesSampleGenerator(Fd, Ts, L)
musisochannel = MuChannel(N=(1, 3), fading_generator=jakes, channel_profile=COST259_TUx)

Perform the transmission

First we need to prepare the input data from our shifted Zadoff-Chu sequences.

To makes things clear, let's start transmiting a single sequence and we won't include the white noise. Since we use a comb to transmit the SRS sequence, we will use Nsc/2 subcarriers from the Nsc subcarriers from a comb like pattern.



In [9]:

    
comb_indexes = np.arange(0, Nsc, 2)
data = np.vstack([r1.seq_array(),r2.seq_array(),r3.seq_array()])
Y = musisochannel.corrupt_data_in_freq_domain(data, Nsc, comb_indexes)
Y = Y[0]  # We only have one receiver



In [10]:

    
impulse_response0 = musisochannel.get_last_impulse_response(0, 0)
impulse_response1 = musisochannel.get_last_impulse_response(0, 1)
impulse_response2 = musisochannel.get_last_impulse_response(0, 2)

H1 = impulse_response0.get_freq_response(Nsc)[:, 0]
H2 = impulse_response1.get_freq_response(Nsc)[:, 0]
H3 = impulse_response2.get_freq_response(Nsc)[:, 0]

h1 = np.fft.ifft(H1)
h2 = np.fft.ifft(H2)
h3 = np.fft.ifft(H3)

According to the paper,

... the received frequency-domain sequence Y is element-wise multiplied with the complex conjugate of the expected root sequence X before the IDFT. This provides in one shot the concatenated CIRs of all UEs multiplexed on the same root sequence.

Now let's get the plot of the signal considering that all users transmitted. Notice how the part due to user 1 in the plot is the same channel when only user 1 transmitted. This indicates that Zadoff-chu 0 cross correlation is indeed working.



In [17]:

    
y = np.fft.ifft(np.conj(a_u1) * Y, 150)
plt.figure(figsize=(12,6))
plt.stem(np.abs(y), use_line_collection=True)
plt.show()

Estimate the channels

Since we get a concatenation of the impulse response of the different users, we need to know for each users we need to know the first and the last sample index corresponding to the particular user's impulse response.

Since we have Nsc subcarriers, from which we will use $Nsc/2$, and we have 3 users, we can imagine that each user can have up to $Nsc/(2*3)$ samples, which for $Nsc=300$ corresponds to 50 subcarriers.

Now let's estimate the channel of the first user. First let's check again what is the shift used by the first user.



In [12]:

    
m_u1









    Out[12]:





1



In [13]:

    
estimator1 = CazacBasedChannelEstimator(r1)
estimator2 = CazacBasedChannelEstimator(r2)
estimator3 = CazacBasedChannelEstimator(r3)

For an index equal to 1 the starting sample of the first user will be 101 and the ending sample will be 101+50-1=150.



In [24]:

    
def plot_channel_responses(h, tilde_h):
    """Plot the estimated and true channel responses

    Parameters
    ----------
    h : numpy complex array
        The true channel impulse response
    tilde_h : numpy complex array
        The estimated channel impulse response
    """
    H = np.fft.fft(h)
    tilde_H = np.fft.fft(tilde_h, Nsc)

    plt.figure(figsize=(16,12))

    # Plot estimated impulse response
    ax1 = plt.subplot2grid((3,2), (0,0))
    ax1.stem(np.abs(tilde_h[0:20]), use_line_collection=True)
    plt.xlabel("Time sample")
    plt.ylabel("Amplitude (abs)")
    plt.title("Estimated Impulse Response")
    plt.grid()

    # Plot TRUE impulse response
    ax2 = plt.subplot2grid((3,2), (0,1))
    ax2.stem(np.abs(h[0:20]),linefmt='g', use_line_collection=True)
    plt.xlabel("Time sample")
    plt.ylabel("Amplitude (abs)")
    plt.xlabel("Time sample")
    plt.title("True Impulse Response")
    plt.grid()

    # Plot estimated frequency response (absolute value)
    ax3 = plt.subplot2grid((3,2), (1,0), colspan=2)
    plt.plot(np.abs(tilde_H))
    #plt.xlabel("Subcarrier")
    plt.ylabel("Amplitude (abs)")
    plt.title("Frequency Response (abs)")

    # Plot TRUE frequency response (absolute value)
    #plt.subplot(3,2,4)
    ax3.plot(np.abs(H), 'g')
    plt.grid()
    plt.legend(["Estimated Value", "True Value"], loc='upper left')

    # Plot estimated frequency response (angle)
    ax4 = plt.subplot2grid((3,2), (2,0), colspan=2)
    ax4.plot(np.angle(tilde_H))
    plt.xlabel("Subcarrier")
    plt.ylabel("Angle (phase)")
    plt.title("Frequency Response (phase)")

    # Plot TRUE frequency response (angle)
    ax4.plot(np.angle(H),'g')
    plt.grid()
    plt.legend(["Estimated Value", "True Value"], loc='upper left')

    # Show the plots
    plt.show()



In [15]:

    
def plot_normalized_squared_error(H, tilde_H):
    """Plot the normalized squared error (in dB).

    Parameters
    ----------
    H : numpy complex array
        The true channel frequency response
    tilde_H : numpy complex array
        The estimated channel frequency response
    """
    plt.figure(figsize=(12,8))
    error = np.abs(tilde_H - H)**2 / (np.abs(H)**2)
    plt.plot(linear2dB(error))
    plt.title("Normalized Squared Error")
    plt.xlabel("Subcarrier")
    plt.ylabel("Normalized Squared Error (in dB)")
    plt.grid()
    plt.show()



In [25]:

    
# y = np.fft.ifft(np.conj(r1) * Y, 150)
# tilde_h1 = y[0:20]
# tilde_H1 = np.fft.fft(tilde_h1, Nsc)
# tilde_Y1 = tilde_H1[comb_indexes] * r1

tilde_H1 = estimator1.estimate_channel_freq_domain(Y, num_taps_to_keep=20)
tilde_h1 = np.fft.ifft(tilde_H1)[0:20]

plot_channel_responses(h1, tilde_h1)

Now we will compute the squared error in each subcarrier.



In [26]:

    
tilde_H1 = np.fft.fft(tilde_h1, Nsc)
plot_normalized_squared_error(H1, tilde_H1)



In [27]:

    
# y = np.fft.ifft(np.conj(r2) * (Y), 150)
# tilde_h2 = y[0:20]
# tilde_H2 = np.fft.fft(tilde_h2, Nsc)
# tilde_Y2 = tilde_H2[comb_indexes] * r2

tilde_H2 = estimator2.estimate_channel_freq_domain(Y, num_taps_to_keep=20)
tilde_h2 = np.fft.ifft(tilde_H2)[0:20]

plot_channel_responses(h2, tilde_h2)



In [28]:

    
tilde_H2 = np.fft.fft(tilde_h2, Nsc)
plot_normalized_squared_error(H2, tilde_H2)



In [29]:

    
# y = np.fft.ifft(np.conj(r3) * (Y), 150)
# tilde_h3 = y[0:11]
# tilde_H3 = np.fft.fft(tilde_h3, Nsc)
# tilde_Y3 = tilde_H3[comb_indexes] * r3

tilde_H3 = estimator3.estimate_channel_freq_domain(Y, num_taps_to_keep=20)
tilde_h3 = np.fft.ifft(tilde_H3)[0:20]

plot_channel_responses(h3, tilde_h3)



In [30]:

    
tilde_H3 = np.fft.fft(tilde_h3, Nsc)
plot_normalized_squared_error(H3, tilde_H3)

Estimated the channels from corrupted (white noise) signal

Now we will add some white noise to Y



In [31]:

    
# Add white noise
noise_var = 1e-2
Y_noised = Y + np.sqrt(noise_var/2.) * (np.random.randn(Nsc//2) + 1j * np.random.randn(Nsc//2))



In [32]:

    
# y_noised = np.fft.ifft(np.conj(r2) * (Y_noised), 150)
# tilde_h2_noised = y_noised[0:20]

tilde_H2_noised = estimator2.estimate_channel_freq_domain(Y_noised, num_taps_to_keep=20)
tilde_h2_noised = np.fft.ifft(tilde_H2_noised)[0:20]

plot_channel_responses(h2, tilde_h2_noised)



In [33]:

    
tilde_H2_noised = np.fft.fft(tilde_h2_noised, Nsc)
plot_normalized_squared_error(H2, tilde_H2_noised)