Amarisoft

IQ Decode - Python

The purpose of this tutorial is to show how to decode signal from the IQ data. Of course, the purpose is not to provide examples for decoding everything. Decoding every signal and every channel involves the implementation of entire PHY stack which is definitely out of the scope of a tutorial. What I am trying to do in this tutorial is to provide a simple entry point from whichyou can expand and expand.

In this tutorial, I assume that you already have some I/Q data to work with. If you are not familiar with how to capture I/Q, check out following tutorial.

Table of Contents

Introduction

In the field of wireless communications and signal processing, understanding and decoding signals from I/Q (In-phase and Quadrature) data is a foundational skill that enables engineers and researchers to analyze, interpret, and process a wide variety of radio frequency transmissions. I/Q data represents the raw, digitized form of electromagnetic signals captured by software-defined radios (SDRs) or other digital receivers, where the signal is decomposed into its in-phase (I) and quadrature (Q) components. This representation preserves both amplitude and phase information, making it possible to reconstruct the original transmitted waveform and extract meaningful information. The process of decoding signals from I/Q data typically involves several steps, such as filtering, demodulation, symbol decoding, and error correction—each corresponding to layers of the physical (PHY) and data link stack. While the full implementation of a PHY stack for every signal type is a complex endeavor reserved for advanced development, this tutorial provides a practical and accessible entry point for those seeking to begin their journey in signal decoding. By focusing on the fundamentals and essential techniques, learners will gain the foundational knowledge required to interpret I/Q data, paving the way for more advanced exploration and experimentation in the field of digital communications and SDR applications.

Summary of the Tutorial

This tutorial provides a detailed example of Python-based post-processing of LTE downlink I/Q data, specifically captured using Amarisoft SDR hardware. The tutorial includes two main test cases: one for SISO (Single Input Single Output) and one for MIMO 2x2 (Multiple Input Multiple Output), each demonstrating a sequence of signal processing tasks for LTE physical layer analysis and validation. The focus is on step-by-step methodologies for synchronizing, demodulating, and channel estimating LTE signals using captured I/Q data.

Methodology Highlights:

Disclaimer and Special NOTEs!!!

NOTE : The purpose of this tutorial is just to provide an example of Python script for processing a I/Q data. We don't provide any technical support for Python itself.

NOTE : The script explained is not a part of Amarisoft product. It is just to show a proof of concept example showing how you can utilize Amarisoft SDR card and IQ capture functionaly and write your own post processing software. We DO NOTprovide any technical support for this script itself.

NOTE : This script would work with only for specific IQ Data with following requirement. (It means this script is NOT a general purpose decoder. It is writen for a very specific case just to show that a user can utilize Amarisoft I/Q capture data for the input to verify their own decoder algorithm)

NOTE : If you want to try with the script, USE the script that is hyperlinked at the beginning of each test case, DO NOT COPY the code in the script description since it may be a little outdated. You can get the test I/Q data as well at the beginning of each test.

Amarisoft as Study and Development tool for PHY

In this tutorial I wanted to propose a new use case of Amarisoft solution which is a study and development tool for physical layer. Probably in most case, the very beginning of PHY layer development would start from purely software based simulation solutions. At this stage, we usually simulate PHY processing (e.g, generation and decode of PHY signal) purely with software.

However, it eventually reach a point where you need to verify our solution with real signal or real hardware. For example, if you are UE PHY reciever developer, you would need real physical gNB transmitter signal for the test. But generating the proper gNB transmitter signal is not easy (especially when you don't have experties on eNB/gNB transmitter chain) and very tedious job to prepare such a signal just for testing. On top of it, the test signal should come from a well verified solution that has proved to be working. If the test signal is not well proved, the entire test result with your own development become unreliable. So in most case we do this (assuming that we want to develope UE side reciever PHY at baseband)

Based on my experience, I struggled with most of these steps. Some of the difficulties that I had were :

Why Amarisoft ?

Test Setup

You can use any kind of test setup that allow you to capture IQ in spectrum Analyzer mode

Python

Following is the version of the python and package that are used for this tutorial. (I tested this on Windows 11 Home Edition). For python, first try with whatever version you are using as long as it is ver 3.x and upgrade it to latest version if it does not work.

Python Packages

For the script in this tutorial, you need to install following packages.

pip install numpy

pip install matplotlib

pip install scipy

Python Script - Command Line

In this tutorial, I will show you an example of Python script to plot IQ data from a file. This is just an example script intended as POC(Proof of Concept) example. You can revise in any way as you like if this does not suit your purpose.

Test 1 : SISO

By running the script, you would get the test result shown in this section. For the original Python file, click here. For the sample IQ data that you may try with this code, you can download it from here (Unzip the file and change the file name and path in the source code).

Step 1 : Detect LTE PSS and Display (20 Mhz BW)

This step starts with detecting the LTE PSS from captured IQ data and displaying the result for a 20 MHz LTE cell. In this test, the LTE cell is running with 20 MHz bandwidth using enb.default.cfg. The TX signal is captured by another SDR card that is not used by the LTE cell itself. In this example, the SDR card on UEsim is used for capture, but any unused SDR card on the Callbox can also be used. The IQ data is captured with sdr_spectrum using center frequency 2680 MHz, sampling rate 30.72 Msps, capture duration 0.01 second, and output path /tmp. The sampling rate 30.72 Msps matches the standard LTE 20 MHz sampling rate, so the captured IQ file can be directly used for LTE PSS detection and basic frequency-domain analysis.

./sdr_spectrum -c 0 -rx_freq 2680e6 -rate 30.72e6 -save-and-exit -duration 0.01 -save_path "/tmp"

To capture the signal with better quality, you can add rx_gain and save_delay to the sdr_spectrum command. In this example, rx_gain is set to -1, which lets the SDR use automatic gain control instead of forcing a fixed gain value. This can help avoid saturation when the received signal is strong. The save_delay is set to 5 seconds, so the tool waits for a short time before saving the IQ data. This gives the receiver more time to settle before the actual capture starts. The rest of the settings are the same as before. The center frequency is 2680 MHz, the sampling rate is 30.72 Msps, the capture duration is 0.01 second, and the IQ file is saved under /tmp.

./sdr_spectrum -c 0 -rx_freq 2680e6 -rate 30.72e6 -rx_gain -1 -save-and-exit -duration 0.01 -save_delay 5 -save_path "/tmp"

In this test, the main routine is to detect LTE PSS from the captured IQ data. The DetectPSS() function does this by sliding through the IQ samples one sample at a time. For each position, it takes one IQ chunk whose size depends on the LTE channel bandwidth. This chunk is still in time domain, so FFT is applied first to convert it into frequency domain. Then the subcarriers around the PSS location are extracted. The extracted PSS area is compared with the three possible LTE PSS sequences defined by 3GPP. These three sequences correspond to the three possible N_ID_2 values. For each IQ position, the function calculates the correlation with all three generated PSS sequences and selects the largest correlation result. If this result is greater than the current maximum value, the function updates the maximum value and saves this position as the best PSS candidate. By repeating this process while moving forward one IQ sample at a time, the function finds the time position where the received signal matches the generated LTE PSS sequence most strongly.

IQ Decode LTE PSS Ex01 FrequencyDomainDetection

When the code is executed, the first result is the text summary of the PSS detection process. In this example, the strongest PSS correlation is found for N_ID_2 = 0. The maximum correlation value is about 1.685, and this is higher than the correlation values for PSS1 and PSS2. This means that the received signal matches the generated PSS0 sequence most strongly. The maximum correlation position is 269166, which indicates the sample position where the PSS match is strongest in the captured IQ data. The estimated frequency offset is about 127.12 kHz, so the received signal is not exactly centered at the expected frequency and the code needs to compensate for this offset before building the resource grid.

After PSS detection, the code also reports the resource grid size as 14 by 1199. This means one LTE subframe is represented with 14 OFDM symbols in time and 1199 subcarriers in frequency. The detected subframe is 0. Then SSS detection gives N_ID_1 = 84. Since LTE PCI is calculated from N_ID_1 and N_ID_2 as PCI = 3 x N_ID_1 + N_ID_2, the detected PCI becomes 252. In this example, N_ID_1 is 84 and N_ID_2 is 0, so the final detected PCI is 252.

Detected PSS Sequence (N_ID_2): 0

Maximum correlation: 1.6852206708698385

Maximum correlation for PSS0: 1.6852206708698385

Maximum correlation for PSS1: 1.2416934348611477

Maximum correlation for PSS2: 1.1363440216733818

Position of maximum correlation: 269166

Estimated Frequency Offset: 127120.00185358223 Hz

Dimension of resourceGrid = (14, 1199)

Detected Subframe : 0

Detected N_ID_2 (PSS) : 0

Detected N_ID_1 (SSS) : 84

Detected PCI: 252

After PSS detection is completed, ShowPSS() displays the detection result in graphical form. The first plot shows the frequency spectrum of the captured IQ samples. Since this test uses a 20 MHz LTE cell with 30.72 Msps sampling rate, the occupied LTE signal appears as a wide flat spectrum in the middle of the capture bandwidth. The lower plots on the left show the detected PSS constellation, the generated reference PSS constellation for N_ID_2 = 0, and the time-domain IQ waveform with the detected PSS position marked. This makes it easier to confirm that the detected PSS is not only a text result, but also visible in the captured signal.

The small plots on the right show the resource grid after synchronization and frequency correction. Each row represents one OFDM symbol in the LTE subframe. The strong narrow peaks around the center subcarriers indicate synchronization signals and other center-frequency components. In symbol 5 and symbol 6, the strong center peaks correspond to the detected PSS/SSS area. Other symbols show wider occupied bandwidth when normal LTE downlink data or reference signals are present. So this figure is useful as a quick visual check. It confirms that the captured IQ file contains a valid LTE downlink signal, the center frequency offset was estimated, the PSS was detected, and the received samples were converted into an LTE resource-grid style view.

IQ Decode LTE PSS Ex01 01

(A) : This is frequency spectrum (plot in frequency) of the entire IQ data given to the program.

(B) : This shows I/Q data in 2D parameteric plot detected from the IQ data. The blue plot shows the I/Q data with magnitue normalized to 1. The plot in orange plot shows the detected PSS sequence with magnitude of all data mapped to 1 reserving angle as it is.

(C) : This shows the expected I/Q data (i.e, The sequence generated by 3GPP spec) corresponding to the detected PSS sequence.

(D) : This is time doman plotof the entire IQ data given to the program. The orange part is the where the detected PSS is located. You can magnify a specific portion of the plot using toolbox at the top of the window.

(0)~(13) : Shows the frequency spectrum for each OFDM symbols within the subframe. (6) is the OFDM symbol where PSS is located.

After PSS detection, ShowGrid() displays the received LTE signal in resource-grid format. In this plot, the horizontal axis is subcarrier index and the vertical axis is OFDM symbol index. Since this is a normal LTE subframe, there are 14 OFDM symbols from Symbol 0 to Symbol 13. The color represents the received power of each resource element. A brighter or redder area means stronger received power, and a darker blue area means weaker or almost empty power.

In this example, the center area around subcarrier 600 is clearly visible. This is expected because the LTE synchronization signals are located around the center of the LTE carrier. The plot marks SSS on Symbol 5 and PSS on Symbol 6. This matches the LTE FDD synchronization signal position in subframe 0. Around this center region, strong power is observed because PSS, SSS, PBCH, and related center-channel signals are located there. Other symbols show different power patterns across the bandwidth depending on reference signals, control/data allocation, and unused resource elements.

This view is useful because it converts raw IQ samples into a form that is closer to the LTE physical resource grid. From this plot, we can visually confirm that the PSS and SSS were detected at the expected location, the center frequency alignment is reasonable, and the captured IQ data contains a valid LTE downlink signal.

IQ Decode LTE PSS Ex01 02

Step 2: Detect LTE SSS and Display (20 Mhz BW)

This step detects LTE SSS from the resource grid that was already time aligned by the detected PSS. First, the script extracts the SSS area from the resource grid. Since the PSS position was already found in Step 1, the code can estimate where the SSS should be located in the same subframe. Then the extracted SSS samples are compensated with the frequency offset estimated during PSS detection. In this example, the SSS data after this compensation is treated as SSS_r.

After frequency compensation, the script generates all possible LTE SSS sequences for the detected N_ID_2 value. At this point, N_ID_2 is already known from PSS detection, so the remaining unknown value is N_ID_1. The script generates 168 possible SSS candidates, corresponding to N_ID_1 values from 0 to 167. For each N_ID_1, a pair of SSS sequences is generated because LTE SSS has two possible sequence forms depending on the subframe. Then the received SSS is correlated with these generated SSS candidates. The sequence that gives the highest correlation is selected as the detected SSS result.

This step completes the PCI detection. PSS gives N_ID_2, and SSS gives N_ID_1. Once both values are known, the LTE PCI can be calculated as 3 x N_ID_1 + N_ID_2. So the purpose of this step is not only to detect SSS itself, but also to finalize the LTE cell identity from the captured IQ data.

IQ Decode LTE SSS Ex02 DetectionAlgorithm

After PSS detection, the script continues with SSS detection and displays the result using ShowSSS(). This figure compares the received SSS before correction, the received SSS after frequency correction, and the generated reference SSS for the detected N_ID_1. In plot A, the SSS constellation before correction is not well aligned. This is because the captured IQ still includes the frequency offset estimated during PSS detection. Due to this offset, the phase of the SSS samples rotates and the constellation does not match the expected reference sequence clearly.

In plot B, the same SSS samples are shown after frequency correction. After applying the estimated frequency offset compensation, the SSS constellation becomes much better aligned. The real and imaginary sample patterns also become more similar to the generated SSS sequence. This shows that the frequency offset estimated from PSS is useful not only for PSS detection, but also for improving the following SSS detection step.

In plot C, the generated reference SSS sequence for the detected N_ID_1 is shown. The corrected received SSS is compared with this generated sequence, and the sequence that gives the highest correlation is selected as the detected SSS. From this result, the script obtains N_ID_1. Together with N_ID_2 from PSS detection, the final LTE PCI can be calculated. So this plot is a visual confirmation that the received SSS was corrected properly and matched against the correct reference SSS sequence.

IQ Decode LTE SSS Ex01

Step 3: Detect CRS, Estimate Channel Coefficientand Display (20 Mhz BW)

This step detects LTE CRS and uses it to estimate the channel coefficient. Once the PCI is detected correctly, the CRS resource element location can be calculated inside the LTE resource grid. CRS position depends on the PCI, antenna port, OFDM symbol, and subcarrier index. So after PSS and SSS detection, the script already has the most important information required to find the CRS locations.

First run crs_rx, crs_ex, h, crs_k = GetCRSandChannelCoef(resourceGrid, PCI, nRB = 100, method = "generate", subframe=0, debug=False) and run following functions for visualization. Followings are description of arguments and returns.

GetCRSandChannelCoef() takes the resource grid and PCI as input. In this example, nRB is set to 100, which means 20 MHz LTE bandwidth. The method is set to generate, so the expected CRS sequence is generated from the 3GPP rule. The script then extracts the received CRS values from the resource grid and compares them with the expected CRS values. For a simple SISO case and low-noise condition, the channel coefficient can be estimated by dividing the received CRS value by the expected CRS value.

resourceGrid is the 2D array that contains the received IQ values for one LTE subframe. PCI is the detected physical cell ID from Step 2. nRB is the number of LTE resource blocks, and 100 means 20 MHz LTE. method controls how the expected CRS is prepared. generate creates the CRS sequence from the 3GPP specification, while testvector uses a predefined CRS reference for PCI 252, nRB 100, and subframe 0. subframe indicates the subframe number within the radio frame. debug controls whether debug messages are printed.

The function returns four main results. crs_k stores the resource element index of CRS on OFDM symbols 0, 4, 7, and 11. crs_rx stores the received CRS values extracted from the resource grid. crs_ex stores the expected CRS values generated by the script. h stores the estimated channel coefficient, calculated from crs_rx and crs_ex. This gives a simple channel estimate at the CRS locations and becomes the basis for checking channel quality, phase rotation, amplitude distortion, and possible equalization.

Followings are constellations and I/Q sequence plot for CRS (Cell Reference Signal) from the detected subframe. (In the constellation, the solid blue indicates the CRS at the lowest frequency and the solid red indcates the CRS at the highest frequency). This is the output of ShowCRS().

ShowCRS() displays the CRS result from the detected LTE subframe. In this example, the CRS is extracted from OFDM symbols 0, 4, 7, and 11, which are the CRS-bearing symbols for this LTE configuration. Each row of the figure corresponds to one of these OFDM symbols.

Column A shows the constellation of the received CRS values extracted from resourceGrid. These are the actual CRS samples captured from the IQ data. The points do not stay at the ideal QPSK positions because they include the radio channel effect, frequency offset residue, phase rotation, gain variation, noise, and possible RF distortion. In the constellation, the solid blue point indicates the CRS at the lowest frequency side, and the solid red point indicates the CRS at the highest frequency side. This coloring helps show how the CRS phase and amplitude change across the channel bandwidth.

Columns B and C show the I and Q sequences of the received CRS. These plots are another way to inspect the same CRS samples in sequence order. Instead of looking at the constellation shape, these plots show how the real and imaginary parts vary across CRS resource elements. If the received signal has frequency-selective fading or phase rotation across the bandwidth, it can be observed as a changing pattern in these I and Q sequences.

Column D shows the expected CRS constellation generated by the 3GPP algorithm. Since CRS uses QPSK-like reference symbols, the generated CRS points appear at clean reference positions. These are not received samples. They are the ideal CRS values that the receiver expects for the detected PCI, subframe number, and CRS position.

Column E shows the estimated channel coefficient. It is calculated by dividing the received CRS value by the expected CRS value, h = crs_rx / crs_ex. This removes the known CRS modulation and leaves the estimated channel response at the CRS locations. The resulting constellation shows how the channel changes over frequency. In this example, the channel coefficient forms a curved shape rather than a single fixed point, which means the received signal has frequency-dependent phase and amplitude variation across the LTE bandwidth. This is the basic idea of CRS-based channel estimation. The known CRS is used as a reference, and the difference between the received CRS and expected CRS is interpreted as the radio channel.

IQ Decode LTE CRS Ex01

Column (A) : The constellation of CRS retrieved from resourceGrid (these are the received CRS). Rows indicates symbol 0, 4, 7, 11

Column (B), (C): I and Q sequence of the CRS retrieved from resourceGrid. Rows indicates symbol 0, 4, 7, 11

Column(D) : The constellation of CRS generated by 3GPP algorithm. Rows indicates symbol 0, 4, 7, 11

Column(E) : The constellation of channel coefficient estimated from (A) and (D). Rows indicates symbol 0, 4, 7, 11

Step 4: Creation of Channel Coefficient Grid

At this step, I will create a resource grid filled with channel coefficient at each and every resource elements. This is done in a sequence of multiple steps listed below.

This step creates a channel coefficient grid. In previous, the channel coefficient was estimated only at the CRS resource elements. However, for equalization or further decoding, it is more useful to have a channel coefficient value for every resource element in the LTE resource grid.

So this step first creates a new grid called hResourceGrid with the same dimension as resourceGrid. All elements are initialized with 0.0 + j0.0. Then the estimated CRS-based channel coefficients are inserted into this grid at the correct CRS positions. For OFDM symbol 0, the CRS subcarrier positions are stored in crs_k[0], and the corresponding channel coefficients are stored in h[0]. These values are written into symbol 0 of hResourceGrid. The same operation is repeated for OFDM symbols 4, 7, and 11 using crs_k[1]/h[1], crs_k[2]/h[2], and crs_k[3]/h[3]. After this operation, hResourceGrid contains valid channel coefficient values only at the CRS locations, while all other resource elements are still zero.

This step does not yet estimate the channel coefficient for every resource element. It only places the known CRS-based channel estimates into the correct positions inside the full resource grid. The remaining empty locations will be filled later by interpolation in time and frequency direction.

Then the channel coefficient grid is expanded from CRS-only positions to the full LTE resource grid by interpolation. At the beginning of this step, hResourceGrid has channel coefficient values only at the CRS resource elements on OFDM symbols 0, 4, 7, and 11. All other resource elements are still empty. The interpolateGrid() function fills these empty resource elements using the CRS-based channel coefficient values as reference points.

In this function, frequency-domain interpolation is performed first. This means the channel coefficient values between CRS subcarriers are estimated along the subcarrier direction within each CRS-bearing OFDM symbol. After that, time-domain interpolation is performed. This means the channel coefficient values for the OFDM symbols between the CRS-bearing symbols are estimated using the already interpolated channel values from nearby symbols. As a result, hResourceGrid_interpolated becomes a full channel coefficient grid where every resource element has an estimated channel value. In this example, hResourceGrid is the CRS-based channel coefficient grid created in the previous step. mode = 0 selects the interpolation mode used by the script. td_method = "interpolate" means the missing channel coefficients in the time direction are estimated by interpolation rather than simple copy or fixed assignment. direction = "ccw" specifies the direction used by the interpolation logic. After this step, the generated channel coefficient grid can be used for visualization and also for basic channel compensation or equalization of the received LTE resource grid.

Followings are the descriptions on arguments and returns :

Then plotGridImage() is used to visualize the interpolated channel coefficient grid. In this example, the input is hResourceGrid_interpolated, which already has channel coefficient values for every resource element after frequency-domain and time-domain interpolation. The plot shows the magnitude of the complex channel coefficient at each resource element. The horizontal axis is the subcarrier index, and the vertical axis is the OFDM symbol index from Symbol 0 to Symbol 13.

The plot shown below is the output of plotGridImage(hResourceGrid_interpolated). This shows the magnitude of complex numbers (channel coefficient) in every resource elements of the resource grid.

The color indicates the magnitude of the estimated channel coefficient. A brighter color means a larger channel gain, and a darker color means a smaller channel gain. In this result, the channel magnitude changes mainly along the frequency direction. This means the estimated channel response is not flat across the 20 MHz bandwidth. Some subcarrier regions have stronger channel gain, while other regions have weaker gain. The variation is smoothly extended across all OFDM symbols because the CRS-based estimates were interpolated over the whole resource grid.

This plot is useful because it gives a full-subframe view of the estimated channel response. Instead of checking only CRS positions, we can now see the estimated channel coefficient over all resource elements. This channel grid can be used later for equalization by compensating each received resource element with its corresponding estimated channel coefficient.

IQ Decode LTE H Grid Ex01

Output of ShowResourceGridIQ(hResourceGrid_interpolated, figTitle="Interpolated [h] Grid I/Q", figTitleColor="green").

This output shows the interpolated channel coefficient grid in constellation form. ShowResourceGridIQ() plots the complex value of hResourceGrid_interpolated for each OFDM symbol. Each small plot corresponds to one OFDM symbol, from symbol 0 to symbol 13. The horizontal axis is the real part of the channel coefficient, and the vertical axis is the imaginary part.

Unlike the previous image plot, this view shows both magnitude and phase of the estimated channel coefficient. The points form a curved arc rather than a compact cluster. This means that the channel coefficient changes its phase across the subcarriers. Since the points also have different distance from the origin, the channel gain also changes across frequency. So this plot gives a more detailed view of the channel response than the magnitude-only image.

The shape is similar across most OFDM symbols because the full h grid was created by interpolation from CRS symbols 0, 4, 7, and 11. Symbols that are close to CRS-bearing symbols follow the CRS-based estimate more directly, while intermediate symbols are filled by time-domain interpolation. This is why the constellation shape changes smoothly from symbol to symbol.

This figure is useful for checking whether the interpolated channel coefficient looks continuous and reasonable over the whole subframe. If the interpolation works properly, the channel coefficient should not jump randomly between neighboring OFDM symbols. In this example, the channel coefficient pattern is mostly continuous, so the interpolated h grid can be used as the channel estimate for later equalization.

IQ Decode LTE H Grid Ex02

Step 5: Equalization

This step performs equalization using the interpolated channel coefficient grid. At this point, the received LTE resource grid is already created from the captured IQ data, and the channel coefficient grid hResourceGrid_interpolated has already been estimated from CRS and expanded to all resource elements. So each received resource element has a corresponding estimated channel coefficient.

The basic equalization is done by compensating the received resource element with the estimated channel coefficient. In simple form, the equalized value is calculated as received resource element divided by channel coefficient. This removes the estimated channel effect from the received signal. After this operation, phase rotation and amplitude distortion caused by the channel should be reduced, and the constellation of the received symbols should become closer to the expected modulation points.

This step is important because the raw received resource grid still includes the effect of the radio channel. Even if the transmitted signal uses clean QPSK, 16QAM, or 64QAM symbols, the received constellation may look rotated, stretched, or curved due to the channel. By using hResourceGrid_interpolated, the script applies channel compensation over the whole resource grid. The result can then be used for visual inspection of the equalized constellation and for further decoding steps.

Before equalization, let's plot the constellation of the received signal as it is using ShowResourceGridIQ(resourceGrid, clip = 20, figTitle="Recieved Resource Grid I/Q : Before Equalization", figTitleColor="green")

ShowResourceGridIQ() is used to plot the received resource grid as it is. This means the constellation is drawn directly from resourceGrid without applying channel compensation. Each small plot corresponds to one OFDM symbol from symbol 0 to symbol 13. The horizontal axis is the I component, and the vertical axis is the Q component. The clip value is set to 20, so very large values are limited in the plot range to make the constellation easier to inspect.

In this result, the points do not form clean modulation clusters. Instead, they are widely spread and rotated in a circular or spiral-like pattern. This is expected because the received resource grid still includes the radio channel effect. The channel changes both amplitude and phase across the subcarriers, so the received symbols are distorted before equalization. Some OFDM symbols also include different LTE physical signals and unused resource elements, so the shape is not exactly the same for every symbol.

This plot is useful as a before-equalization reference. It shows how the raw received LTE resource grid looks when no channel correction is applied. After equalization using hResourceGrid_interpolated, the same type of plot can be generated again. By comparing the before and after plots, we can visually check whether the channel compensation reduces the phase rotation and makes the received constellation closer to the expected LTE modulation points.

IQ Decode LTE Equalization Ex01

Then the received resource grid is equalized by dividing it with the interpolated channel coefficient grid. This is done by running ResourceGridDiv() with resourceGrid and hResourceGrid_interpolated as inputs. The operation is applied to every resource element in the grid.

This means that each received IQ value is divided by the estimated channel coefficient at the same resource element position.

The purpose of this operation is to remove the estimated channel effect from the received signal. Before this step, each resource element includes both the transmitted signal and the channel distortion. After division by the channel coefficient, the phase rotation and amplitude variation caused by the channel are reduced. The output, resourceGrid_corrected, becomes the equalized resource grid. This corrected grid can then be visualized again in constellation form to check how much the signal quality improved after equalization.

After equalization, plot the result of the operation using ShowResourceGridIQ(resourceGrid_corrected, clip = 10, figTitle="Corrected Resource Grid I/Q : After Equalization", figTitleColor="green")

ShowResourceGridIQ() is used again to display the corrected resource grid. In this plot, the input is resourceGrid_corrected, which is the result of dividing the received resource grid by hResourceGrid_interpolated. The clip value is set to 10, so the constellation range is smaller than the before-equalization plot and the corrected modulation points can be seen more clearly.

Compared with the before-equalization result, the constellation is much more organized. The previous circular or spiral-like distortion is mostly removed, and the points now form a regular grid-like pattern. This means that the CRS-based channel estimate successfully compensated for much of the phase rotation and amplitude variation caused by the radio channel. Each OFDM symbol still shows some spreading because the signal includes noise, interpolation error, residual frequency/phase error, and different LTE physical channels. But overall, the corrected constellation is much closer to the expected LTE modulation structure.

This plot is the main visual confirmation of the equalization step. Before equalization, the received IQ points were strongly distorted by the channel. After equalization, the points become aligned into clearer modulation clusters. This shows that the detected CRS, estimated channel coefficient grid, and interpolation-based channel compensation are working reasonably well for this captured 20 MHz LTE downlink signal.

IQ Decode LTE Equalization Ex02

Test 2: MIMO 2x2

This test is to show on how to decode 2x2 MIMO IQ data. In terms of fundamental algorithm, the log is same as SISO case. But the complexity of implementation, especially in terms of channel coefficient and equalization. In case of SISO, the channel coefficient is 1x1 (e.g, complex scalar), but in 2x2 MIMO they becomes 2x2 MIMO.

By running the script, you would get the test result shown in this section. For the original Python file, click here. For the sample IQ data that you may try with this code, you can download it from here (Unzip the file and change the file name and path in the source code).

Step 1: Detect LTE PSS/SSS and Display (20 Mhz BW) for Antenna 0

This step detects LTE PSS and SSS using the IQ data from antenna 0. In this implementation, synchronization and PCI detection are performed only on antenna 0. The detected timing, frequency offset, N_ID_2, N_ID_1, and PCI are then reused when processing antenna 1. This is reasonable because both antenna 0 and antenna 1 are captured from the same LTE cell and should share the same frame timing and cell identity.

The script for this step is as follows.

The Amarisoft Remote API capture creates two separate IQ files for the 2x2 MIMO signal. One file contains the IQ samples for antenna 0, and the other file contains the IQ samples for antenna 1. The script reads these two files separately and stores them into iq_samples_p0 and iq_samples_p1. In this test, the IQ data was captured from UEsim. For LTE 20 MHz, the UEsim capture sample rate is fixed to 23.04 MHz. Even though this data can still contain the 20 MHz LTE signal, the script resamples it to 30.72 MHz because 30.72 Msps is the standard LTE 20 MHz sampling rate and fits better with the FFT size and resource-grid construction used in this tutorial.

After resampling, the script sets sample_rate to 30.72e6, fft_bin_size to 2048, and NRB to 100. NRB 100 means LTE 20 MHz bandwidth. Then DetectPSS() is executed using iq_samples_p0 only. The detected PSS result is passed to ShowPSS(), which displays the spectrum, PSS constellation, detected PSS position, and resource-grid related view for antenna 0. At this stage, the main purpose is to get reliable synchronization and cell identity information from antenna 0 before applying the same timing reference to antenna 1.

Following plot is the ShowPSS() output for antenna 0 in the 2x2 MIMO test. The left side shows the captured IQ signal and the detected PSS result. The frequency spectrum shows a 20 MHz LTE signal after the antenna 0 IQ data was resampled to 30.72 Msps. The occupied LTE bandwidth is clearly visible in the middle of the spectrum, so the captured data contains the expected LTE downlink signal.

The PSS constellation plot compares the detected PSS samples with the normalized reference PSS. The detected PSS is not perfectly aligned with the reference because it still includes channel effect, residual frequency error, phase rotation, and noise. The generated PSS plot shows the ideal PSS sequence for the detected N_ID_2. In this example, N_ID_2 is detected as 0. The time-domain plot at the bottom marks the detected PSS position in the captured antenna 0 IQ stream.

The plots on the right show the OFDM-symbol level spectrum after synchronization. Each row corresponds to one OFDM symbol. The LTE signal appears across the expected bandwidth, and the center area shows the synchronization-related components. This confirms that antenna 0 has been synchronized properly and can be used as the timing reference for the rest of the MIMO processing. After this step, the same timing and PCI information from antenna 0 can be reused to construct the resource grid for antenna 1 as well.

IQ Decode LTE Test2 PSS 01

Step 2: Construct the Resource Gridfor Antenna 0 and Antenna 1

This step constructs the LTE resource grid for both antenna 0 and antenna 1. The synchronization information was already obtained from antenna 0 in the previous step. This includes PSS timing, frequency offset, N_ID_2, N_ID_1, and PCI. The same information is then used to process both antenna IQ streams because both antennas are receiving the same LTE cell.

AnalyzeGrid() is called once for antenna 0 and once for antenna 1. For antenna 0, it uses ofdm_symbols_p0. For antenna 1, it uses ofdm_symbols_p1. In both cases, the same PSS detection result and antenna port information are used. The function first extracts the required parameters from the PssInfo dictionary and the input arguments. Then it calls getResourceGrid() to build the LTE resource grid from the OFDM symbols. After the grid is created, it detects SSS using DetectSSS() and prints the detected subframe, N_ID_2, N_ID_1, and PCI.

AnalyzeGrid() also performs CRS-based channel processing. It calls GetCRSandChannelCoef() to locate the CRS resource elements, extract the received CRS, generate the expected CRS, and estimate the channel coefficient. Then it creates a channel coefficient resource grid, places the CRS-based channel estimates at the CRS positions, and interpolates the missing values over the whole grid. After that, the function equalizes the received resource grid using the interpolated channel coefficient grid. It also checks for invalid values such as NaN or Inf and clips them if needed.

So AnalyzeGrid() is not only a display function. It is the main processing function for each antenna branch. It constructs the received resource grid, estimates the CRS-based channel coefficient, creates the interpolated channel grid, performs equalization, and returns all important results in a dictionary. For the 2x2 MIMO case, this same operation is applied separately to antenna 0 and antenna 1, so each receive antenna gets its own resource grid, CRS result, channel coefficient grid, and corrected resource grid.

Actually this function is an important one and it does many other things in addition to display the resource grid. High level functionality of AnalyzeGrid() function is as follows. (For more details, look into the source code).

Following two plots show the constructed LTE resource grid for antenna 0 and antenna 1. Both grids are generated using the synchronization and PCI information detected from antenna 0. The horizontal axis is the subcarrier index, and the vertical axis is the OFDM symbol index. Since this is a 20 MHz LTE signal, the grid spans about 1200 subcarriers, corresponding to 100 RBs. Each grid has 14 OFDM symbols for one LTE subframe.

The upper plot is the resource grid for antenna 0, and the lower plot is the resource grid for antenna 1. Both plots show a similar LTE downlink structure, which is expected because both antennas are receiving the same LTE cell. The center region around subcarrier 600 is clearly visible in symbols 5 and 6. This is where SSS and PSS are located in the detected subframe. The fact that the center synchronization area appears in both antenna grids confirms that the timing detected from antenna 0 can also be applied to antenna 1.

Even though the overall structure is similar, the detailed color pattern is not exactly the same between antenna 0 and antenna 1. This is normal in a MIMO capture. Each receive antenna observes the same transmitted LTE signal through a different RF path and channel condition. Therefore, the received power, phase, fading pattern, and noise can be different for each antenna. These two resource grids become the starting point for separate CRS extraction and channel estimation for each receive antenna branch.

IQ Decode LTE Test2 ResourceGrid 01

IQ Decode LTE Test2 ResourceGrid 02

Step 3: Construct channel matrix

This step constructs the 2x2 MIMO channel matrix. In the SISO case, the channel coefficient for each resource element is just one complex value, so equalization can be done by simple division. In the 2x2 MIMO case, the receiver has two antenna inputs and the transmitter uses two CRS antenna ports. Therefore, each resource element needs a 2x2 channel matrix instead of one scalar channel coefficient.

The channel matrix is built from the CRS-based channel estimates from both received antenna grids. h11 is the channel coefficient estimated from antenna port 0 CRS received on antenna 0. h21 is the channel coefficient estimated from antenna port 0 CRS received on antenna 1. h12 is the channel coefficient estimated from antenna port 1 CRS received on antenna 0. h22 is the channel coefficient estimated from antenna port 1 CRS received on antenna 1.

So the channel matrix for each resource element can be interpreted as follows.

This means the received signal at each antenna is a mixture of the signals transmitted from both antenna ports. The purpose of the channel matrix is to describe this mixture. Once this 2x2 channel matrix is estimated for each resource element, MIMO equalization can be performed by using matrix-based compensation instead of scalar division.

Following figure illustrates how the 2x2 MIMO channel matrix is constructed from CRS. In SISO, one received CRS value can be compared with one expected CRS value, and the result becomes one channel coefficient. In 2x2 MIMO, there are two transmit antenna ports and two receive antennas, so the receiver needs four channel coefficients.

The red CRS pattern represents CRS from antenna port 0, and the blue CRS pattern represents CRS from antenna port 1. These CRS resource elements are transmitted at different RE positions, so the receiver can separate which CRS belongs to which antenna port. Antenna 0 receives both port 0 CRS and port 1 CRS. Antenna 1 also receives both port 0 CRS and port 1 CRS. From these received CRS values, the script estimates h11, h12, h21, and h22.

h11 is estimated from port 0 CRS received at antenna 0. h21 is estimated from port 0 CRS received at antenna 1. h12 is estimated from port 1 CRS received at antenna 0. h22 is estimated from port 1 CRS received at antenna 1. These four values form the 2x2 channel matrix. This matrix explains how the two transmitted symbols x1 and x2 are mixed together and received as y1 and y2 at the two receive antennas.

The equation in the figure shows this relationship. y1 and y2 are the received signals at antenna 0 and antenna 1. x1 and x2 are the transmitted signals from antenna port 0 and antenna port 1. The channel matrix contains h11, h12, h21, and h22. n1 and n2 represent noise. So for each resource element, the receiver sees a mixture of the two transmitted signals through four different channel paths. This is why 2x2 MIMO equalization cannot be done by simple scalar division. It needs matrix-based equalization using the estimated 2x2 channel matrix.

IQ Decode LTE Test2 Hmatrix 01

image source : sharetechnote

Once the channel coefficient matrix is extracted at the CRS positions, the next step is to create four channel coefficient resource grids for visualization and later MIMO equalization. In 2x2 MIMO, the channel is represented by four components: h11, h21, h12, and h22. Each one is created as a separate resource grid. h11 represents the channel from TX port 0 to RX antenna 0. h21 represents the channel from TX port 0 to RX antenna 1. h12 represents the channel from TX port 1 to RX antenna 0. h22 represents the channel from TX port 1 to RX antenna 1.

First, the script processes antenna 0 with AntPort = 0. AnalyzeGrid() extracts the antenna 0 resource grid, detects or confirms the cell information, extracts CRS for port 0, generates the expected CRS, calculates the CRS-based channel coefficient, creates the h resource grid, and interpolates it. The result is saved as h11 and h11_interpolated. h11 contains channel coefficient values only at CRS positions before interpolation. h11_interpolated contains the channel coefficient values after interpolation over the full resource grid.

Then the script processes antenna 1 with AntPort = 1. This produces resourceGrid_p1 and the channel coefficient grid for port 1 observed on antenna 1. The result is saved as h22 and h22_interpolated. In the same way, h22 contains the CRS-only channel coefficient grid before interpolation, and h22_interpolated contains the full interpolated channel coefficient grid.

After h11 and h22 are obtained directly from AnalyzeGrid(), the script creates the cross-channel terms h21 and h12. h21 is calculated by using the CRS positions and expected CRS values for port 0, but extracting the received CRS values from antenna 1 resource grid. In other words, the script checks how the port 0 CRS is received at antenna 1. For each CRS-bearing OFDM symbol 0, 4, 7, and 11, getREs() extracts the received resource elements from resourceGrid_p1 at the port 0 CRS locations. Then the received values are divided by the expected port 0 CRS values. This gives the h21 channel coefficient at CRS positions. These values are inserted into a blank h21 grid, and then interpolateGrid() creates h21_interpolated.

h12 is created in the opposite way. The script uses the CRS positions and expected CRS values for port 1, but extracts the received CRS values from antenna 0 resource grid. This checks how the port 1 CRS is received at antenna 0. Again, CRS values are extracted at OFDM symbols 0, 4, 7, and 11, divided by the expected CRS values, inserted into a blank h12 grid, and interpolated to create h12_interpolated.

After this step, we have four full channel coefficient grids: h11_interpolated, h21_interpolated, h12_interpolated, and h22_interpolated. These four grids form the 2x2 MIMO channel matrix for every resource element. This is the key difference from the SISO case. In SISO, one interpolated h grid is enough. In 2x2 MIMO, four interpolated h grids are required because each received antenna contains a mixture of signals from both transmit ports.

Now the four channel coefficient grids are plotted before interpolation. At this stage, h11, h21, h12, and h22 contain valid channel coefficient values only at the CRS resource element positions. The other resource elements are still empty. This is why the plot shows sparse horizontal CRS patterns instead of a fully filled resource grid.

gridList contains the four channel coefficient grids, and titleList defines the label for each plot. plotGridTile() then displays all four grids together so that we can compare the channel coefficient pattern for each MIMO path.

h11 shows the channel from TX port 0 to RX antenna 0. h21 shows the channel from TX port 0 to RX antenna 1. h12 shows the channel from TX port 1 to RX antenna 0. h22 shows the channel from TX port 1 to RX antenna 1. In each plot, the bright points appear only on CRS-bearing OFDM symbols and CRS-bearing subcarriers. This confirms that the CRS extraction and channel coefficient placement were done at the expected resource element positions.

This plot is mainly a CRS-position check before interpolation. It verifies that all four MIMO channel components were created correctly and that each h grid has channel estimates at the proper CRS locations. After this, interpolation is applied to each grid so that h11, h21, h12, and h22 can have channel coefficient values for every resource element in the LTE resource grid.

IQ Decode LTE Test2 Hmatrix 02

Then the same interpolation algorithm used in the SISO case is applied to each MIMO channel coefficient grid. Before interpolation, h11, h21, h12, and h22 have valid channel coefficient values only at CRS positions. After interpolation, each grid is filled over the whole LTE resource grid, so every resource element has an estimated channel coefficient for each MIMO path.

In this step, h11_interpolated, h21_interpolated, h12_interpolated, and h22_interpolated are collected into gridList and displayed together with plotGridTile(). Each plot shows the magnitude of one interpolated channel coefficient grid. h11_interpolated represents TX port 0 to RX antenna 0. h21_interpolated represents TX port 0 to RX antenna 1. h12_interpolated represents TX port 1 to RX antenna 0. h22_interpolated represents TX port 1 to RX antenna 1.

Compared with the previous CRS-only plot, the resource grids are now fully filled. The color changes smoothly across subcarriers and OFDM symbols because the missing channel coefficients were estimated by frequency-domain and time-domain interpolation. The detailed pattern is different for each h component, which is expected in 2x2 MIMO because each transmit-to-receive antenna path experiences a different channel response.

These four interpolated grids now form the full 2x2 channel matrix for every resource element. This is the channel information required for MIMO equalization. Instead of using one scalar channel coefficient as in SISO, the receiver can now use the four values h11, h21, h12, and h22 at each resource element to compensate the mixed signals received on antenna 0 and antenna 1.

IQ Decode LTE Test2 Hmatrix 03

For another view, the four interpolated channel coefficient grids are plotted in constellation form. In this plot, the horizontal direction represents the 14 OFDM symbols in one subframe, and the vertical direction represents the four channel components: h11, h21, h12, and h22.

Each small plot shows the complex channel coefficient values for one OFDM symbol. The horizontal axis is the real part, and the vertical axis is the imaginary part. Since these are channel coefficients, the points do not represent user data modulation such as QPSK or QAM. Instead, they represent the complex channel response across subcarriers for each MIMO path.

The curved cluster shape means that the channel coefficient changes in both amplitude and phase across frequency. The shape is relatively compact because the plotted values are the interpolated channel coefficients, not the raw received data. The pattern also changes gradually from symbol to symbol, which indicates that the time-domain interpolation created a smooth channel estimate across the subframe.

This view is useful because it shows the phase behavior of each MIMO channel component more clearly than the magnitude-only grid image. From this plot, we can compare h11, h21, h12, and h22 and confirm that all four MIMO channel paths have been estimated over the full resource grid. These four interpolated channel grids will be used together as the 2x2 channel matrix for MIMO equalization.

IQ Decode LTE Test2 Hmatrix 04

Step 4: Equalization

This step performs MIMO equalization using the 2x2 H matrix. At this point, the four interpolated channel coefficient grids are already available as h11_interpolated, h21_interpolated, h12_interpolated, and h22_interpolated. These four grids form the channel matrix for every resource element. So the received data from antenna 0 and antenna 1 can now be equalized together.

Before applying equalization, it is useful to first plot the received resource grids as they are. resourceGrid_p0 is the received resource grid from RX antenna 0, and resourceGrid_p1 is the received resource grid from RX antenna 1. These two grids still include the MIMO channel effect, so the constellation is expected to be rotated, spread, and mixed. This is different from the SISO case because each received antenna may contain a mixture of signals from both transmit antenna ports.

The two ShowResourceGridIQ() calls display the before-equalization constellation for each receive antenna. The first plot shows rx_p0, which is the RX signal before equalization at antenna 0. The second plot shows rx_p1, which is the RX signal before equalization at antenna 1. The clip value is set to 25 to keep the plot range limited and make the constellation easier to inspect. These plots are used as the reference before applying the H-matrix based MIMO equalizer.

Following two plots show the received resource grid constellation before MIMO equalization. The first figure is rx_p0, which is the signal received at antenna 0. The second figure is rx_p1, which is the signal received at antenna 1. Each small plot represents one OFDM symbol from symbol 0 to symbol 13.

Before equalization, the constellation does not show clear QAM points. Instead, the samples form a circular or ring-shaped cloud around the center. This is because the received signal still includes the effect of the 2x2 MIMO channel. Each receive antenna observes a mixed version of the transmitted signals from both antenna ports, and the signal is also affected by phase rotation, amplitude variation, noise, and channel selectivity.

The two receive antennas show similar overall constellation shapes, but they are not identical. This is expected because rx_p0 and rx_p1 are captured from different receive antenna paths. Each antenna has its own channel response to TX port 0 and TX port 1. Therefore, the received constellation before equalization can have different spreading and rotation for each antenna.

These plots are used as the baseline before MIMO equalization. At this stage, the received data is still represented as y1 and y2. After applying the H matrix using h11, h21, h12, and h22, the receiver tries to separate the mixed signals and recover the transmitted streams x1 and x2. So the next plot after equalization should show whether the H-matrix based compensation makes the constellation cleaner and more structured.

IQ Decode LTE Test2 Eq 01

IQ Decode LTE Test2 Eq 02

Before doing MIMO equalization, the script prepares the variables used for H matrix processing, precoding, and Large CDD handling. First, the interpolated channel coefficient grids are renamed to shorter names. h11_interpolated becomes h11, h21_interpolated becomes h21, h12_interpolated becomes h12, and h22_interpolated becomes h22. This renaming is not mandatory. It is only done to make the following equations and code easier to read. The received resource grids are also renamed as rx_p0 and rx_p1, which represent the received signal at antenna 0 and antenna 1.

Then several copy grids are created from the original received resource grids. These grids are used to store intermediate results during the equalization process. In a real implementation, all of these intermediate resource grids may not be necessary. However, in this tutorial, they are useful because each processing step can be plotted separately. rx_p0_eq and rx_p1_eq are used for the final equalized result. rx_p0_W and rx_p1_W are used to show the result after handling the W matrix. rx_p0_D and rx_p1_D are used for the D matrix step. rx_p0_U and rx_p1_U are used for the U matrix step.

Next, the script defines the matrices used for LTE TM3 2x2 MIMO with Large CDD. The W matrix is the precoding normalization matrix. In this example, it has 1/sqrt(2) on the diagonal and zero on the off-diagonal elements. This means each transmit branch is scaled by 1/sqrt(2). The U matrix is also a fixed matrix. It combines the two layers with a fixed phase relationship. Since exp(-j 2π/2) is equal to -1, the U matrix becomes similar to a 2x2 Hadamard-like matrix with normalization by 1/sqrt(2).

The D matrix is different from W and U because it changes depending on the sequence index i. In the D(i) function, the first diagonal element is always 1, but the second diagonal element is exp(-j 2πi/2). This creates a phase change depending on i. This is the Large CDD part. It applies a cyclic-delay-related phase rotation to one branch, and the phase alternates depending on the data position.

So this preparation step creates all required inputs for the next equalization stage. The receiver already has the interpolated H matrix from CRS-based channel estimation. It also defines W, D, and U according to the LTE TM3 Large CDD structure. With these matrices ready, the script can compensate the received signal by considering not only the radio channel H, but also the transmit-side precoding and Large CDD operation.

The first equalization step is done by applying the inverse of the estimated H matrix to each received resource element. For every OFDM symbol and every subcarrier, the script creates a received signal vector RX using the two receive antenna values, rx_p0 and rx_p1. Then it creates the 2x2 channel matrix H using h11, h12, h21, and h22 at the same resource element position.

The main operation is simple.

This operation tries to remove the channel mixing effect from the received signal. In mathematical form, the received signal can be written as RX = H X. So if H is known, the transmitted signal estimate can be obtained as X = H^-1 RX. In the code, eq_re[0] is saved into rx_p0_eq, and eq_re[1] is saved into rx_p1_eq.

If the H matrix is singular at a certain resource element, the inverse cannot be calculated. In that case, the script prints the singular matrix location and uses the identity matrix as a fallback. This is only a simple protection for the tutorial code. In a real receiver, more robust methods such as MMSE equalization or pseudo-inverse based processing would usually be preferred.

After this operation, the result is plotted using ShowResourceGridIQ(). Now the points start to look more like a constellation. However, the result may still not look like the final 4QAM constellation that we expect. In this test, the original data is 4QAM, but most OFDM symbols still show unusual constellation shapes. This is because the signal has not yet been fully compensated for LTE TM3 precoding and Large CDD. The H-matrix equalization removes the radio channel effect, but the remaining signal still includes the transmit-side W, D, and U matrix effects. These will be handled in the following steps.

Following two plots show the result after applying the inverse H matrix to the received antenna signals. The first plot is Equalized rx_p0, and the second plot is Equalized rx_p1. Each small plot represents one OFDM symbol in the subframe.

Compared with the signal before equalization, the circular channel distortion is mostly removed. The points are no longer forming the same ring-shaped pattern. Instead, they are more concentrated around several regions, so the result starts to look more like a constellation. This means the 2x2 H matrix inversion is working and the receiver is partially separating the mixed signal observed at antenna 0 and antenna 1.

However, the constellation is still not the final expected 4QAM shape. The points are still spread and do not clearly collapse into four clean QAM clusters. This is expected in this step. The H matrix equalization compensates the radio channel between the transmit antenna ports and receive antennas, but the transmitted LTE TM3 signal still includes precoding and Large CDD processing. So the output after H inverse is not yet the original modulation symbol. It is still affected by the W, D, and U matrices.

Symbol 0 may look closer to a recognizable constellation because the effective precoding and reference signal structure can make that symbol appear simpler. But for the other OFDM symbols, the remaining precoding and Large CDD effect makes the constellation look unusual. Therefore, this result should be interpreted as an intermediate equalization result, not the final decoded 4QAM constellation. The next steps will compensate W, D, and U so that the constellation can be moved closer to the expected QAM points.

IQ Decode LTE Test2 Eq 03

IQ Decode LTE Test2 Eq 04

Step 5: Undo Precoding

This step starts the process of undoing the LTE MIMO precoding. In LTE spatial multiplexing, the original layer data is not mapped directly to antenna ports. The layer data is first processed by Large CDD using the U and D matrices, and then it is processed by the precoding matrix W. After that, the result is mapped to antenna port 0 and antenna port 1.

At the receiver side, the decoding process should be done in the reverse order. After the radio channel H is compensated, the next remaining effect is the transmit-side precoding. So the receiver first needs to undo W. After W is undone, the receiver then needs to undo the Large CDD effect from D and U.

The figure summarizes this transmit-side processing flow. The layer mapper creates layer 0 and layer 1 data. For LTE TM3 with Large CDD, the layer data is processed by U and D first, and then W is applied. The output becomes the signal transmitted from antenna port 0 and antenna port 1. This is why the constellation after H-matrix equalization still does not look like clean 4QAM. H equalization only removes the radio channel. It does not remove W, D, and U.

So from this point, the script performs additional reverse processing. First it compensates the W matrix. Then it compensates the D matrix, which changes depending on the sequence index. Finally it compensates the U matrix. After these steps, the constellation should become closer to the original layer-domain QAM symbols.

IQ Decode LTE Test2 Precoding 01

image source : sharetechnote

The result after undoing the precoding matrix W is shown in this step. The operation is simple. For each resource element, the script takes the two equalized values from rx_p0_eq and rx_p1_eq, makes them into a 2x1 vector, and multiplies the vector by the inverse of the W matrix. This removes the precoding matrix effect from the signal.

In the code, RX is the 2x1 vector created from the two equalized receive streams. Winv is the inverse of the W matrix. The result, eq_re, is calculated as inv(W) multiplied by RX. Then eq_re[0] is saved into rx_p0_W, and eq_re[1] is saved into rx_p1_W. This is done for every OFDM symbol and every subcarrier.

The code skips this operation for OFDM symbol 0 and simply copies the previous equalized values into rx_p0_W and rx_p1_W. This is done because symbol 0 is treated differently in this tutorial flow. For the other symbols, the inverse W operation is applied.

After this operation, the constellation may not look very different from the previous equalized constellation. This is expected. Even though many individual points move to new positions, the overall distribution can still fall into a similar set of clusters. So if all points are plotted with the same color, the visual difference may look small. However, the symbol positions are actually changed internally. If different symbols or resource elements are plotted with different colors, the redistribution before and after applying inv(W) can be seen more clearly.

So this step should be understood as an intermediate correction step. H matrix equalization removes the radio channel effect. inv(W) removes the precoding matrix effect. But the signal is still affected by Large CDD, so the constellation is not expected to fully collapse into the final 4QAM constellation yet. The remaining D and U matrix effects will be handled in the following steps.

Following two plots show the result after applying inv(W) to the H-equalized signal. The first plot is inv(W).rx_p0_eq, and the second plot is inv(W).rx_p1_eq. Each plot still shows 14 OFDM symbols separately.

Compared with the previous H-equalized result, the overall constellation shape does not change dramatically. This is expected because W in this example is a simple diagonal scaling matrix with 1/sqrt(2). Undoing W mainly changes the amplitude scaling of the two streams. It does not remove the Large CDD phase rotation. So the points may move internally, but the overall cluster pattern still looks similar when all points are plotted with the same color.

Now the constellation starts to show more structured groups, but it is still not the final 4QAM constellation. Many symbols show multiple cluster-like regions instead of four clean QAM points. This means that the radio channel H and the precoding matrix W have been compensated, but the Large CDD effect is still remaining. The remaining distortion is mainly from the D and U matrices, especially the D matrix because it changes with the data sequence index.

So this result should be treated as the output after undoing precoding only. It confirms that the inv(W) step has been applied, but it is not expected to fully recover the original layer-domain symbols yet. The next step is to undo the Large CDD operation by compensating the D matrix and then the U matrix.

IQ Decode LTE Test2 Precoding 02

IQ Decode LTE Test2 Precoding 03

Step 6: Undo Large CDD

This step undoes Large CDD. This is the last major correction step after H-matrix equalization and W-matrix compensation. In LTE TM3 spatial multiplexing, Large CDD is applied by using two matrices, U and D. The U matrix distributes the signal energy across the layers, and the D matrix applies a phase shift that depends on the sequence index i. After that, the W matrix is applied for precoding before the signal is mapped to antenna ports.

At the receiver side, the processing should be done in reverse order. The radio channel H is removed first by H inverse. Then the precoding matrix W is removed by inv(W). After that, Large CDD should be undone. Since the transmitter applies U first and then D during Large CDD processing, the receiver should undo D first and then undo U.

The important point is that U is a fixed matrix for the two-layer case, but D(i) changes depending on the data position index i. For two layers, D(i) is a diagonal matrix. The first layer keeps the same phase, while the second layer gets a phase rotation of exp(-j 2πi/2). This means the second layer alternates its phase depending on the sequence position. If this phase rotation is not compensated, the constellation can still look strange even after H and W compensation.

So this step prepares the final reverse processing. First, the script applies inv(D(i)) for each resource element using the proper index i. Then it applies inv(U). After these two operations, the Large CDD effect is removed, and the resulting constellation should become closer to the original 4QAM layer symbols.

IQ Decode LTE Test2 LargeCDD 01

image source : sharetechnote

Undoing the D matrix is the first part of undoing Large CDD. The input to this step is the result from the previous step, where the H matrix and W matrix effects have already been compensated. However, the signal still includes the phase rotation introduced by the D matrix. So for each resource element, the script multiplies the two-stream data vector by the inverse of D.

In the code, RX is created from rx_p0_W and rx_p1_W at the same OFDM symbol and subcarrier position. Then D(k) is generated using the subcarrier index k. This is important because the D matrix is not a fixed matrix like W. It changes with the sequence position. After D(k) is created, the script calculates inv(D(k)) and multiplies it with RX. The result is saved into rx_p0_D and rx_p1_D.

For OFDM symbol 0, the script does not apply the D inverse operation and simply copies the previous equalized values. For all other OFDM symbols, the inverse D operation is applied resource element by resource element.

This step mainly compensates the phase rotation caused by Large CDD. Because D is a diagonal matrix, it does not mix the two streams together. It mainly rotates the phase of one stream depending on k. After this operation, some constellation points may move significantly, but the final clean 4QAM constellation may still not appear yet. This is because only D has been undone. The U matrix effect is still remaining and will be removed in the next step.

Following two plots show the result after undoing the D matrix. The first plot is inv(D).inv(W).rx_p0_eq, and the second plot is inv(D).inv(W).rx_p1_eq. This means the signal has already gone through H-matrix equalization, W-matrix compensation, and now D-matrix compensation.

The overall constellation still looks similar to the previous inv(W) result. This can happen because the D matrix mainly applies a phase rotation that depends on the resource element index. Since the plotted points are all shown with the same color, the movement of individual points is not easy to recognize visually. Many points are moved to different positions, but the overall cloud can still look similar.

At this stage, the Large CDD effect is only partially removed. The D matrix phase rotation has been compensated, but the U matrix effect is still remaining. Because U mixes the two layers, the output after inv(D) still does not need to look like clean 4QAM yet. The constellation can still show multiple clusters and spreading.

So this result should be treated as another intermediate output. It confirms that the inverse D operation has been applied after inverse W. The final Large CDD recovery requires one more step, which is applying inv(U). After that, the signal should move closer to the original layer-domain QAM constellation.

IQ Decode LTE Test2 LargeCDD 02

IQ Decode LTE Test2 LargeCDD 03

As the last step, the script undoes the U matrix. This completes the reverse processing of Large CDD. The input to this step is rx_p0_D and rx_p1_D, where the H matrix, W matrix, and D matrix effects have already been compensated. The remaining operation is to remove the U matrix effect.

For each resource element, the script creates a 2x1 vector from rx_p0_D and rx_p1_D. Then it calculates the inverse of U and multiplies it with this vector. The result is saved into rx_p0_U and rx_p1_U.

The processing order is now complete.

This is the reverse order of the transmitter-side processing. At the transmitter, the layer data is processed by U, then D, then W, and then it passes through the radio channel H. At the receiver, the order must be reversed.

After applying inv(U), the constellation finally becomes close to the expected 4QAM shape. This means the receiver has not only compensated the radio channel, but also undone the LTE TM3 precoding and Large CDD processing. So rx_p0_U and rx_p1_U can be interpreted as the recovered layer-domain signals after MIMO equalization and reverse precoding processing.

 

Following two plots show the final result after undoing the U matrix. The first plot is inv(U).inv(D).inv(W).rx_p0_eq, and the second plot is inv(U).inv(D).inv(W).rx_p1_eq. This means the received signal has gone through the full reverse processing chain: H matrix equalization, inverse W, inverse D, and inverse U.

Now the constellation looks much closer to the expected 4QAM result. Instead of a circular cloud or many unclear clusters, most OFDM symbols show four main groups of points. These four groups correspond to the expected QPSK or 4QAM constellation. The points are still spread around each cluster because the signal still includes noise, residual channel estimation error, interpolation error, and possible residual synchronization or RF impairment. But the overall structure is now clearly recognizable as 4QAM.

This result also shows why the previous intermediate plots did not look like the expected constellation. H matrix equalization alone only removes the radio channel. In LTE TM3 with Large CDD, the received signal also includes the transmit-side processing from W, D, and U. Therefore, the receiver has to undo these matrices in reverse order. After inv(W), inv(D), and inv(U) are all applied, the layer-domain symbols can be recovered much more clearly.

So this plot is the final visual confirmation of the 2x2 MIMO decoding flow in this tutorial. The script first estimates the 2x2 channel matrix from CRS, then removes the channel effect using inv(H), then removes precoding using inv(W), then removes Large CDD using inv(D) and inv(U). After all these steps, the recovered signal becomes close to the original 4QAM constellation.

IQ Decode LTE Test2 LargeCDD 04

IQ Decode LTE Test2 LargeCDD 05

Script - Barebone

Following is the source script without any comments. I didn't put any comments on purpose for simplicity. For the original Python file, click here. For the sample IQ data that you may try with this code, you can download it from here (Unzip the file and change the file name and path in the source code).

NOTE : The code shown here is early and simplified version of the implementation. As more functions are added, the script may get longer and more complicated. For simplicity and readability, I will keep the early implementation as it is and just update the source code Python file.

import numpy as np

import matplotlib.pyplot as plt

import struct

import matplotlib.gridspec as gridspec

from scipy.signal import butter, lfilter

 

def generate_pss(n_id_2):

 

    u_shift = [25, 29, 34] # Root Zadoff-Chu sequence numbers

 

    d_u = []

 

    for n in range(62):  

        u = u_shift[n_id_2]

 

        if n <= 30:

            d = np.exp(-1j * np.pi * u * n * (n + 1) / 63)  

        else:

            d = np.exp(-1j * np.pi * u * (n + 1) * (n + 2) / 63)  

        

        d_u.append(d)

 

    return np.array(d_u)

 

 

def frequency_offset_estimation(received_pss, expected_pss, sample_rate = 30.72e6):

 

    phase_difference = np.angle(np.dot(received_pss, np.conj(expected_pss)))

    time_duration_pss = 62 / sample_rate

    frequency_offset = (phase_difference / (2 * np.pi)) / time_duration_pss

 

    return frequency_offset

 

 

def correct_frequency_offset(signal, frequency_offset, sample_rate, debug = False):

 

    if debug :

        print("[Debug Info] Frequency Offset ==> ",frequency_offset)

        print("[Debug Info] sample_rate ==> ",sample_rate)

 

    signal = np.array(signal, dtype=np.complex128)

    time = len(signal) / sample_rate

    correction = np.exp(-1j * 2 * np.pi * frequency_offset * time)

    corrected_signal = signal * correction

 

    # just for debugging

    if debug :

        original_signal_magnitude = np.abs(signal)

        corrected_signal_magnitude = np.abs(corrected_signal)

        # Check if the original signal magnitude and corrected signal magnitude are the same

        print("comparison of mag : ", np.allclose(original_signal_magnitude, corrected_signal_magnitude))

        print("correction : ", correction)

 

    return corrected_signal

 

def generate_sss(NID1,NID2):

 

    q_prime = np.floor(NID1/30)

    q = np.floor(((NID1+q_prime*(q_prime+1)/2))/30)

    m_prime = NID1 + q *(q+1)/2

    m0 = m_prime % 31

    m1 = (m0 + np.floor(m_prime/31) + 1) % 31

 

    # generate d_even() sequence

    x_s = np.zeros(31)

    x_s[:5] = [0, 0, 0, 0, 1]

    for i in range(26):

        x_s[i+5] = (x_s[i+2] + x_s[i]) % 2

 

    x_c = np.zeros(31)

    x_c[:5] = [0, 0, 0, 0, 1]

    for i in range(26):

        x_c[i+5] = (x_c[i+3] + x_c[i]) % 2

 

    s_tilda = 1 - 2*x_s

    c_tilda = 1 - 2*x_c

    s0_m0_even = np.zeros(31)

    s1_m1_even = np.zeros(31)

    for n in range(31):

        s0_m0_even[n] = s_tilda[int((n+m0)%31)]

        s1_m1_even[n] = s_tilda[int((n+m1)%31)]

 

    c0_even = np.zeros(31)

    for n in range(31):

        c0_even[n] = c_tilda[int((n+NID2)%31)]

 

    d_even_sub0 = s0_m0_even * c0_even

    d_even_sub5 = s1_m1_even * c0_even

 

    # generate d_odd() sequence

    x_z = np.zeros(31)

    x_z[:5] = [0, 0, 0, 0, 1]

    for i in range(26):

        x_z[i+5] = (x_z[i+4] + x_z[i+2] + x_z[i+1] + x_z[i]) % 2

 

    z_tilda = 1 - 2*x_z

    s1_m1_odd = np.zeros(31)

    s0_m0_odd = np.zeros(31)

    for n in range(31):

        s1_m1_odd[n] = s_tilda[int((n+m1)%31)]

        s0_m0_odd[n] = s_tilda[int((n+m0)%31)]

 

    c1_odd = np.zeros(31)

    for n in range(31):

        c1_odd[n] = c_tilda[int((n+NID2+3)%31)]

 

    z1_m0_odd = np.zeros(31)

    z1_m1_odd = np.zeros(31)

    for n in range(31):

        z1_m0_odd[n] = z_tilda[int((n+m0%8)%31)]

        z1_m1_odd[n] = z_tilda[int((n+m1%8)%31)]

 

    d_odd_sub0 = s1_m1_odd * c1_odd * z1_m0_odd

    d_odd_sub5 = s0_m0_odd * c1_odd * z1_m1_odd

 

    # calculate d_sub0

    d_sub0 = np.zeros(62)

    d_sub0[::2] = d_even_sub0

    d_sub0[1::2] = d_odd_sub0

    sss0 = d_sub0

 

    d_sub5 = np.zeros(62)

    d_sub5[::2] = d_even_sub5

    d_sub5[1::2] = d_odd_sub5

    sss5 = d_sub5

 

    return sss0, sss5

 

def cross_correlate(a, b):

    return np.sum(a * np.conj(b))

 

 

def DetectPSS(iq_samples, fft_bin_size=2048):

 

    pss_seqs = [generate_pss(i) for i in range(3)]

 

    max_correlation = 0

    max_window_samples = None

    max_correlation_index = None

    max_correlation_position = None

 

    max_correlation_pss0 = max_correlation_pss1 = max_correlation_pss2 = 0

 

    for j in range(len(iq_samples) - fft_bin_size):

        iq_window = iq_samples[j : j + fft_bin_size]

        frequency_samples = np.fft.fft(iq_window, fft_bin_size)

        frequency_samples = np.fft.fftshift(frequency_samples)

 

        for k in range(3):

            window_start = (fft_bin_size - 62) // 2

            window_end = window_start + 62

            window_samples = frequency_samples[window_start:window_end]

 

            correlation = np.sum(window_samples * np.conj(pss_seqs[k]))

 

            if k == 0 and abs(correlation) > max_correlation_pss0:

                max_correlation_pss0 = abs(correlation)

            elif k == 1 and abs(correlation) > max_correlation_pss1:

                max_correlation_pss1 = abs(correlation)

            elif k == 2 and abs(correlation) > max_correlation_pss2:

                max_correlation_pss2 = abs(correlation)

 

            if abs(correlation) > max_correlation:

                max_correlation = abs(correlation)

                max_correlation_index = k

                max_correlation_position = j

                max_window_samples = window_samples

                pss_start = j

                pss_time_domain = np.abs(iq_window)

 

    iq_window = iq_samples[max_correlation_position : max_correlation_position + fft_bin_size]

    frequency_samples = np.fft.fft(iq_window, fft_bin_size)

    frequency_samples = np.fft.fftshift(frequency_samples)

    received_pss = frequency_samples[993:1055]

    expected_pss = pss_seqs[max_correlation_index]

    frequency_offset = frequency_offset_estimation(received_pss, expected_pss)

    

    cp_lengths = [160, 144, 144, 144, 144, 144, 144, 160, 144, 144, 144, 144, 144, 144]

    pss_slot_index = 6

 

    first_symbol_start = (pss_start - cp_lengths[pss_slot_index]) - 5 * (cp_lengths[pss_slot_index]+fft_bin_size) - (cp_lengths[0]+fft_bin_size)

 

    if pss_start >= first_symbol_start + 6 * fft_bin_size and pss_start + 8 * fft_bin_size <= len(iq_samples):

        ofdm_symbols = []

        start_index = first_symbol_start

        for i in range(14):

            start_index = start_index + cp_lengths[i]

            end_index = start_index + fft_bin_size

            ofdm_symbols.append(iq_samples[start_index : end_index])

            start_index = end_index

        ofdm_symbols = np.array(ofdm_symbols)

    else:

        print("Not enough samples to extract 14 OFDM symbols.")

 

 

    return {

        'max_correlation': max_correlation,

        'max_correlation_pss0': max_correlation_pss0,

        'max_correlation_pss1': max_correlation_pss1,

        'max_correlation_pss2': max_correlation_pss2,

        'max_correlation_index': max_correlation_index,

        'max_correlation_position': max_correlation_position,

        'frequency_offset': frequency_offset,

        'max_window_samples': max_window_samples,

        'ofdm_symbols': ofdm_symbols,

        'pss_start' : pss_start,

        'pss_time_domain' : pss_time_domain ,      

    }

 

 

def getREs(ofdm_symbols, symNo, scList):

 

    try:

        return [ofdm_symbols[symNo][sc] for sc in scList]

    except IndexError:

        print("Invalid symbol number or subcarrier index")

        return []

 

 

def ShowPSS(iq_samples, detected_pss_freq_domain, max_correlation_index, ofdm_symbols, pss_start, pss_time_domain, fft_bin_size = 2048):

 

    import matplotlib.pyplot as plt

    import numpy as np

 

    fig = plt.figure(figsize=(15, 9))

 

    ax1 = plt.subplot2grid((7, 3), (0, 0), rowspan=2, colspan=1)

    ax3 = plt.subplot2grid((7, 3), (2, 0), rowspan=2, colspan=1)

    ax4 = plt.subplot2grid((7, 3), (4, 0), rowspan=2, colspan=1)

    ax5 = plt.subplot2grid((7, 3), (6, 0), rowspan=1, colspan=1)

 

    fft_iq_samples = np.fft.fft(iq_samples)

    fft_iq_samples = np.fft.fftshift(fft_iq_samples)

    fft_iq_samples_db = 20 * np.log10(np.abs(fft_iq_samples))

    ax1.plot(fft_iq_samples_db)

    ax1.set_title('Frequency Spectrum of IQ samples (dB scale)', fontsize=8)

    ax1.grid(True)

    ax1.set_xlabel('Frequency', fontsize=8)

    ax1.set_ylabel('Amplitude (dB)', fontsize=8)

    ax1.tick_params(axis='both', labelsize=8)

 

    fft_detected_pss = np.fft.fft(detected_pss_freq_domain)

    fft_detected_pss_normalized = fft_detected_pss / np.abs(fft_detected_pss).max()

 

    ax3.scatter(fft_detected_pss_normalized.real, fft_detected_pss_normalized.imag, label='Detected PSS')

    normalized_data = fft_detected_pss / np.abs(fft_detected_pss)

    ax3.scatter(normalized_data.real, normalized_data.imag, label='Normalized PSS', alpha=0.5)

    ax3.legend()

    ax3.grid(True)

    ax3.set_title('IQ Constellation of Detected PSS and Normalized PSS', fontsize=8)

    ax3.set_xlabel('In-Phase (I)', fontsize=8)

    ax3.set_ylabel('Quadrature (Q)', fontsize=8)

    ax3.legend(fontsize=8)

    ax3.tick_params(axis='both', labelsize=8)

 

    generated_pss = generate_pss(max_correlation_index)

    ax4.scatter(generated_pss.real, generated_pss.imag)

    ax4.grid(True)

    ax4.set_title('IQ Constellation of Generated PSS for NID2='+str(max_correlation_index), fontsize=8)

    ax4.set_xlabel('In-Phase (I)', fontsize=8)

    ax4.set_ylabel('Quadrature (Q)', fontsize=8)

    ax4.tick_params(axis='both', labelsize=8)

 

    mag_iq_samples = np.abs(iq_samples)

    ax5.plot(mag_iq_samples)

    ax5.plot(range(pss_start, pss_start + len(pss_time_domain)), pss_time_domain, label='Detected PSS', linestyle='--')

    ax5.grid(True)

    ax5.set_title('IQ in time doman and Detected PSS', fontsize=8)

    ax5.set_xlabel('Time', fontsize=8)

    ax5.set_ylabel('Amplitude', fontsize=8)

    ax5.legend(fontsize=8)

    ax5.tick_params(axis='both', labelsize=8)

 

    axs = [plt.subplot2grid((7, 3), (i % 7, i // 7 + 1)) for i in range(14)]

 

    for i in range(14):

        ax = axs[i]

        symbol = ofdm_symbols[i]

        freq_domain_symbol = np.fft.fft(symbol, fft_bin_size)

        freq_domain_symbol = np.fft.fftshift(freq_domain_symbol)

        freq_domain_symbol = freq_domain_symbol / np.max(np.abs(freq_domain_symbol))

        freq_domain_symbol_db = 20 * np.log10(np.abs(freq_domain_symbol))

        ax.plot(freq_domain_symbol_db)

        ax.set_ylim([-50, 0])

 

    plt.tight_layout()

    plt.show()

 

 

 

def ShowGrid(ofdm_symbols, nRB = 100, zmin=-25, zmax=0):

 

    frequency_domain_symbols = np.fft.fft(ofdm_symbols, axis=1)

    frequency_domain_symbols = np.fft.fftshift(frequency_domain_symbols, axes=1)

 

    magnitude = np.abs(frequency_domain_symbols)

    magnitude = magnitude / np.max(magnitude, axis=1, keepdims=True)

 

    magnitude_db = 20 * np.log10(magnitude + np.finfo(float).eps)

 

    total_subcarriers = frequency_domain_symbols.shape[1]

 

    start_index = (total_subcarriers - nRB * 12) // 2

    end_index = start_index + nRB * 12

 

    magnitude_db = magnitude_db[:, start_index:end_index]

 

    fig, ax = plt.subplots(figsize=(14, 5))

    

    cax = ax.imshow(magnitude_db, aspect='auto', cmap='RdBu_r', interpolation='none', vmin=zmin, vmax=zmax)

 

    fig.colorbar(cax)

    

    ax.set_xlabel("Subcarriers")

    ax.set_ylabel("OFDM Symbols")

    ax.set_title("Resource Grid")

 

    symbol_labels = ['Symbol {}'.format(i) for i in range(ofdm_symbols.shape[0])]

    ax.set_yticks(np.arange(ofdm_symbols.shape[0]))

    ax.set_yticklabels(symbol_labels)

 

    ax.set_xticks(np.arange(0, nRB * 12 + 1, 12*5))

    ax.set_xticklabels(np.arange(0, nRB * 12 + 1, 12*5))

 

    plt.tight_layout()

    plt.show()

 

 

def getResourceGrid(ofdm_symbols, nRB = 100):

    frequency_domain_symbols = np.fft.fft(ofdm_symbols, axis=1)

    frequency_domain_symbols = np.fft.fftshift(frequency_domain_symbols, axes=1)

 

    total_subcarriers = frequency_domain_symbols.shape[1]

 

    start_index = (total_subcarriers - nRB * 12) // 2

    end_index = start_index + nRB * 12

 

    resource_grid = frequency_domain_symbols[:, start_index:end_index]

 

    return resource_grid

 

 

def DetectSSS(resourceGrid, N_ID_2, frequency_offset, nRB = 100, samplingRate = 30.72e6, debug = False):

    if debug :

        print("[Debug Info] Frequency Offset ==> ",frequency_offset)

        print("[Debug Info] N_ID_2 ==> ",N_ID_2)

 

    noOfRE = nRB * 12

    sss_sc_indices = list(range(noOfRE-noOfRE//2-31, noOfRE-noOfRE//2+31))

 

    sss_no_correction = getREs(resourceGrid, 5, sss_sc_indices)

    sss = correct_frequency_offset(sss_no_correction, frequency_offset, samplingRate)

 

    max_corr = -np.inf

    best_N_ID_1 = None

    for N_ID_1 in range(168):

        expected_sss0, expected_sss5 = generate_sss(N_ID_1, N_ID_2)

        corr0 = np.abs(cross_correlate(sss, expected_sss0))

        corr5 = np.abs(cross_correlate(sss, expected_sss5))

        total_corr = max(corr0, corr5)

        if debug :

            print("corr0 = ", corr0, ",corr5 = ", corr5, ", N_ID_1 = ", N_ID_1)

 

        if total_corr > max_corr:

            max_corr = total_corr

            best_N_ID_1 = N_ID_1

            if debug :

                print("corr0 = ", corr0, ", corr5 = ", corr5, ", max_corr = ", max_corr, ", N_ID_1 = ", N_ID_1)

 

    if debug :

        print('Detected SSS (N_ID_1): ', best_N_ID_1)

 

    detected_subframe = 0 if cross_correlate(sss, generate_sss(best_N_ID_1, N_ID_2)[0]) > cross_correlate(sss, generate_sss(best_N_ID_1, N_ID_2)[1]) else 5

 

    return {

        'sss_no_correction': sss_no_correction,

        'sss_corrected': sss,

        'N_ID_1': best_N_ID_1,

        'detected_subframe': detected_subframe,

    }

 

 

def getPCI(N_ID_1, N_ID_2) :

    return 3*N_ID_1 + N_ID_2

 

 

def ShowSSS(sss_no_correction, sss_corrected, N_ID_1, N_ID_2,detected_subframe):

 

    sss_generated = np.array(generate_sss(N_ID_1, N_ID_2))

    if detected_subframe == 0 :

        sss_generated = sss_generated[0]

    else :

        sss_generated = sss_generated[1]    

 

    sss_no_correction = np.array(sss_no_correction, dtype=np.complex128)

    sss_corrected = np.array(sss_corrected, dtype=np.complex128)

    

    sss_no_correction = sss_no_correction / np.mean(np.abs(sss_no_correction))

    sss_corrected = sss_corrected / np.mean(np.abs(sss_corrected))

 

    data_list = [sss_no_correction, sss_corrected, sss_generated ]

    titles = ['sss_no_correction', 'sss_corrected', 'sss from N_ID_1']

    fig, axs = plt.subplots(3, 3, figsize=(9, 9))

 

    for i, data in enumerate(data_list):

        data = np.array(data, dtype=np.complex128)

 

        axs[0, i].scatter(data.real, data.imag, color='b')  

        axs[0, i].set_title(f'{titles[i]}', fontsize=8)

        axs[0, i].set_xlabel('Real', fontsize=8)

        axs[0, i].set_ylabel('Imaginary', fontsize=8)

        axs[0, i].tick_params(axis='both', which='major', labelsize=8)

        axs[0, i].set_xlim(-2,2)

        axs[0, i].set_ylim(-2,2)

 

        axs[1, i].plot(data.real)

        axs[1, i].set_title(f'{titles[i]}', fontsize=8)

        axs[1, i].set_xlabel('Data Index', fontsize=8)

        axs[1, i].set_ylabel('Real', fontsize=8)

        axs[1, i].tick_params(axis='both', which='major', labelsize=8)

        axs[1, i].set_ylim(-2,2)

 

        axs[2, i].plot(data.imag)

        axs[2, i].set_title(f'{titles[i]}', fontsize=8)

        axs[2, i].set_xlabel('Data Index', fontsize=8)

        axs[2, i].set_ylabel('Imaginary', fontsize=8)

        axs[2, i].tick_params(axis='both', which='major', labelsize=8)

        axs[2, i].set_ylim(-2,2)

 

    plt.tight_layout()

    plt.show()

 

 

 

def main():

    # Load IQ data from binary file

    with open('IQ/lte_20Mhz_rate30.72Mhz_dur_10ms_pci252.bin', 'rb') as f:  # expected PSS Seq = 0       

        iq_samples = np.fromfile(f, dtype=np.complex64)

 

    resultsPSS = DetectPSS(iq_samples)

 

    print(f'Detected PSS Sequence (N_ID_2): {resultsPSS["max_correlation_index"]}')

    print(f'Maximum correlation: {resultsPSS["max_correlation"]}')

    print(f'Maximum correlation for PSS0: {resultsPSS["max_correlation_pss0"]}')

    print(f'Maximum correlation for PSS1: {resultsPSS["max_correlation_pss1"]}')

    print(f'Maximum correlation for PSS2: {resultsPSS["max_correlation_pss2"]}')

    print(f'Position of maximum correlation: {resultsPSS["max_correlation_position"]}')

    print(f'Estimated Frequency Offset: {resultsPSS["frequency_offset"]} Hz')

 

    max_correlation_index = resultsPSS["max_correlation_index"]

    pss_start = resultsPSS["pss_start"]

    pss_time_domain = resultsPSS["pss_time_domain"]

    ofdm_symbols = resultsPSS["ofdm_symbols"]

    detected_pss_freq_domain = resultsPSS["max_window_samples"]

    frequency_offset = resultsPSS["frequency_offset"]

 

    resourceGrid = getResourceGrid(ofdm_symbols, nRB = 100)

    print("Dimension of resourceGrid = ",resourceGrid.shape)

 

    N_ID_2 = max_correlation_index

    resultsSSS = DetectSSS(resourceGrid, N_ID_2, frequency_offset, nRB = 100, samplingRate = 30.72e6)

 

    N_ID_1 = resultsSSS["N_ID_1"]

    sss_no_correction = resultsSSS["sss_no_correction"]

    sss_corrected = resultsSSS["sss_corrected"]

    detected_subframe = resultsSSS["detected_subframe"]

 

    print('Detected Subframe : ', detected_subframe)

 

    print('Detected N_ID_2 (PSS) : ', N_ID_2)

    print('Detected N_ID_1 (SSS) : ', N_ID_1)

 

    PCI = getPCI(N_ID_1, N_ID_2)

    print('Detected PCI: ', PCI)

 

    ShowPSS(iq_samples, detected_pss_freq_domain, max_correlation_index, ofdm_symbols, pss_start, pss_time_domain)

    ShowGrid(ofdm_symbols, nRB = 100, zmin = -25, zmax = 0)

    ShowSSS(sss_no_correction, sss_corrected, N_ID_1, N_ID_2,detected_subframe)  

 

 

if __name__ == "__main__":

    main()   

Script - Comments

For the readers who is not familiar with python script and packages that are used in the code. I put the comments for each lines.For the original Python file, click here. For the sample IQ data that you may try with this code, you can download it from here (Unzip the file and change the file name and path in the source code).

NOTE : The code shown here is early and simplified version of the implementation. As more functions are added, the script may get longer and more complicated. For simplicity and readability, I will keep the early implementation as it is and just update the source code Python file.

import numpy as np

import matplotlib.pyplot as plt

import struct

import matplotlib.gridspec as gridspec

from scipy.signal import butter, lfilter

 

 

# The function generate_pss(n_id_2) is designed to generate the Primary Synchronization Signal (PSS) for LTE systems. This generate 3 predefined PSS sequence specified by 3GPP.

# The three PSS signal will be used in the process of correlation check with the received signal

def generate_pss(n_id_2):

 

    u_shift = [25, 29, 34]

 

    d_u = []

 

    for n in range(62):  

        u = u_shift[n_id_2]

 

        if n <= 30:

            d = np.exp(-1j * np.pi * u * n * (n + 1) / 63)  

        else:

            d = np.exp(-1j * np.pi * u * (n + 1) * (n + 2) / 63)  

        

        d_u.append(d)

 

    return np.array(d_u)

 

 

# The function `frequency_offset_estimation` estimates the frequency offset in a received PSS.

# The frequency offset is the difference between the received signal's frequency and the transmitter's frequency.

 

# Parameters:

# `received_pss`: This is the Primary Synchronization Signal that was received.

# `expected_pss`: This is the expected or original Primary Synchronization Signal transmitted from the base station. This is used as a reference to estimate the frequency offset.

 

# The function assumes that both received_pss and expected_pss are numpy arrays of complex numbers representing the signal in the frequency domain.

 

# The function calculates the frequency offset by comparing the phase difference between the received_pss  and the expected_pss. This is done by taking the phase of the complex cross-correlation

# between the two signals.Then, this phase difference is used to estimate the frequency offset.

def frequency_offset_estimation(received_pss, expected_pss):

    phase_difference = np.angle(np.dot(received_pss, np.conj(expected_pss)))

    frequency_offset = phase_difference / (2 * np.pi * 62 / 30.72e6)  # Convert phase difference to frequency offset

    return frequency_offset

 

 

# The function `correct_frequency_offset` compensate the frequency offset for a given signal and returns the corrected signal.

 

# Parameters:

#  signal: The original signal to be corrected.

#   frequency_offset : The estimated frequency offset, in Hz, to be corrected.

#   sample_rate: The sample rate of the signal, in samples per second (Hz).

#   debug: If True, additional debug information is printed. Default is False.

 

def correct_frequency_offset(signal, frequency_offset, sample_rate, debug = False):

 

    if debug :

        print("[Debug Info] Frequency Offset ==> ",frequency_offset)

        print("[Debug Info] sample_rate ==> ",sample_rate)

 

    signal = np.array(signal, dtype=np.complex128)

    time = len(signal) / sample_rate

    correction = np.exp(-1j * 2 * np.pi * frequency_offset * time)

    corrected_signal = signal * correction

 

    # just for debugging

    if debug :

        original_signal_magnitude = np.abs(signal)

        corrected_signal_magnitude = np.abs(corrected_signal)

        # Check if the original signal magnitude and corrected signal magnitude are the same

        print("comparison of mag : ", np.allclose(original_signal_magnitude, corrected_signal_magnitude))

        print("correction : ", correction)

 

    return corrected_signal

 

 

# The function `generate_sss` generate the SSS signal for a given NID1 and NID2. This function returns two sequence. One for subframe 0 SSS and the other one for subframe 5 

# Parameters:

#   NID1 :  SSS Sequence Number.

#   NID2 : PSS Sequence Number.

 

def generate_sss(NID1,NID2):

 

    q_prime = np.floor(NID1/30)

    q = np.floor(((NID1+q_prime*(q_prime+1)/2))/30)

    m_prime = NID1 + q *(q+1)/2

    m0 = m_prime % 31

    m1 = (m0 + np.floor(m_prime/31) + 1) % 31

 

    # generate d_even() sequence

    x_s = np.zeros(31)

    x_s[:5] = [0, 0, 0, 0, 1]

    for i in range(26):

        x_s[i+5] = (x_s[i+2] + x_s[i]) % 2

 

    x_c = np.zeros(31)

    x_c[:5] = [0, 0, 0, 0, 1]

    for i in range(26):

        x_c[i+5] = (x_c[i+3] + x_c[i]) % 2

 

    s_tilda = 1 - 2*x_s

    c_tilda = 1 - 2*x_c

    s0_m0_even = np.zeros(31)

    s1_m1_even = np.zeros(31)

    for n in range(31):

        s0_m0_even[n] = s_tilda[int((n+m0)%31)]

        s1_m1_even[n] = s_tilda[int((n+m1)%31)]

 

    c0_even = np.zeros(31)

    for n in range(31):

        c0_even[n] = c_tilda[int((n+NID2)%31)]

 

    d_even_sub0 = s0_m0_even * c0_even

    d_even_sub5 = s1_m1_even * c0_even

 

    # generate d_odd() sequence

    x_z = np.zeros(31)

    x_z[:5] = [0, 0, 0, 0, 1]

    for i in range(26):

        x_z[i+5] = (x_z[i+4] + x_z[i+2] + x_z[i+1] + x_z[i]) % 2

 

    z_tilda = 1 - 2*x_z

    s1_m1_odd = np.zeros(31)

    s0_m0_odd = np.zeros(31)

    for n in range(31):

        s1_m1_odd[n] = s_tilda[int((n+m1)%31)]

        s0_m0_odd[n] = s_tilda[int((n+m0)%31)]

 

    c1_odd = np.zeros(31)

    for n in range(31):

        c1_odd[n] = c_tilda[int((n+NID2+3)%31)]

 

    z1_m0_odd = np.zeros(31)

    z1_m1_odd = np.zeros(31)

    for n in range(31):

        z1_m0_odd[n] = z_tilda[int((n+m0%8)%31)]

        z1_m1_odd[n] = z_tilda[int((n+m1%8)%31)]

 

    d_odd_sub0 = s1_m1_odd * c1_odd * z1_m0_odd

    d_odd_sub5 = s0_m0_odd * c1_odd * z1_m1_odd

 

    # calculate d_sub0

    d_sub0 = np.zeros(62)

    d_sub0[::2] = d_even_sub0

    d_sub0[1::2] = d_odd_sub0

    sss0 = d_sub0

 

   # calculate d_sub5

    d_sub5 = np.zeros(62)

    d_sub5[::2] = d_even_sub5

    d_sub5[1::2] = d_odd_sub5

    sss5 = d_sub5

 

    return sss0, sss5

 

 

# The function `cross_correlate` compute the cross correlation between the two sequence 

# Parameters:

#   a, b : sequence of complex numbers.

 

def cross_correlate(a, b):

    return np.sum(a * np.conj(b))

 

 

# The function `DetectPSS` is responsible for detecting the presence of the Primary Synchronization Signal (PSS) in a given set of IQ samples.

# Parameters:

#     `iq_samples`: These are the In-phase and Quadrature-phase samples of the received signal.

#      `fft_bin_size`: This is the size of the FFT bin used to transform the signal into the frequency domain. By default, the fft_bin_size is 2048, which is a common size in many LTE systems.

#                           However, it can be changed according to the system's specifications.

 

# The function works by first converting the time-domain signal (iq_samples) into the frequency domain using FFT.

# Then, it attempts to detect the PSS by cross-correlating the received frequency-domain signal with the expected PSS in the frequency domain.

# If the PSS is present, this cross-correlation should produce a strong peak in its output.

def DetectPSS(iq_samples, fft_bin_size=2048):

 

    pss_seqs = [generate_pss(i) for i in range(3)]

 

    max_correlation = 0

    max_window_samples = None

    max_correlation_index = None

    max_correlation_position = None

 

    max_correlation_pss0 = max_correlation_pss1 = max_correlation_pss2 = 0

 

    for j in range(len(iq_samples) - fft_bin_size):

        iq_window = iq_samples[j : j + fft_bin_size]

        frequency_samples = np.fft.fft(iq_window, fft_bin_size)

        frequency_samples = np.fft.fftshift(frequency_samples)

 

        for k in range(3):

            # Calculate the start and end of the window

            window_start = (fft_bin_size - 62) // 2

            window_end = window_start + 62

            window_samples = frequency_samples[window_start:window_end]

 

            correlation = np.sum(window_samples * np.conj(pss_seqs[k]))

 

            if k == 0 and abs(correlation) > max_correlation_pss0:

                max_correlation_pss0 = abs(correlation)

            elif k == 1 and abs(correlation) > max_correlation_pss1:

                max_correlation_pss1 = abs(correlation)

            elif k == 2 and abs(correlation) > max_correlation_pss2:

                max_correlation_pss2 = abs(correlation)

 

            if abs(correlation) > max_correlation:

                max_correlation = abs(correlation)

                max_correlation_index = k

                max_correlation_position = j

                max_window_samples = window_samples

                pss_start = j

                pss_time_domain = np.abs(iq_window)

 

    # Calculate Frequency Offset

    iq_window = iq_samples[max_correlation_position : max_correlation_position + fft_bin_size]

    frequency_samples = np.fft.fft(iq_window, fft_bin_size)

    frequency_samples = np.fft.fftshift(frequency_samples)

    received_pss = frequency_samples[993:1055]

    expected_pss = pss_seqs[max_correlation_index]

    frequency_offset = frequency_offset_estimation(received_pss, expected_pss)

    

    # Specify the cyclic prefix lengths for each OFDM symbol in a slot. This is CP length for 20Mhz BW. CP removal for other BW will be implemented later release

    cp_lengths = [160, 144, 144, 144, 144, 144, 144, 160, 144, 144, 144, 144, 144, 144]

    pss_slot_index = 6

 

    # Calculate the start position of the first OFDM symbol by adding up the CP lengths

    first_symbol_start = (pss_start - cp_lengths[pss_slot_index]) - 5 * (cp_lengths[pss_slot_index]+fft_bin_size) - (cp_lengths[0]+fft_bin_size)

 

    # Check if there are enough samples before and after PSS

    if pss_start >= first_symbol_start + 6 * fft_bin_size and pss_start + 8 * fft_bin_size <= len(iq_samples):

        # Extract 14 OFDM symbols where PSS is located, taking into account the CP

        ofdm_symbols = []

        start_index = first_symbol_start

        for i in range(14):

            start_index = start_index + cp_lengths[i]

            end_index = start_index + fft_bin_size

            ofdm_symbols.append(iq_samples[start_index : end_index])

            start_index = end_index

        ofdm_symbols = np.array(ofdm_symbols)

    else:

        print("Not enough samples to extract 14 OFDM symbols.")

 

 

# This function returns various values in the form of a dictionary.

    return {

        'max_correlation': max_correlation,   # max_correlation value detected from the process

        'max_correlation_pss0': max_correlation_pss0,   # max_correlation value for pss0 detected from the process (this is just for the reference)

        'max_correlation_pss1': max_correlation_pss1,  # max_correlation value for pss1 detected from the process  (this is just for the reference)

        'max_correlation_pss2': max_correlation_pss2,  # max_correlation value for pss2 detected from the process  (this is just for the reference)

        'max_correlation_index': max_correlation_index,# the pss index (i.e, pss sequence number that gives max_correlation). this is the detected PSS sequence

       'max_correlation_position': max_correlation_position,# this is position of the sliding window that gives the max correlation. It is the start of the sliding window on the received I/Q data

        'frequency_offset': frequency_offset,

        'max_window_samples': max_window_samples,# this is the frequency domain I/Q of the detected PSS

        'ofdm_symbols': ofdm_symbols,# 2D array of one radio frame OFDM symbols (14 symbols). Frequency domain I/Q data after CP removal

        'pss_start' : pss_start,# this is same as max_correlation_position

        'pss_time_domain' : pss_time_domain ,  # this is the sliding window I/Q data that carries the detected PSS    

    }

 

# The function `getREs` extracts the value from the given resource grid pointed by OFDM symbol number and the list of subcarriers.

# Parameters:

#   `resourceGrid`: OFDM resource Grid.

#   `symNo`: OFDM symbol number within the resource grid. This starts from 0.

#   'scList' : the list of subcarrier index (the index starting from 0) .

def getREs(resourceGrid, symNo, scList):

 

    try:

        # Retrieve the REs corresponding to the given symbol number and subcarrier indices

        return [ofdm_symbols[symNo][sc] for sc in scList]

    except IndexError:

        print("Invalid symbol number or subcarrier index")

        return []

 

 

# The function is used to visualize the detection of the Primary Synchronization Signal (PSS) in a given set of IQ samples. The function displays several plots to visualize the PSS detection process

 

# Parameters:

#     `iq_samples`: These are the In-phase and Quadrature-phase samples of the received signal.

#      `detected_pss_freq_domain`: This is the detected PSS in the frequency domain, obtained from the `DetectPSS` function.

#      `max_correlation_index`: This is the index in the IQ samples where the correlation between the expected and received PSS is the maximum.

#                                            This indicates the location of the PSS in the received signal.

#      `ofdm_symbols`: These are the OFDM symbols extracted from the received signal.

#      `pss_start`: This is the starting index of the PSS in the IQ samples, obtained from the `DetectPSS` function.

#      `pss_time_domain`: this is the time-domain IQ sequence in which the PSS is located.

#      `fft_bin_size`: This is the size of the FFT bin used to transform the signal into the frequency domain. this is same as the length of the sliding window searching for PSS.

#      By default, the fft_bin_size is 2048, which is a common size in many LTE systems. However, it can be changed according to the system's specifications.

def ShowPSS(iq_samples, detected_pss_freq_domain, max_correlation_index, ofdm_symbols, pss_start, pss_time_domain, fft_bin_size = 2048):

 

    # Create a single figure

    fig = plt.figure(figsize=(15, 9))

 

    # Create the left side subplots

    ax1 = plt.subplot2grid((7, 3), (0, 0), rowspan=2, colspan=1)

    ax3 = plt.subplot2grid((7, 3), (2, 0), rowspan=2, colspan=1)

    ax4 = plt.subplot2grid((7, 3), (4, 0), rowspan=2, colspan=1)

    ax5 = plt.subplot2grid((7, 3), (6, 0), rowspan=1, colspan=1)

 

    # Frequency Spectrum of IQ samples

    fft_iq_samples = np.fft.fft(iq_samples)

    fft_iq_samples = np.fft.fftshift(fft_iq_samples)

    fft_iq_samples_db = 20 * np.log10(np.abs(fft_iq_samples))

    ax1.plot(fft_iq_samples_db)

    ax1.set_title('Frequency Spectrum of IQ samples (dB scale)', fontsize=8)

    ax1.grid(True)

    ax1.set_xlabel('Frequency', fontsize=8)

    ax1.set_ylabel('Amplitude (dB)', fontsize=8)

    ax1.tick_params(axis='both', labelsize=8)

 

    # Frequency Spectrum of Detected PSS

    fft_detected_pss = np.fft.fft(detected_pss_freq_domain)

    # Normalize by dividing all elements by the maximum magnitude

    fft_detected_pss_normalized = fft_detected_pss / np.abs(fft_detected_pss).max()

 

    # Plot the Detected PSS

    ax3.scatter(fft_detected_pss_normalized.real, fft_detected_pss_normalized.imag, label='Detected PSS')

    # Normalized PSS data

    normalized_data = fft_detected_pss / np.abs(fft_detected_pss)

    # Plot the Normalized data

    ax3.scatter(normalized_data.real, normalized_data.imag, label='Normalized PSS', alpha=0.5)

    ax3.legend()

    ax3.grid(True)

    ax3.set_title('IQ Constellation of Detected PSS and Normalized PSS', fontsize=8)

    ax3.set_xlabel('In-Phase (I)', fontsize=8)

    ax3.set_ylabel('Quadrature (Q)', fontsize=8)

    ax3.legend(fontsize=8)

    ax3.tick_params(axis='both', labelsize=8)

 

    # Generate the PSS for the detected PCI

    generated_pss = generate_pss(max_correlation_index)

    # IQ plot of the generated PSS

    ax4.scatter(generated_pss.real, generated_pss.imag)

    ax4.grid(True)

    ax4.set_title('IQ Constellation of Generated PSS for NID2='+str(max_correlation_index), fontsize=8)

    ax4.set_xlabel('In-Phase (I)', fontsize=8)

    ax4.set_ylabel('Quadrature (Q)', fontsize=8)

    ax4.tick_params(axis='both', labelsize=8)

 

    # Frequency Spectrum of IQ samples

    mag_iq_samples = np.abs(iq_samples)

    ax5.plot(mag_iq_samples)

    ax5.plot(range(pss_start, pss_start + len(pss_time_domain)), pss_time_domain, label='Detected PSS', linestyle='--')

    ax5.grid(True)

    ax5.set_title('IQ in time doman and Detected PSS', fontsize=8)

    ax5.set_xlabel('Time', fontsize=8)

    ax5.set_ylabel('Amplitude', fontsize=8)

    ax5.legend(fontsize=8)

    ax5.tick_params(axis='both', labelsize=8)

 

    # Create the right side subplots

    axs = [plt.subplot2grid((7, 3), (i % 7, i // 7 + 1)) for i in range(14)]

 

    # Code to generate the plots for axs

    for i in range(14):

        ax = axs[i]

        symbol = ofdm_symbols[i]

        freq_domain_symbol = np.fft.fft(symbol, fft_bin_size)

        freq_domain_symbol = np.fft.fftshift(freq_domain_symbol)

        freq_domain_symbol = freq_domain_symbol / np.max(np.abs(freq_domain_symbol))

        freq_domain_symbol_db = 20 * np.log10(np.abs(freq_domain_symbol))

        ax.plot(freq_domain_symbol_db)

        ax.set_ylim([-50, 0])

 

    plt.tight_layout()

    plt.show()

 

 

# The function is used to visualize the OFDM grid, representing the Resource Blocks (RBs) of an LTE signal in Resource Grid. This grid is particularly useful for visualizing  the allocation of resources

#  in the frequency-time space in a LTE system.

# The function generates a heatmap style plot of the OFDM grid. The x-axis represents the time in terms of  the OFDM symbol number, and the y-axis represents the frequency in terms of

#  the sub-carrier number. The color of each cell in the grid corresponds to the amplitude of the corresponding sub-carrier at the corresponding time, with warmer colors indicating higher amplitudes.

 

# Parameters:

#     `ofdm_symbols`: These are OFDM symbols extracted from the received signal. These are 2D array comprised of OFDM symbols after removing CP (Cyclic Prefix).

#      `nRB`: This parameter represents the number of Resource Blocks in the LTE signal.  The default value is 100,  which corresponds to a 20 MHz bandwidth in LTE, but it can be adjusted

#                according to the system's specifications.

 

#      `zmin`, `zmax`: These parameters define the color range for the plot. Any amplitude that falls within  this range is represented with a color from the colormap.

#                            `zmin` is the minimum dB value and `zmax`  is the maximum dB value for the color range. By default, `zmin` is -25 dB and `zmax` is 0 dB.

 

def ShowGrid(ofdm_symbols, nRB = 100, zmin=-25, zmax=0):

 

    # Calculate the magnitude of the OFDM symbols

    frequency_domain_symbols = np.fft.fft(ofdm_symbols, axis=1)

    frequency_domain_symbols = np.fft.fftshift(frequency_domain_symbols, axes=1)

 

    # Compute the magnitude of the frequency domain symbols

    magnitude = np.abs(frequency_domain_symbols)

    magnitude = magnitude / np.max(magnitude, axis=1, keepdims=True)

 

    # Convert to decibels and normalize

    magnitude_db = 20 * np.log10(magnitude + np.finfo(float).eps)

 

    # Number of subcarriers

    total_subcarriers = frequency_domain_symbols.shape[1]

 

    # Calculate the start and end index for the desired number of RBs

    start_index = (total_subcarriers - nRB * 12) // 2

    end_index = start_index + nRB * 12

 

    # Slice the magnitude_db array to include only the desired RBs

    magnitude_db = magnitude_db[:, start_index:end_index]

 

    # Create a figure and an axes object

    fig, ax = plt.subplots(figsize=(14, 5))

    

    # Use the imshow function to create the 2D graph

    cax = ax.imshow(magnitude_db, aspect='auto', cmap='RdBu_r', interpolation='none', vmin=zmin, vmax=zmax)

 

    # Create a color bar to show the mapping between values and colors

    fig.colorbar(cax)

    

    # Add labels and title

    ax.set_xlabel("Subcarriers")

    ax.set_ylabel("OFDM Symbols")

    ax.set_title("Resource Grid")

 

    # Set y-axis labels for every symbol

    symbol_labels = ['Symbol {}'.format(i) for i in range(ofdm_symbols.shape[0])]

    ax.set_yticks(np.arange(ofdm_symbols.shape[0]))

    ax.set_yticklabels(symbol_labels)

 

    # Adjust x-axis ticks and labels to start from 0 and go up to nRB*12 with step size of 12

    ax.set_xticks(np.arange(0, nRB * 12 + 1, 12*5))

    ax.set_xticklabels(np.arange(0, nRB * 12 + 1, 12*5))

 

    # Display the graph

    plt.tight_layout()

    plt.show()

 

 

# The function constructs 2D resource grid (ofdm symbol index vs subcarrier) from the 2D time domain array (ofdm symbol index vs ofdm symbol in time domain)

 

# Parameters:

#     `ofdm_symbols`: These are OFDM symbols extracted from the received signal. These are 2D array comprised of OFDM symbols after removing CP (Cyclic Prefix).

#      `nRB`: This parameter represents the number of Resource Blocks in the LTE signal.  The default value is 100,  which corresponds to a 20 MHz bandwidth in LTE, but it can be adjusted

 

def getResourceGrid(ofdm_symbols, nRB = 100):

    # Calculate the magnitude of the OFDM symbols

    frequency_domain_symbols = np.fft.fft(ofdm_symbols, axis=1)

    frequency_domain_symbols = np.fft.fftshift(frequency_domain_symbols, axes=1)

 

    # Number of subcarriers

    total_subcarriers = frequency_domain_symbols.shape[1]

 

    # Calculate the start and end index for the desired number of RBs

    start_index = (total_subcarriers - nRB * 12) // 2

    end_index = start_index + nRB * 12

 

    # Slice the frequency_domain_symbols array to include only the desired RBs

    resource_grid = frequency_domain_symbols[:, start_index:end_index]

 

    return resource_grid

 

 

 

# The function detects SSS embedded in a given resourceGrid

 

# Parameters:

#    resourceGrid : The resource grid containing the received signal.

#    N_ID_2: The second cell identity index.

#    frequency_offset: The estimated frequency offset to be corrected.

#    nRB : The number of Resource Blocks. Default is 100.

#    samplingRate: The sampling rate of the signal, in samples per second (Hz). Default is 30.72e6.

#    debug (bool, optional): If True, additional debug information is printed. Default is False.

 

# This function returns various values in the form of dictionary as follows :

#    'sss_no_correction': The extracted SSS without frequency offset correction.

#    'sss_corrected': The extracted SSS with frequency offset correction.

#    'N_ID_1': The detected first cell identity index. (The detected SSS Sequence Number)

#    'detected_subframe': The detected subframe number (0 or 5).

 

def DetectSSS(resourceGrid, N_ID_2, frequency_offset, nRB = 100, samplingRate = 30.72e6, debug = False):

    if debug :

        print("[Debug Info] Frequency Offset ==> ",frequency_offset)

        print("[Debug Info] N_ID_2 ==> ",N_ID_2)

 

    # Generate the list of subcarrier indices for SSS

    noOfRE = nRB * 12

    sss_sc_indices = list(range(noOfRE-noOfRE//2-31, noOfRE-noOfRE//2+31))

 

    # Extract SSS for the first and second slots

    sss_no_correction = getREs(resourceGrid, 5, sss_sc_indices)

 

    # Apply frequency offset correction to the extracted SSS

    sss = correct_frequency_offset(sss_no_correction, frequency_offset, samplingRate)

 

    # Calculate correlations for all possible N_ID_1

    max_corr = -np.inf

    best_N_ID_1 = None

    for N_ID_1 in range(168):

        # Generate expected SSS

        expected_sss0, expected_sss5 = generate_sss(N_ID_1, N_ID_2)

        # Calculate correlation

        corr0 = np.abs(cross_correlate(sss, expected_sss0))

        corr5 = np.abs(cross_correlate(sss, expected_sss5))

        total_corr = max(corr0, corr5)

        if debug :

            print("corr0 = ", corr0, ",corr5 = ", corr5, ", N_ID_1 = ", N_ID_1)

 

        # If this correlation is the highest we've seen, save it and the current N_ID_1

        if total_corr > max_corr:

            max_corr = total_corr

            best_N_ID_1 = N_ID_1

            if debug :

                print("corr0 = ", corr0, ", corr5 = ", corr5, ", max_corr = ", max_corr, ", N_ID_1 = ", N_ID_1)

 

    if debug :

        print('Detected SSS (N_ID_1): ', best_N_ID_1)

 

    # Detect the subframe (0 or 5) based on the correlation results

    detected_subframe = 0 if cross_correlate(sss, generate_sss(best_N_ID_1, N_ID_2)[0]) > cross_correlate(sss, generate_sss(best_N_ID_1, N_ID_2)[1]) else 5

 

    return {

        'sss_no_correction': sss_no_correction,

        'sss_corrected': sss,

        'N_ID_1': best_N_ID_1,

        'detected_subframe': detected_subframe,

    }

 

 

# The function calculate PCI from N_ID_1 (SSS sequence number) and N_ID_2 (PSS Sequence Number)

def getPCI(N_ID_1, N_ID_2) :

    return 3*N_ID_1 + N_ID_2

 

 

# This function is to visualize the Secondary Synchronization Signal (SSS) before and after frequency offset correction, and the generated SSS for a given N_ID_1 and N_ID_2.

 

# Parameters:

#         sss_no_correction : The SSS extracted from the resource grid before frequency offset correction.

#         sss_corrected : The SSS after frequency offset correction.

#         N_ID_1:SSS Sequence Number.

#         N_ID_2: PSS Sequence Number .

#         detected_subframe : The detected subframe number (0 or 5).

 

def ShowSSS(sss_no_correction, sss_corrected, N_ID_1, N_ID_2,detected_subframe):

 

    sss_generated = np.array(generate_sss(N_ID_1, N_ID_2))

    if detected_subframe == 0 :

        sss_generated = sss_generated[0]

    else :

        sss_generated = sss_generated[1]    

 

    sss_no_correction = np.array(sss_no_correction, dtype=np.complex128)

    sss_corrected = np.array(sss_corrected, dtype=np.complex128)

    

    # Normalize sss_no_correction and sss_corrected

    sss_no_correction = sss_no_correction / np.mean(np.abs(sss_no_correction))

    sss_corrected = sss_corrected / np.mean(np.abs(sss_corrected))

 

    data_list = [sss_no_correction, sss_corrected, sss_generated ]

    #print("[Debug Info] ==> " , data_list)

    titles = ['sss_no_correction', 'sss_corrected', 'sss from N_ID_1']

    fig, axs = plt.subplots(3, 3, figsize=(9, 9))

 

    for i, data in enumerate(data_list):

        # Ensure data is a numpy array with complex numbers

        data = np.array(data, dtype=np.complex128)

 

        axs[0, i].scatter(data.real, data.imag, color='b')  

        axs[0, i].set_title(f'{titles[i]}', fontsize=8)

        axs[0, i].set_xlabel('Real', fontsize=8)

        axs[0, i].set_ylabel('Imaginary', fontsize=8)

        axs[0, i].tick_params(axis='both', which='major', labelsize=8)

        axs[0, i].set_xlim(-2,2)

        axs[0, i].set_ylim(-2,2)

 

        axs[1, i].plot(data.real)

        axs[1, i].set_title(f'{titles[i]}', fontsize=8)

        axs[1, i].set_xlabel('Data Index', fontsize=8)

        axs[1, i].set_ylabel('Real', fontsize=8)

        axs[1, i].tick_params(axis='both', which='major', labelsize=8)

        axs[1, i].set_ylim(-2,2)

 

        axs[2, i].plot(data.imag)

        axs[2, i].set_title(f'{titles[i]}', fontsize=8)

        axs[2, i].set_xlabel('Data Index', fontsize=8)

        axs[2, i].set_ylabel('Imaginary', fontsize=8)

        axs[2, i].tick_params(axis='both', which='major', labelsize=8)

        axs[2, i].set_ylim(-2,2)

 

    plt.tight_layout()

    plt.show()

 

 

 

# Test Code 

 

def main():

    # Load IQ data from binary file 

    with open('IQ/lte_20Mhz_rate30.72Mhz_dur_10ms_pci252.bin', 'rb') as f:  # expected PSS Seq = 0s        

        iq_samples = np.fromfile(f, dtype=np.complex64)

 

    # Detect PSS 

    results = DetectPSS(iq_samples)

 

   # Print various outputs returned by DetectPSS function 

    print(f'Detected PSS Sequence (N_ID_2): {results["max_correlation_index"]}')

    print(f'Maximum correlation: {results["max_correlation"]}')

    print(f'Maximum correlation for PSS0: {results["max_correlation_pss0"]}')

    print(f'Maximum correlation for PSS1: {results["max_correlation_pss1"]}')

    print(f'Maximum correlation for PSS2: {results["max_correlation_pss2"]}')

    print(f'Position of maximum correlation: {results["max_correlation_position"]}')

    print(f'Estimated Frequency Offset: {results["frequency_offset"]} Hz')

 

   # Extract the values from DetectPSS return that are necessary for Visualization 

    max_correlation_index = results["max_correlation_index"]

    pss_start = results["pss_start"]

    pss_time_domain = results["pss_time_domain"]

    ofdm_symbols = results["ofdm_symbols"]

    detected_pss_freq_domain = results["max_window_samples"]

 

    resourceGrid = getResourceGrid(ofdm_symbols, nRB = 100)

    print("Dimension of resourceGrid = ",resourceGrid.shape)

 

    # Detect SSS

    N_ID_2 = max_correlation_index

    resultsSSS = DetectSSS(resourceGrid, N_ID_2, frequency_offset, nRB = 100, samplingRate = 30.72e6)

 

    N_ID_1 = resultsSSS["N_ID_1"]

    sss_no_correction = resultsSSS["sss_no_correction"]

    sss_corrected = resultsSSS["sss_corrected"]

    detected_subframe = resultsSSS["detected_subframe"]

 

    print('Detected Subframe : ', detected_subframe)

 

    print('Detected N_ID_2 (PSS) : ', N_ID_2)

    print('Detected N_ID_1 (SSS) : ', N_ID_1)

 

    PCI = getPCI(N_ID_1, N_ID_2)

    print('Detected PCI: ', PCI)

 

    ShowPSS(iq_samples, detected_pss_freq_domain, max_correlation_index, ofdm_symbols, pss_start, pss_time_domain)

    ShowGrid(ofdm_symbols, nRB = 100, zmin = -25, zmax = 0)

    ShowSSS(sss_no_correction, sss_corrected, N_ID_1, N_ID_2,detected_subframe)  

 

 

if __name__ == "__main__":

    main()

© 2026 All rights reserved | Legal notice |