diff --git a/CMakeLists.txt b/CMakeLists.txt index 58d95f281..9c166e60e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,7 @@ PROJECT (LIBLTE) LIST(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules") INCLUDE(libLTEPackage) #setup cpack +INCLUDE(BuildMex) include(CTest) set( CTEST_MEMORYCHECK_COMMAND valgrind ) @@ -140,7 +141,6 @@ MACRO (APPEND_INTERNAL_LIST LIST_NAME VALUE) ENDIF (${LIST_NAME}) ENDMACRO (APPEND_INTERNAL_LIST) - ######################################################################## # Print summary ######################################################################## diff --git a/cmake/modules/BuildMex.cmake b/cmake/modules/BuildMex.cmake new file mode 100644 index 000000000..124e149a8 --- /dev/null +++ b/cmake/modules/BuildMex.cmake @@ -0,0 +1,65 @@ +# BuildMex.cmake +# Author: Kent Williams norman-k-williams at uiowa.edu +# Modified by Ismael Gomez, 2014 + +include(CMakeParseArguments) + +if(NOT MATLAB_FOUND) + find_package(MATLAB) +endif() + +if(NOT OCTAVE_FOUND) + find_package(OCTAVE) +endif() + +# CMake 2.8.12 & earlier apparently don't define the +# Mex script path, so find it. +if(NOT MATLAB_MEX_PATH) + find_program( MATLAB_MEX_PATH mex + HINTS ${MATLAB_ROOT}/bin + PATHS ${MATLAB_ROOT}/bin + DOC "The mex program path" + ) +endif() + +IF (MATLAB_FOUND) + message(STATUS "Found MATLAB in ${MATLAB_ROOT}") +ELSE(MATLAB_FOUND) + message(STATUS "Could NOT find MATLAB. MEX files won't be compiled") +ENDIF(MATLAB_FOUND) + +# +# BuildMex -- arguments +# MEXNAME = root of mex library name +# SOURCE = list of source files +# LIBRARIES = libraries needed to link mex library +FUNCTION(BuildMex) + set(oneValueArgs MEXNAME) + set(multiValueArgs SOURCES LIBRARIES) + cmake_parse_arguments(BuildMex "" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) + if (MATLAB_FOUND) + add_library(${BuildMex_MEXNAME}-mat SHARED ${BuildMex_SOURCES}) + target_include_directories(${BuildMex_MEXNAME}-mat PUBLIC ${MATLAB_INCLUDE_DIR}) + set_target_properties(${BuildMex_MEXNAME}-mat PROPERTIES + SUFFIX "${MATLAB_MEX_EXTENSION}" + PREFIX "liblte_" + OUTPUT_NAME "${BuildMex_MEXNAME}" + COMPILE_FLAGS "-fvisibility=default ${MATLAB_MEX_CFLAGS}" + ) + target_link_libraries(${BuildMex_MEXNAME}-mat ${BuildMex_LIBRARIES} ${MATLAB_MEX_LIBRARY}) + install(TARGETS ${BuildMex_MEXNAME}-mat DESTINATION mex) + endif(MATLAB_FOUND) + if (OCTAVE_FOUND) + add_library(${BuildMex_MEXNAME}-oct SHARED ${BuildMex_SOURCES}) + target_include_directories(${BuildMex_MEXNAME}-oct PUBLIC ${OCTAVE_INCLUDE_DIR}) + set_target_properties(${BuildMex_MEXNAME}-oct PROPERTIES + SUFFIX ".${OCTAVE_MEXFILE_EXT}" + PREFIX "liblte_" + OUTPUT_NAME "${BuildMex_MEXNAME}" + COMPILE_FLAGS "-fvisibility=default ${OCTAVE_MEX_CFLAGS} -DUNDEF_BOOL" + ) + target_link_libraries(${BuildMex_MEXNAME}-oct ${BuildMex_LIBRARIES} ${OCTAVE_LIBRARIES}) + install(TARGETS ${BuildMex_MEXNAME}-oct DESTINATION mex) + endif (OCTAVE_FOUND) +ENDFUNCTION(BuildMex) + diff --git a/cmake/modules/CheckFunctionExists.c b/cmake/modules/CheckFunctionExists.c index 1fbdf9c53..314263dbf 100644 --- a/cmake/modules/CheckFunctionExists.c +++ b/cmake/modules/CheckFunctionExists.c @@ -1,13 +1,13 @@ #ifdef CHECK_FUNCTION_EXISTS -uint8_t CHECK_FUNCTION_EXISTS(); +char CHECK_FUNCTION_EXISTS(); #ifdef __CLASSIC_C__ int main(){ int ac; - uint8_t*av[]; + char*av[]; #else -int main(int ac, uint8_t*av[]){ - +int main(int ac, char*av[]){ + #endif float ac2 = sqrtf(rand()); CHECK_FUNCTION_EXISTS(); diff --git a/cmake/modules/FindMATLAB.cmake b/cmake/modules/FindMATLAB.cmake new file mode 100644 index 000000000..a8218f8ba --- /dev/null +++ b/cmake/modules/FindMATLAB.cmake @@ -0,0 +1,212 @@ +# - this module looks for Matlab +# Defines: +# MATLAB_INCLUDE_DIR: include path for mex.h, engine.h +# MATLAB_LIBRARIES: required libraries: libmex, etc +# MATLAB_MEX_LIBRARY: path to libmex.lib +# MATLAB_MX_LIBRARY: path to libmx.lib +# MATLAB_MAT_LIBRARY: path to libmat.lib # added +# MATLAB_ENG_LIBRARY: path to libeng.lib +# MATLAB_ROOT: path to Matlab's root directory + +# This file is part of Gerardus +# +# This is a derivative work of file FindMatlab.cmake released with +# CMake v2.8, because the original seems to be a bit outdated and +# doesn't work with my Windows XP and Visual Studio 10 install +# +# (Note that the original file does work for Ubuntu Natty) +# +# Author: Ramon Casero , Tom Doel +# Version: 0.2.3 +# $Rev$ +# $Date$ +# +# The original file was copied from an Ubuntu Linux install +# /usr/share/cmake-2.8/Modules/FindMatlab.cmake + +set(MATLAB_FOUND 0) +if(WIN32) + # Search for a version of Matlab available, starting from the most modern one to older versions + foreach(MATVER "7.14" "7.11" "7.10" "7.9" "7.8" "7.7" "7.6" "7.5" "7.4") + if((NOT DEFINED MATLAB_ROOT) + OR ("${MATLAB_ROOT}" STREQUAL "") + OR ("${MATLAB_ROOT}" STREQUAL "/registry")) + get_filename_component(MATLAB_ROOT "[HKEY_LOCAL_MACHINE\\SOFTWARE\\MathWorks\\MATLAB\\${MATVER};MATLABROOT]" + ABSOLUTE) + set(MATLAB_VERSION ${MATVER}) + endif() + endforeach() + + # Directory name depending on whether the Windows architecture is 32 + # bit or 64 bit + set(CMAKE_SIZEOF_VOID_P 8) # Note: For some wierd reason this variable is undefined in my system... + if(CMAKE_SIZEOF_VOID_P MATCHES "4") + set(WINDIR "win32") + elseif(CMAKE_SIZEOF_VOID_P MATCHES "8") + set(WINDIR "win64") + else() + message(FATAL_ERROR "CMAKE_SIZEOF_VOID_P (${CMAKE_SIZEOF_VOID_P}) doesn't indicate a valid platform") + endif() + + # Folder where the MEX libraries are, depending of the Windows compiler + if(${CMAKE_GENERATOR} MATCHES "Visual Studio 6") + set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc60") + elseif(${CMAKE_GENERATOR} MATCHES "Visual Studio 7") + # Assume people are generally using Visual Studio 7.1, + # if using 7.0 need to link to: ../extern/lib/${WINDIR}/microsoft/msvc70 + set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc71") + # set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc70") + elseif(${CMAKE_GENERATOR} MATCHES "Borland") + # Assume people are generally using Borland 5.4, + # if using 7.0 need to link to ../extern/lib/${WINDIR}/microsoft/msvc70 + set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc54") + # set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc50") + # set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc51") + elseif(${CMAKE_GENERATOR} MATCHES "Visual Studio*") + # If the compiler is Visual Studio, but not any of the specific + # versions above, we try our luck with the microsoft directory + set(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/") + else() + message(FATAL_ERROR "Generator not compatible: ${CMAKE_GENERATOR}") + endif() + + # Get paths to the Matlab MEX libraries + find_library(MATLAB_MEX_LIBRARY libmex ${MATLAB_LIBRARIES_DIR} ) + find_library(MATLAB_MX_LIBRARY libm ${MATLAB_LIBRARIES_DIR} ) + find_library(MATLAB_MAT_LIBRARY libmat ${MATLAB_LIBRARIES_DIR} ) + find_library(MATLAB_ENG_LIBRARY libeng ${MATLAB_LIBRARIES_DIR} ) + + # Get path to the include directory + find_path(MATLAB_INCLUDE_DIR "mex.h" "${MATLAB_ROOT}/extern/include" ) + +else() + + if((NOT DEFINED MATLAB_ROOT) + OR ("${MATLAB_ROOT}" STREQUAL "")) + # get path to the Matlab root directory + execute_process( + COMMAND which matlab + COMMAND xargs readlink + COMMAND xargs dirname + COMMAND xargs dirname + COMMAND xargs echo -n + OUTPUT_VARIABLE MATLAB_ROOT + ) + endif() + + # Check if this is a Mac + if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") + + set(LIBRARY_EXTENSION .dylib) + + # If this is a Mac and the attempts to find MATLAB_ROOT have so far failed, + # we look in the applications folder + if((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL "")) + + # Search for a version of Matlab available, starting from the most modern one to older versions + foreach(MATVER "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b") + if((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL "")) + if(EXISTS /Applications/MATLAB_${MATVER}.app) + set(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app) + endif() + endif() + endforeach() + endif() + + else() + set(LIBRARY_EXTENSION .so) + endif() + + # Get path to the MEX libraries + execute_process( + #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmex${LIBRARY_EXTENSION} # Peter + COMMAND find "${MATLAB_ROOT}/bin" -name libmex${LIBRARY_EXTENSION} # standard + COMMAND xargs echo -n + OUTPUT_VARIABLE MATLAB_MEX_LIBRARY + ) + execute_process( + #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmx${LIBRARY_EXTENSION} # Peter + COMMAND find "${MATLAB_ROOT}/bin" -name libmx${LIBRARY_EXTENSION} # Standard + COMMAND xargs echo -n + OUTPUT_VARIABLE MATLAB_MX_LIBRARY + ) + execute_process( + #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmat${LIBRARY_EXTENSION} # Peter + COMMAND find "${MATLAB_ROOT}/bin" -name libmat${LIBRARY_EXTENSION} # Standard + COMMAND xargs echo -n + OUTPUT_VARIABLE MATLAB_MAT_LIBRARY + ) + execute_process( + #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libeng${LIBRARY_EXTENSION} # Peter + COMMAND find "${MATLAB_ROOT}/bin" -name libeng${LIBRARY_EXTENSION} # Standard + COMMAND xargs echo -n + OUTPUT_VARIABLE MATLAB_ENG_LIBRARY + ) + + # Get path to the include directory + find_path(MATLAB_INCLUDE_DIR + "mex.h" + PATHS "${MATLAB_ROOT}/extern/include" + ) + + find_program( MATLAB_MEX_PATH mex + HINTS ${MATLAB_ROOT}/bin + PATHS ${MATLAB_ROOT}/bin + DOC "The mex program path" + ) + + find_program( MATLAB_MEXEXT_PATH mexext + HINTS ${MATLAB_ROOT}/bin + PATHS ${MATLAB_ROOT}/bin + DOC "The mexext program path" + ) + + execute_process( + COMMAND ${MATLAB_MEXEXT_PATH} + OUTPUT_STRIP_TRAILING_WHITESPACE + OUTPUT_VARIABLE MATLAB_MEX_EXT + ) + +endif() + +# This is common to UNIX and Win32: +set(MATLAB_LIBRARIES + ${MATLAB_MEX_LIBRARY} + ${MATLAB_MX_LIBRARY} + ${MATLAB_ENG_LIBRARY} +) + +if(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES) + set(MATLAB_FOUND 1) +endif() + +# 32-bit or 64-bit mex +if(WIN32) + if (CMAKE_CL_64) + SET(MATLAB_MEX_EXTENSION .mexw64) + else(CMAKE_CL_64) + SET(MATLAB_MEX_EXTENSION .mexw32) + endif(CMAKE_CL_64) +else(WIN32) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET(MATLAB_MEX_EXTENSION .mexa64) + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET(MATLAB_MEX_EXTENSION .mexglx) + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") +endif(WIN32) + +SET(MATLAB_MEX_CFLAGS "-DMATLAB_MEX_FILE -DMX_COMPAT_32") + +mark_as_advanced( + MATLAB_LIBRARIES + MATLAB_MEX_LIBRARY + MATLAB_MX_LIBRARY + MATLAB_ENG_LIBRARY + MATLAB_INCLUDE_DIR + MATLAB_FOUND + MATLAB_ROOT + MATLAB_MEX_PATH + MATLAB_MEXEXT_PATH + MATLAB_MEX_EXT +) + diff --git a/cmake/modules/FindOCTAVE.cmake b/cmake/modules/FindOCTAVE.cmake new file mode 100644 index 000000000..63a43d740 --- /dev/null +++ b/cmake/modules/FindOCTAVE.cmake @@ -0,0 +1,138 @@ +# - Try to find a version of Octave and headers/library required by the +# used compiler. It determines the right MEX-File extension and add +# a macro to help the build of MEX-functions. +# +# This module defines: +# OCTAVE_INCLUDE_DIR: include path for mex.h, mexproto.h +# OCTAVE_OCTINTERP_LIBRARY: path to the library octinterp +# OCTAVE_OCTAVE_LIBRARY: path to the library octave +# OCTAVE_CRUFT_LIBRARY: path to the library cruft +# OCTAVE_LIBRARIES: required libraries: octinterp, octave, cruft +# OCTAVE_CREATE_MEX: macro to build a MEX-file +# +# The macro OCTAVE_CREATE_MEX requires in this order: +# - function's name which will be called in Octave; +# - C/C++ source files; +# - third libraries required. + +# Copyright (c) 2009-2013 Arnaud Barré +# Redistribution and use is allowed according to the terms of the BSD license. +# For details see the accompanying COPYING-CMAKE-SCRIPTS file. + +IF(OCTAVE_ROOT AND OCTAVE_INCLUDE_DIR AND OCTAVE_LIBRARIES) + STRING(COMPARE NOTEQUAL "${OCTAVE_ROOT}" "${OCTAVE_ROOT_LAST}" OCTAVE_CHANGED) + IF(OCTAVE_CHANGED) + SET(OCTAVE_USE_MINGW32 OCTAVE_USE_MINGW32-NOTFOUND CACHE INTERNAL "") + SET(OCTAVE_OCTAVE_LIBRARY OCTAVE_OCTAVE_LIBRARY-NOTFOUND CACHE INTERNAL "") + SET(OCTAVE_INCLUDE_DIR OCTAVE_INCLUDE_DIR-NOTFOUND CACHE INTERNAL "") + ELSE(OCTAVE_CHANGED) + # in cache already + SET(Octave_FIND_QUIETLY TRUE) + ENDIF(OCTAVE_CHANGED) +ENDIF(OCTAVE_ROOT AND OCTAVE_INCLUDE_DIR AND OCTAVE_LIBRARIES) + +SET(OCTAVE_MEXFILE_EXT mex) + +IF(WIN32) + SET(OCTAVE_PATHS_L1 ) + SET(OCTAVE_PATHS_L2 ) + # Level 0 + FILE(GLOB OCTAVE_PATHS_L0 "c:/Octave*") + # Level 1 + FOREACH(_file_ ${OCTAVE_PATHS_L0}) + FILE(GLOB OCTAVE_PATHS_TEMP "${_file_}/*") + SET(OCTAVE_PATHS_L1 ${OCTAVE_PATHS_L1};${OCTAVE_PATHS_TEMP}) + ENDFOREACH(_file_ OCTAVE_PATHS_L0) + # Level 2 + FOREACH(_file_ ${OCTAVE_PATHS_L1}) + FILE(GLOB OCTAVE_PATHS_TEMP "${_file_}/*") + SET(OCTAVE_PATHS_L2 ${OCTAVE_PATHS_L2};${OCTAVE_PATHS_TEMP}) + ENDFOREACH(_file_ OCTAVE_PATHS_L1) + # Merge levels + SET(OCTAVE_PATHS ${OCTAVE_PATHS_L0} ${OCTAVE_PATHS_L1} ${OCTAVE_PATHS_L2}) + + FIND_PATH(OCTAVE_ROOT "bin/octave.exe" ${OCTAVE_PATHS}) + FIND_PATH(OCTAVE_USE_MINGW32 "bin/mingw32-make.exe" "${OCTAVE_ROOT}/mingw32") + + IF(MSVC AND OCTAVE_USE_MINGW32) + MESSAGE(FATAL_ERROR + "You must use the generator \"MinGW Makefiles\" as the " + "version of Octave installed on your computer was compiled " + "with MinGW. You should also specify the native compiler " + "(GCC, G++ and GFortan) and add the path of MinGW in the " + "environment variable PATH. Contact the developers of the " + "project for more details") + ENDIF(MSVC AND OCTAVE_USE_MINGW32) + + FILE(GLOB OCTAVE_INCLUDE_PATHS "${OCTAVE_ROOT}/include/octave-*/octave") + FILE(GLOB OCTAVE_LIBRARIES_PATHS "${OCTAVE_ROOT}/lib/octave-*") + IF (NOT OCTAVE_LIBRARIES_PATHS) + FILE(GLOB OCTAVE_LIBRARIES_PATHS "${OCTAVE_ROOT}/lib/octave/*") + ENDIF (NOT OCTAVE_LIBRARIES_PATHS) + + # LIBOCTINTERP, LIBOCTAVE, LIBCRUFT names + SET(LIBOCTAVE "liboctave") + +ELSE(WIN32) + IF(APPLE) + FILE(GLOB OCTAVE_PATHS "/Applications/Octave*") + FIND_PATH(OCTAVE_ROOT "Contents/Resources/bin/octave" ${OCTAVE_PATHS}) + + FILE(GLOB OCTAVE_INCLUDE_PATHS "${OCTAVE_ROOT}/Contents/Resources/include/octave-*/octave") + FILE(GLOB OCTAVE_LIBRARIES_PATHS "${OCTAVE_ROOT}/Contents/Resources/lib/octave-*") + + SET(LIBOCTAVE "liboctave.dylib") + ELSE(APPLE) + FILE(GLOB OCTAVE_LOCAL_PATHS "/usr/local/lib/octave-*") + FILE(GLOB OCTAVE_USR_PATHS "/usr/lib/octave-*") + FILE(GLOB OCTAVE_INCLUDE_PATHS "/usr/include/octave-*") + + SET (OCTAVE_INCLUDE_PATHS + "/usr/local/include" + "/usr/local/include/octave" + "/usr/include" + "${OCTAVE_INCLUDE_PATHS}" + "${OCTAVE_INCLUDE_PATHS}/octave") + SET (OCTAVE_LIBRARIES_PATHS + "/usr/local/lib" + "/usr/local/lib/octave" + ${OCTAVE_LOCAL_PATHS} + "/usr/lib" + "/usr/lib/octave" + ${OCTAVE_USR_PATHS}) + + SET (LIBOCTAVE "octave") + ENDIF(APPLE) +ENDIF(WIN32) + +FIND_LIBRARY(OCTAVE_OCTAVE_LIBRARY + ${LIBOCTAVE} + ${OCTAVE_LIBRARIES_PATHS} + ) +FIND_PATH(OCTAVE_INCLUDE_DIR + "mex.h" + ${OCTAVE_INCLUDE_PATHS} + ) + +SET(OCTAVE_ROOT_LAST "${OCTAVE_ROOT}" CACHE INTERNAL "" FORCE) + +# This is common to UNIX and Win32: +SET(OCTAVE_LIBRARIES + ${OCTAVE_OCTAVE_LIBRARY} + CACHE INTERNAL "Octave libraries" FORCE +) + +INCLUDE(FindPackageHandleStandardArgs) + +# The variable OCTAVE_ROOT is only relevant for WIN32 +IF(WIN32) + FIND_PACKAGE_HANDLE_STANDARD_ARGS(Octave DEFAULT_MSG OCTAVE_ROOT OCTAVE_INCLUDE_DIR OCTAVE_OCTAVE_LIBRARY ) +ELSE(WIN32) + FIND_PACKAGE_HANDLE_STANDARD_ARGS(Octave DEFAULT_MSG OCTAVE_INCLUDE_DIR OCTAVE_OCTAVE_LIBRARY ) +ENDIF(WIN32) + +MARK_AS_ADVANCED( + OCTAVE_OCTAVE_LIBRARY + OCTAVE_LIBRARIES + OCTAVE_INCLUDE_DIR +) \ No newline at end of file diff --git a/cmake/modules/FindVolk.cmake b/cmake/modules/FindVolk.cmake index 5861f1232..eff504408 100644 --- a/cmake/modules/FindVolk.cmake +++ b/cmake/modules/FindVolk.cmake @@ -35,6 +35,7 @@ CHECK_FUNCTION_EXISTS_MATH(volk_32fc_x2_multiply_conjugate_32fc HAVE_VOLK_MULT2_ CHECK_FUNCTION_EXISTS_MATH(volk_32fc_32f_multiply_32fc HAVE_VOLK_MULT_REAL_FUNCTION) CHECK_FUNCTION_EXISTS_MATH(volk_32f_s32f_multiply_32f HAVE_VOLK_MULT_FLOAT_FUNCTION) CHECK_FUNCTION_EXISTS_MATH(volk_32fc_magnitude_32f HAVE_VOLK_MAG_FUNCTION) +CHECK_FUNCTION_EXISTS_MATH(volk_32fc_magnitude_square_32f HAVE_VOLK_MAG_SQUARE_FUNCTION) CHECK_FUNCTION_EXISTS_MATH(volk_32f_x2_divide_32f HAVE_VOLK_DIVIDE_FUNCTION) CHECK_FUNCTION_EXISTS_MATH(volk_32fc_32f_dot_prod_32fc HAVE_VOLK_DOTPROD_FC_FUNCTION) CHECK_FUNCTION_EXISTS_MATH(volk_32fc_x2_conjugate_dot_prod_32fc HAVE_VOLK_DOTPROD_CONJ_FC_FUNCTION) diff --git a/lte/phy/include/liblte/phy/ch_estimation/chest_dl.h b/lte/phy/include/liblte/phy/ch_estimation/chest_dl.h new file mode 100644 index 000000000..c1f3d23ca --- /dev/null +++ b/lte/phy/include/liblte/phy/ch_estimation/chest_dl.h @@ -0,0 +1,105 @@ +/** + * + * \section COPYRIGHT + * + * Copyright 2013-2014 The libLTE Developers. See the + * COPYRIGHT file at the top-level directory of this distribution. + * + * \section LICENSE + * + * This file is part of the libLTE library. + * + * libLTE is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * libLTE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * A copy of the GNU Lesser General Public License can be found in + * the LICENSE file in the top-level directory of this distribution + * and at http://www.gnu.org/licenses/. + * + */ + + + +#ifndef CHEST_DL_ +#define CHEST_DL_ + +#include + +#include "liblte/config.h" + +#include "liblte/phy/resampling/interp.h" +#include "liblte/phy/filter/filter2d.h" +#include "liblte/phy/ch_estimation/refsignal_dl.h" +#include "liblte/phy/common/phy_common.h" + +#define CHEST_RS_AVERAGE_TIME 0 +#define CHEST_RS_AVERAGE_FREQ 3 + + + +/** 3GPP LTE Downlink channel estimator and equalizer. + * Estimates the channel in the resource elements transmitting references and interpolates for the rest + * of the resource grid. + * + * The equalizer uses the channel estimates to produce an estimation of the transmitted symbol. + * + * This object depends on the refsignal_t object for creating the LTE CSR signal. +*/ + +typedef struct { + lte_cell_t cell; + refsignal_cs_t csr_signal; + cf_t *pilot_estimates[MAX_PORTS]; + cf_t *pilot_recv_signal[MAX_PORTS]; + cf_t *pilot_symbol_avg; + + interp_t interp_time[MAX_PORTS]; + interp_t interp_freq[MAX_PORTS]; + + float rssi; + float rsrq; + float rsrp; +} chest_dl_t; + + +LIBLTE_API int chest_dl_init(chest_dl_t *q, + lte_cell_t cell); + +LIBLTE_API void chest_dl_free(chest_dl_t *q); + +LIBLTE_API int chest_dl_estimate(chest_dl_t *q, + cf_t *input, + cf_t *ce[MAX_PORTS], + uint32_t sf_idx); + +LIBLTE_API int chest_dl_estimate_port(chest_dl_t *q, + cf_t *input, + cf_t *ce, + uint32_t sf_idx, + uint32_t port_id); + +LIBLTE_API int chest_dl_equalize_zf(chest_dl_t *q, + cf_t *input, + cf_t *ce[MAX_PORTS], + cf_t *output); + +LIBLTE_API int chest_dl_equalize_mmse(chest_dl_t *q, + cf_t *input, + cf_t *ce[MAX_PORTS], + float *noise_estimate, + cf_t *output); + +LIBLTE_API float chest_dl_get_rssi(chest_dl_t *q); + +LIBLTE_API float chest_dl_get_rsrq(chest_dl_t *q); + +LIBLTE_API float chest_dl_get_rsrp(chest_dl_t *q); + +#endif \ No newline at end of file diff --git a/lte/phy/include/liblte/phy/ch_estimation/refsignal.h b/lte/phy/include/liblte/phy/ch_estimation/refsignal.h index e237c021a..6efb85fa9 100644 --- a/lte/phy/include/liblte/phy/ch_estimation/refsignal.h +++ b/lte/phy/include/liblte/phy/ch_estimation/refsignal.h @@ -87,6 +87,6 @@ LIBLTE_API int refsignal_init_LTEUL_drms_pusch(refsignal_t *q, LIBLTE_API void refsignal_free(refsignal_t *q); LIBLTE_API int refsignal_put(refsignal_t *q, - cf_t *slot_symbols); + cf_t *slot_symbols); #endif diff --git a/lte/phy/include/liblte/phy/ch_estimation/refsignal_dl.h b/lte/phy/include/liblte/phy/ch_estimation/refsignal_dl.h new file mode 100644 index 000000000..6b8374572 --- /dev/null +++ b/lte/phy/include/liblte/phy/ch_estimation/refsignal_dl.h @@ -0,0 +1,87 @@ +/** + * + * \section COPYRIGHT + * + * Copyright 2013-2014 The libLTE Developers. See the + * COPYRIGHT file at the top-level directory of this distribution. + * + * \section LICENSE + * + * This file is part of the libLTE library. + * + * libLTE is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * libLTE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * A copy of the GNU Lesser General Public License can be found in + * the LICENSE file in the top-level directory of this distribution + * and at http://www.gnu.org/licenses/. + * + */ + +#ifndef REFSIGNAL_DL_ +#define REFSIGNAL_DL_ + +/* Object to manage Downlink reference signals for channel estimation. + * + */ + +#include "liblte/config.h" +#include "liblte/phy/common/phy_common.h" + +typedef _Complex float cf_t; + +// Number of references in a subframe: there are 2 symbols for port_id=0,1 x 2 slots x 2 refs per prb +#define REFSIGNAL_NUM_SF(nof_prb, port_id) (((port_id)<2?8:4)*(nof_prb)) +#define REFSIGNAL_MAX_NUM_SF(nof_prb) REFSIGNAL_NUM_SF(nof_prb, 0) + +#define REFSIGNAL_PILOT_IDX(i,l,ns,cell) (2*cell.nof_prb*((l)+2*((ns)%2))+(i)) + + +/** Cell-Specific Reference Signal */ +typedef struct LIBLTE_API { + lte_cell_t cell; + cf_t *pilots[NSUBFRAMES_X_FRAME]; // save 2 reference symbols per slot +} refsignal_cs_t; + + +LIBLTE_API int refsignal_cs_generate(refsignal_cs_t *q, + lte_cell_t cell); + +LIBLTE_API void refsignal_cs_free(refsignal_cs_t *q); + +LIBLTE_API int refsignal_cs_put_sf(lte_cell_t cell, + uint32_t port_id, + uint32_t sf_idx, + cf_t *pilots, + cf_t *sf_symbols); + +LIBLTE_API int refsignal_cs_get_sf(lte_cell_t cell, + uint32_t port_id, + uint32_t sf_idx, + cf_t *sf_symbols, + cf_t *pilots); + +LIBLTE_API uint32_t refsignal_fidx(lte_cell_t cell, + uint32_t ns, + uint32_t l, + uint32_t port_id, + uint32_t m); + +LIBLTE_API uint32_t refsignal_nsymbol(lte_cell_t cell, + uint32_t ns, + uint32_t l); + +LIBLTE_API uint32_t refsignal_cs_v(uint32_t port_id, + uint32_t ns, + uint32_t symbol_id); + +LIBLTE_API uint32_t refsignal_cs_nof_symbols(uint32_t port_id); + +#endif diff --git a/lte/phy/include/liblte/phy/common/phy_common.h b/lte/phy/include/liblte/phy/common/phy_common.h index 5a4b80419..a1d7e6dd8 100644 --- a/lte/phy/include/liblte/phy/common/phy_common.h +++ b/lte/phy/include/liblte/phy/common/phy_common.h @@ -146,6 +146,14 @@ LIBLTE_API bool lte_cell_isvalid(lte_cell_t *cell); LIBLTE_API void lte_cell_fprint(FILE *stream, lte_cell_t *cell); +LIBLTE_API bool lte_cellid_isvalid(uint32_t cell_id); + +LIBLTE_API bool lte_nofprb_isvalid(uint32_t nof_prb); + +LIBLTE_API bool lte_sfidx_isvalid(uint32_t sf_idx); + +LIBLTE_API bool lte_portid_isvalid(uint32_t port_id); + LIBLTE_API bool lte_N_id_2_isvalid(uint32_t N_id_2); LIBLTE_API bool lte_N_id_1_isvalid(uint32_t N_id_1); diff --git a/lte/phy/include/liblte/phy/common/sequence.h b/lte/phy/include/liblte/phy/common/sequence.h index 183d5cddc..a220b230a 100644 --- a/lte/phy/include/liblte/phy/common/sequence.h +++ b/lte/phy/include/liblte/phy/common/sequence.h @@ -44,6 +44,9 @@ LIBLTE_API int sequence_LTE_pr(sequence_t *q, uint32_t len, uint32_t seed); +LIBLTE_API void sequence_set_LTE_pr(sequence_t *q, + uint32_t seed); + LIBLTE_API int sequence_pbch(sequence_t *seq, lte_cp_t cp, uint32_t cell_id); diff --git a/lte/phy/include/liblte/phy/filter/filter2d.h b/lte/phy/include/liblte/phy/filter/filter2d.h index afbd3b6f1..a3978145f 100644 --- a/lte/phy/include/liblte/phy/filter/filter2d.h +++ b/lte/phy/include/liblte/phy/filter/filter2d.h @@ -31,6 +31,7 @@ #define FILTER2D_ #include "liblte/config.h" +#include /* 2-D real filter of complex input * @@ -38,18 +39,42 @@ typedef _Complex float cf_t; typedef struct LIBLTE_API{ - int sztime; // Output signal size in the time domain - int szfreq; // Output signal size in the freq domain - int ntime; // 2-D Filter size in time domain - int nfreq; // 2-D Filter size in frequency domain + uint32_t sztime; // Output signal size in the time domain + uint32_t szfreq; // Output signal size in the freq domain + uint32_t ntime; // 2-D Filter size in time domain + uint32_t nfreq; // 2-D Filter size in frequency domain float **taps; // 2-D filter coefficients + float norm; //normalization factor cf_t *output; // Output signal } filter2d_t; -LIBLTE_API int filter2d_init (filter2d_t* q, float **taps, int ntime, int nfreq, int sztime, int szfreq); -LIBLTE_API int filter2d_init_default (filter2d_t* q, int ntime, int nfreq, int sztime, int szfreq); +LIBLTE_API int filter2d_init (filter2d_t* q, + float **taps, + uint32_t ntime, + uint32_t nfreq, + uint32_t sztime, + uint32_t szfreq); + +LIBLTE_API int filter2d_init_ones (filter2d_t* q, + uint32_t ntime, + uint32_t nfreq, + uint32_t sztime, + uint32_t szfreq); + LIBLTE_API void filter2d_free(filter2d_t *q); + +LIBLTE_API void filter2d_step(filter2d_t *q); + LIBLTE_API void filter2d_reset(filter2d_t *q); -LIBLTE_API void filter2d_add(filter2d_t *q, cf_t h, int time_idx, int freq_idx); +LIBLTE_API void filter2d_add(filter2d_t *q, + cf_t h, + uint32_t time_idx, + uint32_t freq_idx); + +LIBLTE_API void filter2d_add_out(filter2d_t *q, cf_t x, int time_idx, int freq_idx); + +LIBLTE_API void filter2d_get_symbol(filter2d_t *q, + uint32_t nsymbol, + cf_t *output); #endif // FILTER2D_ diff --git a/lte/phy/include/liblte/phy/phy.h b/lte/phy/include/liblte/phy/phy.h index eef9eccc1..ed8152bd5 100644 --- a/lte/phy/include/liblte/phy/phy.h +++ b/lte/phy/include/liblte/phy/phy.h @@ -52,7 +52,9 @@ #include "liblte/phy/common/fft.h" #include "liblte/phy/ch_estimation/chest.h" +#include "liblte/phy/ch_estimation/chest_dl.h" #include "liblte/phy/ch_estimation/refsignal.h" +#include "liblte/phy/ch_estimation/refsignal_dl.h" #include "liblte/phy/resampling/interp.h" #include "liblte/phy/resampling/decim.h" diff --git a/lte/phy/include/liblte/phy/resampling/interp.h b/lte/phy/include/liblte/phy/resampling/interp.h index c39c83602..e2028b22b 100644 --- a/lte/phy/include/liblte/phy/resampling/interp.h +++ b/lte/phy/include/liblte/phy/resampling/interp.h @@ -78,6 +78,10 @@ LIBLTE_API void interp_run_offset(interp_t *q, uint32_t off_st, uint32_t off_end); +LIBLTE_API cf_t interp_linear_onesample(cf_t *input); + +LIBLTE_API cf_t interp_linear_onesample2(cf_t *input); + LIBLTE_API void interp_linear_offset(cf_t *input, cf_t *output, uint32_t M, diff --git a/lte/phy/include/liblte/phy/utils/vector.h b/lte/phy/include/liblte/phy/utils/vector.h index ee59bfad8..27f3fcc19 100644 --- a/lte/phy/include/liblte/phy/utils/vector.h +++ b/lte/phy/include/liblte/phy/utils/vector.h @@ -96,7 +96,10 @@ LIBLTE_API cf_t vec_dot_prod_conj_ccc(cf_t *x, cf_t *y, uint32_t len); LIBLTE_API float vec_dot_prod_fff(float *x, float *y, uint32_t len); /* z=x/y vector division (element-wise) */ -LIBLTE_API void vec_div_ccc(cf_t *x, cf_t *y, cf_t *z, uint32_t len); +LIBLTE_API void vec_div_ccc(cf_t *x, cf_t *y, float *y_mod, cf_t *z, uint32_t len); + +/* z=x/y vector division with mod(y)=1 (element-wise) */ +LIBLTE_API void vec_div_ccc_mod1(cf_t *x, cf_t *y, cf_t *z, uint32_t len); /* conjugate */ LIBLTE_API void vec_conj_cc(cf_t *x, cf_t *y, uint32_t len); @@ -113,6 +116,7 @@ LIBLTE_API void vec_quant_fuc(float *in, uint8_t *out, float gain, float offset, /* magnitude of each vector element */ LIBLTE_API void vec_abs_cf(cf_t *x, float *abs, uint32_t len); +LIBLTE_API void vec_abs_square_cf(cf_t *x, float *abs_square, uint32_t len); /* argument of each vector element */ LIBLTE_API void vec_arg_cf(cf_t *x, float *arg, uint32_t len); diff --git a/lte/phy/lib/ch_estimation/src/chest.c b/lte/phy/lib/ch_estimation/src/chest.c index b34471122..417695899 100644 --- a/lte/phy/lib/ch_estimation/src/chest.c +++ b/lte/phy/lib/ch_estimation/src/chest.c @@ -38,7 +38,7 @@ #define SLOT_SZ(q) (q->nof_symbols * q->symbol_sz) #define SF_SZ(q) (2 * SLOT_SZ(q)) -//#define VOLK_INTERP +#define VOLK_INTERP void chest_fprint(chest_t *q, FILE *stream, uint32_t nslot, uint32_t port_id) { chest_ref_fprint(q, stream, nslot, port_id); diff --git a/lte/phy/lib/ch_estimation/src/chest_dl.c b/lte/phy/lib/ch_estimation/src/chest_dl.c new file mode 100644 index 000000000..d515fc413 --- /dev/null +++ b/lte/phy/lib/ch_estimation/src/chest_dl.c @@ -0,0 +1,317 @@ +/** + * + * \section COPYRIGHT + * + * Copyright 2013-2014 The libLTE Developers. See the + * COPYRIGHT file at the top-level directory of this distribution. + * + * \section LICENSE + * + * This file is part of the libLTE library. + * + * libLTE is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * libLTE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * A copy of the GNU Lesser General Public License can be found in + * the LICENSE file in the top-level directory of this distribution + * and at http://www.gnu.org/licenses/. + * + */ + + + +#include +#include +#include +#include +#include + +#include "liblte/config.h" + +#include "liblte/phy/ch_estimation/chest_dl.h" +#include "liblte/phy/utils/vector.h" + +//#define VOLK_INTERP + +/** 3GPP LTE Downlink channel estimator and equalizer. + * Estimates the channel in the resource elements transmitting references and interpolates for the rest + * of the resource grid. + * + * The equalizer uses the channel estimates to produce an estimation of the transmitted symbol. + * + * This object depends on the refsignal_t object for creating the LTE CSR signal. +*/ + +int chest_dl_init(chest_dl_t *q, lte_cell_t cell) +{ + int ret = LIBLTE_ERROR_INVALID_INPUTS; + if (q != NULL && + lte_cell_isvalid(&cell)) + { + bzero(q, sizeof(chest_dl_t)); + + ret = refsignal_cs_generate(&q->csr_signal, cell); + if (ret != LIBLTE_SUCCESS) { + fprintf(stderr, "Error initializing CSR signal (%d)\n",ret); + goto clean_exit; + } + + q->pilot_symbol_avg = vec_malloc(sizeof(cf_t) * 2*cell.nof_prb); + if (!q->pilot_symbol_avg) { + perror("malloc"); + goto clean_exit; + } + + for (int i=0;ipilot_estimates[i] = vec_malloc(sizeof(cf_t) * REFSIGNAL_MAX_NUM_SF(cell.nof_prb)); + if (!q->pilot_estimates[i]) { + perror("malloc"); + goto clean_exit; + } + q->pilot_recv_signal[i] = vec_malloc(sizeof(cf_t) * REFSIGNAL_MAX_NUM_SF(cell.nof_prb)); + if (!q->pilot_recv_signal[i]) { + perror("malloc"); + goto clean_exit; + } + #ifdef VOLK_INTERP + ret = interp_init(&q->interp_freq[i], LINEAR, 2*cell.nof_prb, RE_X_RB/2); + if (ret == LIBLTE_SUCCESS) { + ret = interp_init(&q->interp_time[i], LINEAR, 2, CP_NSYMB(cell.cp) - 3); + } + #endif + + } + + /* Init buffer for holding CE estimates averages */ + + q->cell = cell; + } + + ret = LIBLTE_SUCCESS; + +clean_exit: + if (ret != LIBLTE_SUCCESS) { + chest_dl_free(q); + } + return ret; +} + +void chest_dl_free(chest_dl_t *q) +{ + refsignal_cs_free(&q->csr_signal); + + if (q->pilot_symbol_avg) { + free(q->pilot_symbol_avg); + } + + for (int i=0;ipilot_estimates[i]) { + free(q->pilot_estimates[i]); + } + if (q->pilot_recv_signal[i]) { + free(q->pilot_recv_signal[i]); + } + #ifdef VOLK_INTERP + interp_free(&q->interp_freq[i]); + interp_free(&q->interp_time[i]); + #endif + } +} + +#define pilot_est(idx) q->pilot_estimates[port_id][REFSIGNAL_PILOT_IDX(idx,l,ns,q->cell)] + +#if CHEST_RS_AVERAGE_TIME > 1 +cf_t timeavg[CHEST_RS_AVERAGE_TIME-1][12]; +int nof_timeavg=0; +#endif + +static void average_pilots(chest_dl_t *q, uint32_t sf_idx, uint32_t port_id) +{ + uint32_t ns, l; + int i; + /* For each symbol with pilots in a slot */ + for (ns=2*sf_idx;ns<2*(sf_idx+1);ns++) { + for (l=0;lpilot_symbol_avg, 2*q->cell.nof_prb); + + /** Frequency average */ +#if CHEST_RS_AVERAGE_FREQ > 1 + const uint32_t M = CHEST_RS_AVERAGE_FREQ; + cf_t xint[CHEST_RS_AVERAGE_FREQ]; + int j, k; + /* Extrapolate first M/2 samples */ + for (i=M/2-1;i>=0;i--) { + k=0; + for (j=i+M/2;j>=0;j--) { + xint[k]=pilot_est(j); + k++; + } + for (;j>=i-M/2;j--) { + if (k>=2) { + xint[k] = interp_linear_onesample(&xint[k-2]); + k++; + } + } + q->pilot_symbol_avg[i] = vec_acc_cc(xint,M)/M; + //q->pilot_symbol_avg[i] = (pilot_est(0)+pilot_est(1))/2; + } + + for (i=M/2;i<2*q->cell.nof_prb-M/2;i++) { + q->pilot_symbol_avg[i] = vec_acc_cc(&pilot_est(i-M/2),M)/M; + } + + /* Extrapolate last M/2 samples */ + for (;i<2*q->cell.nof_prb;i++) { + k=0; + for (j=i-M/2;j<2*q->cell.nof_prb;j++) { + xint[k]=pilot_est(j); + k++; + } + for (;k=2) { + xint[k] = interp_linear_onesample(&xint[k-2]); + } + } + q->pilot_symbol_avg[i] = vec_acc_cc(xint,M)/M; + //q->pilot_symbol_avg[i] = (pilot_est(i)+pilot_est(i+1))/2; + + } +#else + memcpy(q->pilot_symbol_avg, &pilot_est(0), 2*q->cell.nof_prb*sizeof(cf_t)); +#endif + + /* Time average last symbols */ +#if CHEST_RS_AVERAGE_TIME > 1 + if (nof_timeavgpilot_symbol_avg, 2*q->cell.nof_prb * sizeof(cf_t)); + nof_timeavg++; + } else { + bzero(&pilot_est(0),2*q->cell.nof_prb*sizeof(cf_t)); + for (i=0;icell.nof_prb); + } + vec_sum_ccc(q->pilot_symbol_avg,&pilot_est(0),&pilot_est(0),2*q->cell.nof_prb); + vec_sc_prod_cfc(&pilot_est(0), 1.0/CHEST_RS_AVERAGE_TIME, &pilot_est(0), 2*q->cell.nof_prb); + for (i=0;icell.nof_prb*sizeof(cf_t)); + } + memcpy(timeavg[i],q->pilot_symbol_avg,2*q->cell.nof_prb*sizeof(cf_t)); + } +#else + memcpy(&pilot_est(0), q->pilot_symbol_avg, 2*q->cell.nof_prb * sizeof(cf_t)); +#endif + + } + } + +} + +static void interpolate_pilots(chest_dl_t *q, cf_t *ce, uint32_t sf_idx, uint32_t port_id) +{ + /* interpolate the symbols with references in the freq domain */ + uint32_t ns, l, i,j; + cf_t x[2], y[MAX_NSYMB]; + + for (ns=2*sf_idx;ns<2*(sf_idx+1);ns++) { + for (l=0;lcell, ns, l, port_id, 0); +#ifdef VOLK_INTERP + interp_run_offset(&q->interp_freq[port_id], + &q->pilot_estimates[port_id][((ns%2)*2+l)*2*q->cell.nof_prb], + &ce[refsignal_nsymbol(q->cell,ns,l) * q->cell.nof_prb * RE_X_RB], + fidx_offset, RE_X_RB/2-fidx_offset); +#else + interp_linear_offset(&q->pilot_estimates[port_id][((ns%2)*2+l)*2*q->cell.nof_prb], + &ce[refsignal_nsymbol(q->cell,ns,l) * q->cell.nof_prb * RE_X_RB], RE_X_RB/2, + 2*q->cell.nof_prb, fidx_offset, RE_X_RB/2-fidx_offset); +#endif + } + } + /* now interpolate in the time domain */ + for (i=0;icell.nof_prb; i++) { + if (refsignal_cs_nof_symbols(port_id) > 1) { + for (ns=2*sf_idx;ns<2*(sf_idx+1);ns++) { + j=0; + for (l=0;lcell,ns,l) * q->cell.nof_prb * RE_X_RB + i]; + j++; + } + #ifdef VOLK_INTERP + interp_run_offset(&q->interp_time[port_id], x, y, + 0, CP_NSYMB(q->cell.cp) - 4); + #else + interp_linear_offset(x, y, CP_NSYMB(q->cell.cp) - 3, + 2, 0, CP_NSYMB(q->cell.cp) - 4); + #endif + for (j=0;jcell.cp);j++) { + ce[(j+((ns%2)*CP_NSYMB(q->cell.cp))) * q->cell.nof_prb*RE_X_RB + i] = y[j]; + } + } + } else { + fprintf(stderr, "3/4 Ports interpolator not implemented\n"); + exit(-1); + } + } +} + +int chest_dl_estimate_port(chest_dl_t *q, cf_t *input, cf_t *ce, uint32_t sf_idx, uint32_t port_id) +{ + //filter2d_reset(&q->filter); + + /* Get references from the input signal */ + refsignal_cs_get_sf(q->cell, port_id, sf_idx, input, q->pilot_recv_signal[port_id]); + + /* Use the known CSR signal to compute Least-squares estimates */ + vec_div_ccc_mod1(q->pilot_recv_signal[port_id], q->csr_signal.pilots[sf_idx], + q->pilot_estimates[port_id], REFSIGNAL_NUM_SF(q->cell.nof_prb, port_id)); + + /* Average pilot estimates */ + average_pilots(q, sf_idx, port_id); + + /* Interpolate to create channel estimates for all resource grid */ + interpolate_pilots(q, ce, sf_idx, port_id); + + return 0; +} + +int chest_dl_estimate(chest_dl_t *q, cf_t *input, cf_t *ce[MAX_PORTS], uint32_t sf_idx) +{ + uint32_t port_id; + + for (port_id=0;port_idcell.nof_ports;port_id++) { + chest_dl_estimate_port(q, input, ce[port_id], sf_idx, port_id); + } + return LIBLTE_SUCCESS; +} + +int chest_dl_equalize_zf(chest_dl_t *q, cf_t *input, cf_t *ce[MAX_PORTS], cf_t *output) +{ + fprintf(stderr, "Not implemented\n"); + return -1; +} + +int chest_dl_equalize_mmse(chest_dl_t *q, cf_t *input, cf_t *ce[MAX_PORTS], float *noise_estimate, cf_t *output) +{ + fprintf(stderr, "Not implemented\n"); + return -1; +} + +float chest_dl_get_rssi(chest_dl_t *q) { + return q->rssi; +} + +float chest_dl_get_rsrq(chest_dl_t *q) { + return q->rsrq; +} + +float chest_dl_get_rsrp(chest_dl_t *q) { + return q->rsrp; +} + diff --git a/lte/phy/lib/ch_estimation/src/refsignal_dl.c b/lte/phy/lib/ch_estimation/src/refsignal_dl.c new file mode 100644 index 000000000..e9d927bc6 --- /dev/null +++ b/lte/phy/lib/ch_estimation/src/refsignal_dl.c @@ -0,0 +1,246 @@ +/** + * + * \section COPYRIGHT + * + * Copyright 2013-2014 The libLTE Developers. See the + * COPYRIGHT file at the top-level directory of this distribution. + * + * \section LICENSE + * + * This file is part of the libLTE library. + * + * libLTE is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * libLTE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * A copy of the GNU Lesser General Public License can be found in + * the LICENSE file in the top-level directory of this distribution + * and at http://www.gnu.org/licenses/. + * + */ + + +#include +#include +#include +#include +#include + +#include "liblte/phy/common/phy_common.h" +#include "liblte/phy/ch_estimation/refsignal_dl.h" +#include "liblte/phy/utils/vector.h" +#include "liblte/phy/utils/debug.h" +#include "liblte/phy/common/sequence.h" + +uint32_t refsignal_cs_v(uint32_t port_id, uint32_t ns, uint32_t symbol_id) +{ + uint32_t v = 0; + switch (port_id) { + case 0: + if (symbol_id == 0) { + v = 0; + } else { + v = 3; + } + break; + case 1: + if (symbol_id == 0) { + v = 3; + } else { + v = 0; + } + break; + case 2: + v = 3 * (ns % 2); + break; + case 3: + v = 3 + 3 * (ns % 2); + break; + } + return v; +} + +uint32_t refsignal_cs_nof_symbols(uint32_t port_id) +{ + if (port_id < 2) { + return 2; + } else { + return 1; + } +} + +static uint32_t lp(uint32_t l, lte_cp_t cp) { + if (l == 1) { + return CP_NSYMB(cp) - 3; + } else { + return 0; + } +} + +/** Allocates and precomputes the Cell-Specific Reference (CSR) signal for + * the 20 slots in a subframe + */ +int refsignal_cs_generate(refsignal_cs_t * q, lte_cell_t cell) +{ + + uint32_t c_init; + uint32_t i, ns, l; + uint32_t N_cp, mp; + sequence_t seq; + int ret = LIBLTE_ERROR_INVALID_INPUTS; + + if (q != NULL && + lte_cell_isvalid(&cell)) + { + ret = LIBLTE_ERROR; + + bzero(q, sizeof(refsignal_cs_t)); + bzero(&seq, sizeof(sequence_t)); + if (sequence_init(&seq, 2 * 2 * MAX_PRB)) { + goto free_and_exit; + } + + if (CP_ISNORM(cell.cp)) { + N_cp = 1; + } else { + N_cp = 0; + } + + q->cell = cell; + + for (i=0;ipilots[i] = vec_malloc(sizeof(cf_t) * REFSIGNAL_MAX_NUM_SF(q->cell.nof_prb)); + if (!q->pilots[i]) { + perror("malloc"); + goto free_and_exit; + } + } + + for (ns=0;nscell.nof_prb; i++) { + mp = i + MAX_PRB - cell.nof_prb; + /* save signal */ + q->pilots[ns/2][REFSIGNAL_PILOT_IDX(i,l,ns,q->cell)] = + (1 - 2 * (float) seq.c[2 * mp]) / sqrt(2) + + _Complex_I * (1 - 2 * (float) seq.c[2 * mp + 1]) / sqrt(2); + } + } + } + sequence_free(&seq); + ret = LIBLTE_SUCCESS; + } +free_and_exit: + if (ret == LIBLTE_ERROR) { + sequence_free(&seq); + refsignal_cs_free(q); + } + return ret; +} + +/** Deallocates a refsignal_cs_t object allocated with refsignal_cs_init */ +void refsignal_cs_free(refsignal_cs_t * q) +{ + int i; + + for (i=0;ipilots[i]) { + free(q->pilots[i]); + } + } + bzero(q, sizeof(refsignal_cs_t)); +} + + +inline uint32_t refsignal_fidx(lte_cell_t cell, uint32_t ns, uint32_t l, uint32_t port_id, uint32_t m) { + return 6*m + ((refsignal_cs_v(port_id, ns, lp(l,cell.cp)) + (cell.id % 6)) % 6); +} + +inline uint32_t refsignal_nsymbol(lte_cell_t cell, uint32_t ns, uint32_t l) { + return (ns%2)*CP_NSYMB(cell.cp)+lp(l,cell.cp); +} + +/* Maps a reference signal initialized with refsignal_cs_init() into an array of subframe symbols */ +int refsignal_cs_put_sf(lte_cell_t cell, uint32_t port_id, uint32_t sf_idx, + cf_t *pilots, cf_t *sf_symbols) +{ + uint32_t i, l, ns; + uint32_t fidx; + + if (lte_cell_isvalid(&cell) && + lte_sfidx_isvalid(sf_idx) && + lte_portid_isvalid(port_id) && + pilots != NULL && + sf_symbols != NULL) + { + + for (ns=2*sf_idx;ns<2*(sf_idx+1);ns++) { + for (l=0;l. + * This file is part of ALOE++ (http://flexnets.upc.edu/) + * + * ALOE++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ALOE++ is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with ALOE++. If not, see . + */ + +#include +#include "liblte/phy/phy.h" +#ifdef UNDEF_BOOL +#undef bool +#endif +#include "mex.h" + + +/** MEX function to be called from MATLAB to test the channel estimator + * + * [estChannel] = liblte_chest(cell_id, nof_ports, inputSignal, (optional) sf_idx) + * + * Returns a matrix of size equal to the inputSignal matrix with the channel estimates + * for each resource element in inputSignal. The inputSignal matrix is the received Grid + * of size nof_resource_elements x nof_ofdm_symbols_in_subframe. + * + * The sf_idx is the subframe index only used if inputSignal is 1 subframe length. + * + */ + +#define CELLID prhs[0] +#define PORTS prhs[1] +#define INPUT prhs[2] +#define NOF_INPUTS 3 +#define SFIDX prhs[3] + +void help() +{ + mexErrMsgTxt + ("[estChannel] = liblte_chest(cell_id, nof_ports, inputSignal,[sf_idx])\n\n" + " Returns a matrix of size equal to the inputSignal matrix with the channel estimates\n " + "for each resource element in inputSignal. The inputSignal matrix is the received Grid\n" + "of size nof_resource_elements x nof_ofdm_symbols.\n" + "The sf_idx is the subframe index only used if inputSignal is 1 subframe length.\n"); +} + +/* the gateway function */ +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + + int i; + lte_cell_t cell; + chest_dl_t chest; + cf_t *input_signal; + cf_t *ce[MAX_PORTS]; + double *outr0, *outi0, *outr1, *outi1; + + if (nrhs < NOF_INPUTS) { + help(); + return; + } + + if (!mxIsDouble(CELLID) && mxGetN(CELLID) != 1 && + !mxIsDouble(PORTS) && mxGetN(PORTS) != 1 && + mxGetM(CELLID) != 1) { + help(); + return; + } + + cell.id = (uint32_t) *((double*) mxGetPr(CELLID)); + cell.nof_prb = mxGetM(INPUT)/RE_X_RB; + cell.nof_ports = (uint32_t) *((double*) mxGetPr(PORTS)); + if ((mxGetN(INPUT)%14) == 0) { + cell.cp = CPNORM; + } else if ((mxGetN(INPUT)%12)!=0) { + cell.cp = CPEXT; + } else { + help(); + return; + } + + if (chest_dl_init(&chest, cell)) { + mexErrMsgTxt("Error initiating channel estimator\n"); + return; + } + + int nsubframes; + if (cell.cp == CPNORM) { + nsubframes = mxGetN(INPUT)/14; + } else { + nsubframes = mxGetN(INPUT)/12; + } + + uint32_t sf_idx=0; + if (nsubframes == 1) { + if (nrhs != NOF_INPUTS+1) { + help(); + return; + } + sf_idx = (uint32_t) *((double*) mxGetPr(SFIDX)); + } + + double *inr=(double *)mxGetPr(INPUT); + double *ini=(double *)mxGetPi(INPUT); + + /** Allocate input buffers */ + int nof_re = 2*CP_NSYMB(cell.cp)*cell.nof_prb*RE_X_RB; + for (i=0;i= 1) { + plhs[0] = mxCreateDoubleMatrix(1,nof_re * nsubframes, mxCOMPLEX); + outr0 = mxGetPr(plhs[0]); + outi0 = mxGetPi(plhs[0]); + } + if (nlhs == 2) { + plhs[1] = mxCreateDoubleMatrix(REFSIGNAL_MAX_NUM_SF(cell.nof_prb)*nsubframes, cell.nof_ports, mxCOMPLEX); + outr1 = mxGetPr(plhs[1]); + outi1 = mxGetPi(plhs[1]); + } + + for (int sf=0;sf= 1) { + for (i=0;iid < 504 && - cell->nof_ports > 0 && - cell->nof_ports < MAX_PORTS+1 && - cell->nof_prb > 5 && - cell->nof_prb < MAX_PRB+1 - ) { + +bool lte_cellid_isvalid(uint32_t cell_id) { + if (cell_id < 504) { return true; } else { return false; } } +bool lte_nofprb_isvalid(uint32_t nof_prb) { + if (nof_prb >= 6 && nof_prb <= MAX_PRB) { + return true; + } else { + return false; + } +} + +bool lte_cell_isvalid(lte_cell_t *cell) { + return lte_cellid_isvalid(cell->id) && + lte_portid_isvalid(cell->nof_ports) && + lte_nofprb_isvalid(cell->nof_prb); +} + void lte_cell_fprint(FILE *stream, lte_cell_t *cell) { fprintf(stream, "PCI: %d, CP: %s, PRB: %d, Ports: %d\n", cell->id, lte_cp_string(cell->cp), cell->nof_prb, cell->nof_ports); } +bool lte_sfidx_isvalid(uint32_t sf_idx) { + if (sf_idx <= NSUBFRAMES_X_FRAME) { + return true; + } else { + return false; + } +} + +bool lte_portid_isvalid(uint32_t port_id) { + if (port_id <= MAX_PORTS) { + return true; + } else { + return false; + } +} + bool lte_N_id_2_isvalid(uint32_t N_id_2) { if (N_id_2 < 3) { return true; diff --git a/lte/phy/lib/common/src/sequence.c b/lte/phy/lib/common/src/sequence.c index 2367a3e9b..4eac18a74 100644 --- a/lte/phy/lib/common/src/sequence.c +++ b/lte/phy/lib/common/src/sequence.c @@ -40,7 +40,7 @@ * It follows the 3GPP Release 8 (LTE) 36.211 * Section 7.2 */ -void generate_prs_c(sequence_t *q, uint32_t seed) { +void sequence_set_LTE_pr(sequence_t *q, uint32_t seed) { int n; uint32_t *x1, *x2; @@ -79,7 +79,7 @@ int sequence_LTE_pr(sequence_t *q, uint32_t len, uint32_t seed) { return LIBLTE_ERROR; } q->len = len; - generate_prs_c(q, seed); + sequence_set_LTE_pr(q, seed); return LIBLTE_SUCCESS; } @@ -92,6 +92,7 @@ int sequence_init(sequence_t *q, uint32_t len) { if (!q->c) { return LIBLTE_ERROR; } + q->len = len; } return LIBLTE_SUCCESS; } diff --git a/lte/phy/lib/filter/src/filter2d.c b/lte/phy/lib/filter/src/filter2d.c index 26e1af259..3d5428bda 100644 --- a/lte/phy/lib/filter/src/filter2d.c +++ b/lte/phy/lib/filter/src/filter2d.c @@ -32,6 +32,7 @@ #include "liblte/phy/utils/debug.h" +#include "liblte/phy/resampling/interp.h" #include "liblte/phy/filter/filter2d.h" #include "liblte/phy/utils/matrix.h" #include "liblte/phy/utils/vector.h" @@ -40,10 +41,10 @@ /* Useful macros */ #define intceil(X, Y) ((X-1)/Y+1) -#define idx(a, b) ((a)*(q->szfreq)+b) +#define idx(a, b) ((a)*(q->szfreq+q->nfreq)+b) -int filter2d_init(filter2d_t* q, float **taps, int ntime, int nfreq, int sztime, - int szfreq) { +int filter2d_init(filter2d_t* q, float **taps, uint32_t ntime, uint32_t nfreq, uint32_t sztime, + uint32_t szfreq) { int ret = -1; bzero(q, sizeof(filter2d_t)); @@ -54,18 +55,24 @@ int filter2d_init(filter2d_t* q, float **taps, int ntime, int nfreq, int sztime, matrix_copy((void**) q->taps, (void**) taps, ntime, nfreq, sizeof(float)); - q->output = vec_malloc((ntime+sztime)*(szfreq)*sizeof(cf_t)); + q->output = vec_malloc((ntime+sztime)*(szfreq+nfreq)*sizeof(cf_t)); if (!q->output) { goto free_and_exit; } - bzero(q->output, (ntime+sztime)*(szfreq)*sizeof(cf_t)); + bzero(q->output, (ntime+sztime)*(szfreq+nfreq)*sizeof(cf_t)); q->nfreq = nfreq; q->ntime = ntime; q->szfreq = szfreq; q->sztime = sztime; - + + q->norm = 0.0; + for (int i = 0; i < ntime; i++) { + for (int j = 0; j < nfreq; j++) { + q->norm += q->taps[i][j]; + } + } ret = 0; free_and_exit: if (ret == -1) { @@ -83,33 +90,26 @@ void filter2d_free(filter2d_t *q) { bzero(q, sizeof(filter2d_t)); } -int filter2d_init_default(filter2d_t* q, int ntime, int nfreq, int sztime, - int szfreq) { +int filter2d_init_ones(filter2d_t* q, uint32_t ntime, uint32_t nfreq, uint32_t sztime, + uint32_t szfreq) +{ int i, j; - int ret = -1; float **taps; + int ret = -1; if (matrix_init((void***) &taps, ntime, nfreq, sizeof(float))) { goto free_and_exit; } - /* Compute the default 2-D interpolation mesh */ for (i = 0; i < ntime; i++) { for (j = 0; j < nfreq; j++) { - if (j < nfreq / 2) - taps[i][j] = (j + 1.0) / (2.0 * intceil(nfreq, 2)); - - else if (j == nfreq / 2) - taps[i][j] = 0.5; - - else if (j > nfreq / 2) - taps[i][j] = (nfreq - j) / (2.0 * intceil(nfreq, 2)); + taps[i][j] = 1.0/(i+1); } } - INFO("Using default interpolation matrix of size %dx%d\n", ntime, nfreq); - if (verbose >= VERBOSE_DEBUG) { + INFO("Using all-ones interpolation matrix of size %dx%d\n", ntime, nfreq); + if (verbose >= VERBOSE_INFO) { matrix_fprintf_f(stdout, taps, ntime, nfreq); } @@ -126,30 +126,69 @@ free_and_exit: /* Moves the last ntime symbols to the start and clears the remaining of the output. * Should be called, for instance, before filtering each OFDM frame. */ -void filter2d_reset(filter2d_t *q) { +void filter2d_step(filter2d_t *q) { int i; for (i = 0; i < q->ntime; i++) { memcpy(&q->output[idx(i,0)], &q->output[idx(q->sztime + i,0)], - sizeof(cf_t) * (q->szfreq)); + sizeof(cf_t) * (q->szfreq+q->nfreq)); } + for (; i < q->ntime + q->sztime; i++) { - memset(&q->output[idx(i,0)], 0, sizeof(cf_t) * (q->szfreq)); + memset(&q->output[idx(i,0)], 0, sizeof(cf_t) * (q->szfreq+q->nfreq)); } } +void filter2d_reset(filter2d_t *q) { + bzero(q->output, (q->ntime+q->sztime)*(q->szfreq+q->nfreq)*sizeof(cf_t)); +} + /** Adds samples x to the from the given time/freq indexes to the filter * and computes the output. */ -void filter2d_add(filter2d_t *q, cf_t x, int time_idx, int freq_idx) { +void filter2d_add(filter2d_t *q, cf_t x, uint32_t time_idx, uint32_t freq_idx) { int i, j; - int ntime = q->ntime; - int nfreq = q->nfreq; + uint32_t ntime = q->ntime; + uint32_t nfreq = q->nfreq; + + if (freq_idx < q->szfreq && time_idx < q->sztime) { + DEBUG("Adding %f+%fi to %d:%d\n",__real__ x,__imag__ x,time_idx,freq_idx); + + for (i = 0; i < ntime; i++) { + for (j = 0; j < nfreq; j++) { + q->output[idx(i+time_idx, j+freq_idx)] += x * (cf_t)(q->taps[i][j])/q->norm; + } + } + } +} + +void filter2d_add_out(filter2d_t *q, cf_t x, int time_idx, int freq_idx) { + int i, j; + uint32_t ntime = q->ntime; + uint32_t nfreq = q->nfreq; + float norm=0; for (i = 0; i < ntime; i++) { for (j = 0; j < nfreq; j++) { - q->output[idx(i+time_idx, j+freq_idx - nfreq/2)] += x * (cf_t)(q->taps[i][j]); + if (i+time_idx >= 0 && j+freq_idx >= 0 && + i+time_idx < q->ntime+q->sztime && j+freq_idx < q->nfreq + q->szfreq) + { + norm += q->taps[i][j]; + } } - } + } + for (i = 0; i < ntime; i++) { + for (j = 0; j < nfreq; j++) { + if (i+time_idx >= 0 && j+freq_idx >= 0 && + i+time_idx < q->ntime+q->sztime && j+freq_idx < q->nfreq + q->szfreq) + { + q->output[idx(i+time_idx, j+freq_idx)] += x * (cf_t)(q->taps[i][j])/q->norm; + } + } + } +} + +void filter2d_get_symbol(filter2d_t *q, uint32_t nsymbol, cf_t *output) { + memcpy(output, &q->output[idx(nsymbol,q->nfreq/2)], sizeof(cf_t) * (q->szfreq)); } diff --git a/lte/phy/lib/mimo/src/precoding.c b/lte/phy/lib/mimo/src/precoding.c index 5599ea46f..6b2d52481 100644 --- a/lte/phy/lib/mimo/src/precoding.c +++ b/lte/phy/lib/mimo/src/precoding.c @@ -122,6 +122,9 @@ int precoding_type(cf_t *x[MAX_LAYERS], cf_t *y[MAX_PORTS], int nof_layers, return 0; } + +float y_mod[110*12*14]; + /* ZF detector */ int predecoding_single_zf(cf_t *y, cf_t *ce, cf_t *x, int nof_symbols) { for (int i=0;i 0 && M > 0) { ret = LIBLTE_ERROR; q->in_arg = vec_malloc(len * sizeof(float)); @@ -191,6 +191,28 @@ void interp_run(interp_t *q, cf_t *input, cf_t *output) { interp_run_offset(q, input, output, 0, 1); } +cf_t interp_linear_onesample(cf_t *input) { + float mag0=0, mag1=0, arg0=0, arg1=0, mag=0, arg=0; + mag0 = cabsf(input[0]); + mag1 = cabsf(input[1]); + arg0 = cargf(input[0]); + arg1 = cargf(input[1]); + mag = 2*mag1 -mag0; + arg = 2*arg1-arg0; + return mag * cexpf(I * arg); +} + +cf_t interp_linear_onesample2(cf_t *input) { + float re0=0, im0=0, re1=0, im1=0, re=0, im=0; + re0 = crealf(input[0]); + im0 = cimagf(input[1]); + re1 = crealf(input[0]); + im1 = cimagf(input[1]); + re = 2*re1-re0; + im = 2*im1-im0; + return (re+im*_Complex_I); +} + /* Performs 1st order linear interpolation with out-of-bound interpolation */ void interp_linear_offset(cf_t *input, cf_t *output, uint32_t M, uint32_t len, uint32_t off_st, uint32_t off_end) { uint32_t i, j; diff --git a/lte/phy/lib/utils/src/vector.c b/lte/phy/lib/utils/src/vector.c index 11b2c41ea..845743174 100644 --- a/lte/phy/lib/utils/src/vector.c +++ b/lte/phy/lib/utils/src/vector.c @@ -327,13 +327,21 @@ void vec_prod_conj_ccc(cf_t *x,cf_t *y, cf_t *z, uint32_t len) { #endif } -void vec_div_ccc(cf_t *x, cf_t *y, cf_t *z, uint32_t len) { +/* Complex division is conjugate multiplication + real division */ +void vec_div_ccc(cf_t *x, cf_t *y, float *y_mod, cf_t *z, uint32_t len) { + vec_prod_conj_ccc(x,y,z,len); + vec_abs_square_cf(y,y_mod,len); int i; for (i=0;igriddata for types) +% InterpWindow - Interpolation window type: ('Causal','Non-causal', +% 'Centred','Centered'). Note: 'Centred' and 'Centered' +% are equivalent +% InterpWinSize - Interpolation window size: +% 'Causal','Non-causal' - any number >=1 +% 'Centred','Centered' - odd numbers >=1 +% +% The 'TestEVM' pilot averaging will ignore other structure fields in +% CEC, and the method follows that described in TS36.104/TS36.141 Annex +% E/F for the purposes of transmitter EVM testing. +% +% The 'UserDefined' pilot averaging uses a rectangular kernel of size +% CEC.FreqWindow-by-CEC.TimeWindow and performs a 2D filtering operation +% upon the pilots. Note that pilots near the edge of the resource grid +% will be averaged less as they have no neighbors outside of the grid, +% or a limited number of neighbors outside of the grid obtained by the +% creation of virtual pilots. +% +% [HEST NOISEEST] = lteDLChannelEstimate(ENB, CHS, CEC, RXGRID) returns +% the estimated channel given cell-wide settings structure ENB, PDSCH +% transmission configuration CHS and channel estimator configuration CEC. +% CHS must include the following fields: +% TxScheme - Transmission scheme, one of: +% 'Port0' - Single-antenna port, Port 0 +% 'TxDiversity' - Transmit diversity scheme +% 'CDD' - Large delay CDD scheme +% 'SpatialMux' - Closed-loop spatial multiplexing scheme +% 'MultiUser' - Multi-user MIMO scheme +% 'Port5' - Single-antenna port, Port 5 +% 'Port7-8' - Single-antenna port, port 7 (when +% NLayers=1); Dual layer transmission, port 7 +% and 8 (when NLayers=2) +% 'Port8' - Single-antenna port, Port 8 +% 'Port7-14' - Up to 8 layer transmission, ports 7-14 +% PRBSet - A 1- or 2-column matrix, containing the 0-based Physical +% Resource Block indices (PRBs) corresponding to the resource +% allocations for this PDSCH. +% RNTI - Radio Network Temporary Identifier (16-bit) +% +% For the 'Port5', 'Port7-8', 'Port8' and 'Port7-14' transmission +% schemes, the channel estimator configuration structure CEC contains the +% following additional field: +% Reference - Optional. Specifies point of reference (signals to +% internally generate) for channel estimation. +% ('DMRS'(default), 'CSIRS') +% +% For the 'Port5', 'Port7-8', 'Port8' and 'Port7-14' transmission +% schemes, with Reference='DMRS', the channel estimation is performed +% using UE-specific reference signals and the returned channel estimate +% will be of size M-by-N-by-NRxAnts-by-NLayers. Alternatively, with +% Reference='CSIRS' the channel estimation is performed using the CSI +% reference signals (CSI) and the returned channel estimate will be of +% size M-by-N-by-NRxAnts-by-CSIRefP. For other transmission schemes the +% channel estimation is performed using cell-specific reference signals +% and the channel estimate will be of size M-by-N-by-NRxAnts-by-CellRefP. +% Note that CSI-RS based channel estimation and hence Reference='CSIRS' +% is strictly only valid within the standard for the 'Port7-14' +% transmission scheme. The optional CSIRSPeriod parameter controls the +% downlink subframes in which CSI-RS will be present, either always 'On' +% or 'Off', or defined by the scalar subframe configuration index Icsi-rs +% (0...154) or the explicit subframe periodicity and offset pair [Tcsi-rs +% Dcsi-rs] (TS 36.211 Section 6.10.5.3). +% +% For the 'Port7-8' and 'Port7-14' transmission schemes with +% 'UserDefined' pilot averaging, if CEC.TimeWindow = 2 or 4 and +% CEC.FreqWindow=1 the estimator will enter a special case where an +% averaging window of 2 or 4 pilots in time will be used to average the +% pilot estimates; the averaging is always applied across 2 or 4 pilots, +% regardless of their separation in OFDM symbols. This operation ensures +% that averaging is always done on 2 or 4 pilots. This provides the +% appropriate "despreading" operation required for the case of UE RS +% ports / CSI-RS ports which occupy the same time/frequency locations but +% use different orthogonal covers to allow them to be differentiated at +% the receiver. For the CSI-RS and any number of configured CSI-RS ports +% (given by ENB.CSIRefP), the pilot REs occur in pairs, one pair per +% subframe, that require averaging with CEC.TimeWindow=2 and will result +% in a single estimate per subframe. For the UE RS with between 1 and 4 +% layers (given by CHS.NLayers), the pilot REs occur in pairs, repeated +% in each slot, that require averaging with CEC.TimeWindow=2 and will +% result in two estimates per subframe, one for each slot; for between 5 +% and 8 layers, the pairs are distinct between the slots of the subframe +% and the required averaging is CEC.TimeWindow=4, resulting in one +% estimate per subframe. +% +% Example: +% Transmit RMC R.12 (4-antenna transmit diversity), model the propagation +% channel by combining all transmit antennas onto one receive antenna, +% OFDM demodulate and finally channel estimate. +% +% enb = lteRMCDL('R.12'); +% cec = struct('FreqWindow',1,'TimeWindow',1,'InterpType','cubic',... +% 'PilotAverage','UserDefined','InterpWinSize',3,... +% 'InterpWindow','Causal'); +% txWaveform = lteRMCDLTool(enb,[1;0;0;1]); +% rxWaveform = sum(txWaveform,2); +% rxGrid = lteOFDMDemodulate(enb,rxWaveform); +% hest = lteDLChannelEstimate(enb,cec,rxGrid); +% +% See also lteOFDMDemodulate, lteEqualizeMMSE, lteEqualizeZF, +% lteDLPerfectChannelEstimate, griddata. + +% Copyright 2009-2013 The MathWorks, Inc. + +function [H_EST, NoisePowerEst, AvgEstimates, Estimates] = lteDLChannelEstimate2(varargin) + + if(isstruct(varargin{2})) + if(isstruct(varargin{3})) + PDSCH = varargin{2}; + CEC = varargin{3}; + RXGRID = varargin{4}; + else + PDSCH = []; + CEC = varargin{2}; + RXGRID = varargin{3}; + end + else + % If no configuration structure then use the TestEVM method + PDSCH = []; + CEC.PilotAverage = 'TestEVM'; + RXGRID = varargin{2}; + end + + ENB = varargin{1}; + % Get dimensions of resource grid + Dims = lteDLResourceGridSize(ENB); + K = Dims(1); + Lsf = Dims(2); + + % Determine number of subframes + nsfs = size(RXGRID,2)/Dims(2); + + % Determine number of Tx- and RxAntennas + NRxAnts = size(RXGRID,3); + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port5' 'Port7-8' 'Port8' 'Port7-14'}))) + if (~isfield(CEC,'Reference')) + CEC.Reference='DMRS'; + end + if (strcmpi(CEC.Reference,'DMRS')==1) + NTx = PDSCH.NLayers; + elseif (strcmpi(CEC.Reference,'CSIRS')==1) + NTx = ENB.CSIRefP; + else + error('lte:error','Reference must be "DMRS" or "CSIRS", see help for details.'); + end + else + if (~isfield(CEC,'Reference')) + CEC.Reference='CellRS'; + end + NTx = ENB.CellRefP; + end + + % Initialize size of estimated channel grid + H_EST = zeros([size(RXGRID,1) size(RXGRID,2) size(RXGRID,3) NTx]); + + % Preallocate noise power estimate vector for speed + noiseVec = zeros(size(NRxAnts,NTx)); + + if (NTx == 4) + nn=3; + else + nn=NTx; + end + RefEstimates=zeros(nn,2*2*2*ENB.NDLRB*10); + + for rxANT = 1:NRxAnts + for nTxAntLayer = 1:NTx + + % Extract pilot symbols from received grid and calculate least + % squares estimate + [ls_estimates,specialvec,DwPTS]= GetPilotEstimates(ENB,nTxAntLayer-1,RXGRID(:,:,rxANT),PDSCH,CEC); + + Estimates(nTxAntLayer,1:length(ls_estimates(3,:))) = ls_estimates(3,:); + + % Average the pilots as defined by CEC.PilotAverage, this can + % be 'UserDefined' or 'TestEVM'. + % Note: Setting the window size to 1x1 is equivalent to no + % averaging. It is recommended that no averaging should be + % carried out when a high SNR is present, as this would have an + % adverse effect on the least squares estimates + if strcmpi(CEC.PilotAverage,'UserDefined')&&(CEC.FreqWindow==1)&&(CEC.TimeWindow==1) + P_EST = ls_estimates; + ScalingVec = 0; + else + [P_EST, ScalingVec]= PilotAverage(PDSCH,CEC,H_EST,ls_estimates); + end + + AvgEstimates(nTxAntLayer,1:length(P_EST(3,:))) = P_EST(3,:); + + if strcmpi(CEC.PilotAverage,'TestEVM') + % Channel estimation as defined in TS36.141 Annex F.3.4 is + % performed, the pilots are averaged in time and frequency + % and the resulting vector is linearly interpolated. + %---------------------- + % Interpolate eq coefficients extrapolating to account for + % missing pilots at the edges + interpEqCoeff = interp1(find(P_EST~=0),P_EST(P_EST~=0),(1:length(P_EST)).','linear','extrap'); + % --------------------- + % The DC carrier needs to be accounted for during + % interpolation + % Get pilot symbol index values + ind = lteCellRSIndices(ENB,0); + % Using initial index value calculate position of pilot + % symbols within interpolated coefficients + if ind(1)-3>0 + vec = (ind(1)-3):3:K; + else + vec = ind(1):3:K; + end + % Determine the index values of the pilots either side of + % the DC carrier + preDC = vec(length(vec)/2); + postDC = vec(1+length(vec)/2); + % Interpolate the section using the DC carrier and insert + % this into the interpolated vector containing the + % coefficients + interpEqCoeffTemp = interp1(1:4:5,[interpEqCoeff(preDC) interpEqCoeff(postDC)],1:5,'linear'); + interpEqCoeff(preDC:K/2) = interpEqCoeffTemp(1:length(interpEqCoeff(preDC:K/2))); + interpEqCoeff(1+K/2:postDC) = interpEqCoeffTemp(2+length(interpEqCoeff(preDC:K/2)):end); + %---------------------- + % Generate grid by replicating estimated equalizer channel + % coefficients + H_EST(:,:,rxANT,nTxAntLayer) = repmat(interpEqCoeff,1,Lsf*nsfs); + % The value of the noise averaged pilot symbols are placed + % into P_EST(3,:) matrix and these are used to determine + % the noise power present on the pilot symbol estimates of + % the channel. + P_EST = ls_estimates; + tempGrid = squeeze(H_EST(:,:,rxANT,nTxAntLayer)); + P_EST(3,:) = tempGrid(sub2ind(size(H_EST),P_EST(1,:),P_EST(2,:))); + + else + % Channel estimation is performed by analyzing up to + % 'CEC.InterpWinSize' subframes together and estimating + % virtual pilots outwith the bounds of the subframe. + % Averaging window can be set to Causal, Non-causal or + % Centred, Centered. + + % Using the desired subframeAveraging type the initial + % position of the window is determined + if(CEC.InterpWinSize>=1) + if strcmpi(CEC.InterpWindow,'Centred')||strcmpi(CEC.InterpWindow,'Centered') + % For Centered averaging the window size must be odd + if ~mod(CEC.InterpWinSize,2) + error('lte:error','Window size must be odd for centered window type'); + end + % For Centered windowing both current,future and past + % data are used to estimate current channel coefficients + x = floor(CEC.InterpWinSize/2); + y = floor(CEC.InterpWinSize/2); + elseif strcmpi(CEC.InterpWindow,'Causal') + % For causal windowing current and past data are used + % to estimate current subframe channel coefficients + x = 0; + y = CEC.InterpWinSize-1; + elseif strcmpi(CEC.InterpWindow,'Non-causal') + % For non-causal windowing current and future data are + % used to estimate current subframe channel + % coefficients + x = CEC.InterpWinSize-1; + y = 0; + else + error('lte:error','Channel estimation structure field InterpWindow must be one of the following: Causal,Non-causal,Centred,Centered'); + end + else + error('lte:error','InterpWinSize cannot be less than 1'); + end + + % Interpolate using averaged pilot estimates and defined + % interpolation settings + for sf = 0:nsfs-1 + if specialvec(sf+1)=='U' + H_EST(:,sf*Lsf+1:(sf+1)*Lsf,rxANT,nTxAntLayer) = NaN; + else + % Extract the pilots from required subframes + p_use = P_EST; + p_use(:,p_use(4,:)>(sf+1)+x)=[]; + p_use(:,p_use(4,:)<(sf+1)-y)=[]; + + % Account for DC offset + p_use(1,(p_use(1,:)>K/2)) = p_use(1,(p_use(1,:)>K/2))+1; + + if strcmpi(CEC.InterpType,'Cubic')||strcmpi(CEC.InterpType,'Linear') + % VPVEC is used to determine if virtual pilots are + % needed at the beginning or end of the + % interpolation window, if the current subframe,sf, + % is located at beginning or end of the window then + % virtual pilots are created accordingly + vpvec = unique(p_use(4,:)); + + if (strcmp(CEC.Reference,'CellRS')==1) + % Create virtual pilots and append to vector + % containing pilots using cell-specific RS + % methodology. + vps = createVirtualPilots(ENB,p_use(1:3,:),1,1,(sf+1==min(vpvec)),(sf+1==max(vpvec))); + p_use = [p_use(1:3,:) vps]; + else + % Create edge virtual pilots suitable for + % UE RS which can have partial bandwidth, + % or CSI-RS. + if (~isempty(p_use)) + vps = createEdgeVirtualPilots(ENB,p_use(1:3,:)); + p_use = [p_use(1:3,:) vps]; + end + end + end + % Perform 2D interpolation + % Interpolation is carried out on a (K+1)-by-L + % matrix to account for DC offset being added in + if (~isempty(p_use)) + Htemp = griddata(p_use(2,:)-Lsf*sf,p_use(1,:),p_use(3,:),1:Lsf,(1:K+1)',CEC.InterpType); %#ok + % Remove DC offset + Htemp(1+(K/2),:) = []; + + if specialvec(sf+1)=='S' + Htemp(:,DwPTS+1:Lsf) = NaN; + end + H_EST(:,sf*Lsf+1:(sf+1)*Lsf,rxANT,nTxAntLayer) = Htemp; + + if isnan(H_EST) + error('lte:error','H_EST NaN'); + end + end + end + end + end + if nTxAntLayer<3 || (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port7-8' 'Port7-14'}))) + % The noise level present can be determined using the noisy + % least squares estimates of the channel at pilot symbol + % locations and the noise averaged pilot symbol estimates of the + % channel + noise = ls_estimates(3,~isnan(ls_estimates(3,:))) - P_EST(3,~isnan(P_EST(3,:))); + if strcmpi(CEC.PilotAverage,'UserDefined') + noise = sqrt(ScalingVec./(ScalingVec+1)).*noise; + end + + % Additional averaging for noise estimation in LTE-A case, + % to suppress interference from orthogonal sequences on + % other antennas in same time-frequency locations. + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port7-8' 'Port7-14'})) &&... + (CEC.TimeWindow==2 || CEC.TimeWindow==4) && CEC.FreqWindow==1) + if (~isempty(ls_estimates)) + if (strcmpi(CEC.Reference,'DMRS')==1) + temp=[]; + pilotSC = unique(ls_estimates(1,:)); + for i = pilotSC + x=find(ls_estimates(1,:)==i); + x=[x(1) x(end)]; + temp=[temp mean(noise(x))*2]; %#ok + end + noise=temp; + end + end + end + + % Taking the variance of the noise present on the pilot symbols + % results in a value of the noise power for each transmit and + % receive antenna pair + if (isempty(noise)) + noiseVec(rxANT,nTxAntLayer)=NaN; + else + noiseVec(rxANT,nTxAntLayer) = mean(noise.*conj(noise)); + end + end + end + end + + % The mean of the noise power across all the transmit/receive antenna + % pairs is used as the estimate of the noise power + NoisePowerEst = mean(mean(noiseVec)); + +end + +% GetPilotEstimates Obtain the least squares estimates of the reference +% signals +% [ls_estimates] = GetPilotEstimates(ENB,nTxAntLayer,RXGRID) Extracts the +% reference signals for a specific transmit/receive antenna pair and +% calculates their least squares estimate. The results are placed in a +% 3xNp matrix containing the subcarrier and OFDM symbol location, row and +% column subscripts, and value. Np is the number of cell specific +% reference (pilot) symbols per resource grid +% +% RXGRID is an MxN array where M is the number of subcarriers, +% N is the number of OFDM symbols. Dimension M must be +% 12*NDLRB where NDLRB must be {6,15,25,50,75,100}. +% Dimension N must be a multiple of number of symbols in a subframe L, +% where L=14 for normal cyclic prefix and L=12 for extended cyclic +% prefix. +% +% nTxAntLayer defines which transmit antennas' or layers' pilot symbols to extract +% +% ENB is a structure and must contain the following fields: +% NDLRB - Number of downlink resource blocks +% NCellID - Physical layer cell identity +% CellRefP - Number of transmit antenna ports {1,2,4} +% CyclicPrefix - Optional. Cyclic prefix length{'Normal'(default),'Extended'} +% +% +% Example +% Return the least squares estimate of the pilots symbols in a 3xNp array +% where the first row is the subcarrier index,k the second row is the OFDM +% symbol,l and the third row defines the least squares estimate of the +% pilot located at that position. +% +% enb=struct('NDLRB',6,'CellRefP',1,'NCellID',0,'CyclicPrefix','Normal'); +% rxGrid=ones(lteDLResourceGridSize(enb)); +% nTxAntLayer=0; +% ls_estimates = GetPilotEstimates(enb,nTxAntLayer,rxGrid).' +% ans = +% 1.0000 1.0000 -0.7071 - 0.7071i +% 7.0000 1.0000 -0.7071 - 0.7071i +% 13.0000 1.0000 -0.7071 - 0.7071i +% . . . +% . . . +% . . . +% 58.0000 12.0000 -0.7071 + 0.7071i +% 64.0000 12.0000 0.7071 - 0.7071i +% 70.0000 12.0000 0.7071 - 0.7071i + +% Copyright 2009-2010 The MathWorks, Inc. +function [ls_estimates,specialvec,DwPTS] = GetPilotEstimates(ENB,nTxAntLayer,RXGRID,PDSCH,CEC) + +% Get dimensions of resource grid +Dims = lteDLResourceGridSize(ENB); +nsfs = size(RXGRID,2)/Dims(2); +K = Dims(1); +L = Dims(2); + +if (rem(nsfs,1)~=0) + error('lte:error','The received grid input must contain a whole number of subframes.'); +end + +if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port5' 'Port7-8' 'Port8' 'Port7-14'}))) + if (strcmp(CEC.Reference,'DMRS')==1) + % create UE RS indices + PDSCH.NTxAnts=0; + PDSCH.W=[]; + linearIndSF = lteDMRSIndices(ENB,PDSCH,'mat'); + else + % create CSI-RS indices and symbols and deal with zero removal + linearIndSF = lteCSIRSIndices(ENB,'mat'); + CsiRS = lteCSIRS(ENB,'mat'); + CsiRS = CsiRS(:,nTxAntLayer+1); + linearIndSF(CsiRS==0,:)=[]; + CsiRS(CsiRS==0)=[]; + end + % If more than one TxAntenna is used we need to adjust indices values so that + % they start in first antenna plane + linearIndSF = linearIndSF(:,nTxAntLayer+1) - (K*L*nTxAntLayer); +else + % Extract linear indices of reference signals for particular TxAntenna + dumENB = ENB; + dumENB.NSubframe = 0; + linearIndSF = lteCellRSIndices(dumENB,nTxAntLayer); + % If more than one TxAntenna is used we need to adjust indices values so that + % they start in first antenna plane + linearIndSF = linearIndSF - (K*L*nTxAntLayer); +end + +% Create grid for pilot symbols +linearIndGrid = zeros(size(linearIndSF,1),nsfs); +idealPilotGrid = zeros(size(linearIndSF,1),nsfs); +offset = double(ENB.NSubframe); +specialvec = char(zeros(1,0)); +DwPTS = 0; +for sf = ENB.NSubframe:ENB.NSubframe+nsfs-1 + ENB.NSubframe = mod(sf,10); + linearIndGrid(:,sf-offset+1) = (linearIndSF+(sf*K*L)); + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port5' 'Port7-8' 'Port8' 'Port7-14'}))) + if (strcmp(CEC.Reference,'DMRS')==1) + % create UE RS symbols + DMRS = lteDMRS(ENB,PDSCH,'mat'); + idealPilotGrid(1:length(DMRS),sf-offset+1) = DMRS(:,nTxAntLayer+1); + else + % use CSI-RS symbols created earlier + idealPilotGrid(1:length(CsiRS),sf-offset+1) = CsiRS; + end + else + CellRS = lteCellRS(ENB, nTxAntLayer); + idealPilotGrid(1:length(CellRS),sf-offset+1) = CellRS; + end + dupinfo = lteDuplexingInfo(ENB); + if strcmpi(dupinfo.SubframeType,'Special') + specialvec = [specialvec 'S']; %#ok + DwPTS = dupinfo.NSymbolsDL; + elseif strcmpi(dupinfo.SubframeType,'Uplink') + specialvec = [specialvec 'U']; %#ok + elseif strcmpi(dupinfo.SubframeType,'Downlink') + specialvec = [specialvec 'D']; %#ok + else + error('lte:error','Invalid SubframeType, in structure dupinfo'); + end +end + +if (~isempty(idealPilotGrid)) + ulvec = find(idealPilotGrid(size(idealPilotGrid,1),:)==0); + linearIndGrid(:,ulvec) = []; + idealPilotGrid(:,ulvec) = []; +else + ulvec = []; +end + +linearIndGrid = linearIndGrid - (offset*K*L); +% Extract the row and column subscripts of the pilot symbols for entire +% grid +[p_estSC, p_estSym] = ind2sub(size(RXGRID),linearIndGrid); + +% Calculate least squares channel estimates at pilot locations +p_est = RXGRID(linearIndGrid)./idealPilotGrid; + +% Create a vector - [k;l;p_est] +sfref = repmat((1:nsfs),length(linearIndSF),1); +sfref(:,ulvec) = []; +ls_estimates = [double(p_estSC(:).') ; double(p_estSym(:).') ; p_est(:).';double(sfref(:)).']; + +end + +%PilotAverage Average reference signals +% [P_EST] = PilotAverage(PDSCH,CEC,H_EST,LS_EST) performs a moving average of pilot +% symbols +% +% LS_EST is a 3xNp matrix containing the least square estimates of the +% pilots symbols and their column and row indices within the received +% grid. +% LS_EST = [k;l;p_est] +% H_EST is an MxN matrix and defines the size of the grid that the +% averaging will be performed on +% +% CEC is a structure which defines the type of channel estimation +% performed. CEC must contain a set of the following fields: +% PilotAverage - Type of pilot averaging {'TestEVM', 'UserDefined'} +% FreqWindow - Size of window used to average in frequency in +% resource elements. +% TimeWindow - Size of window used to average in time in resource +% elements +% InterpType - Type of 2D interpolation used(see griddata for types) +% +% The dimensions of the averaging window are defined in structure CEC. +% The window is defined in terms of Resource Elements, and depending on +% the size of the averaging window, averaging will be performed in either +% the time or frequency direction only, or a combination of both creating +% a square/rectangular window. The pilot to be averaged will always be +% placed at the center of the window, therefore the window size must be +% an odd number. +% +% Frequency Direction: 9x1 Time Direction: 1x9 +% +% x +% x +% x +% x +% P x x x x P x x x x +% x +% x +% x +% x +% +% +% Square window: 9x9 +% x x x x x x x x x +% P x x x x x x P x +% x x x x x x x x x +% x x x x x x x x x +% x x x x P x x x x +% x x x x x x x x x +% x x x x x x x x x +% P x x x x x x P x +% x x x x x x x x x +% +% Performing EVM compliance testing as per TS36.141 AnnexF.3.4, requires +% time averaging be done over 10 subframes(i.e. 1 Frame) across each +% pilot symbol carrying subcarrier creating a TotalNumberPilotsx1 vector. +% This is then frequency averaged using a moving window average with a window +% size of 19. This type of averaging can be performed by setting the +% PilotAverage type in structure CEC to 'TestEVM'. + +% Copyright 2009-2013 The MathWorks, Inc. + +function [P_EST, scalingVec] = PilotAverage(PDSCH,CEC,H_EST,P_EST) + + switch lower(CEC.PilotAverage) + + case 'userdefined' + if(CEC.FreqWindow<1)||(CEC.TimeWindow<1) + error('lte:error','Frequency and time averaging window size cannot be less than 1'); + end + + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port7-8' 'Port7-14'})) &&... + (CEC.TimeWindow==2 || CEC.TimeWindow==4) && CEC.FreqWindow==1) + + % Perform averaging in time direction using a moving window + % of size N = CEC.TimeWindow + N = CEC.TimeWindow; + + if (~isempty(P_EST)) + % Average only subcarriers which contain RS symbols. Extract RS + % symbols on a per subcarrier basis and use window to average. + for subcarrier = unique(P_EST(1,:)) + + % Store DRS symbols from relevant slot in vector for easy access + symbVec = (P_EST(3,P_EST(1,:)==subcarrier)).'; + + % Define temporary vector to store the averaged values + avgVec = zeros(size(symbVec)); + + % Check length of input at least equal to size of window; + % if not, reduce window size. Despreading will be incomplete + % here + if (length(symbVec) < N) + N = length(symbVec); + end + + % Determines position of window w.r.t element being averaged, + % and performs the averaging + for symbNo = 1:length(symbVec) + if (symbNo-(N/2+1)<=0) + avgVec(symbNo)= sum(symbVec(1:N))/N; + else + avgVec(symbNo) = sum(symbVec(end-(N-1):end))/N; + end + end + + % Update P_EST with averaged values + P_EST(3,P_EST(1,:)==subcarrier)= avgVec; + end + + end + + % This vector is used to scale the noise by the number of averaging + % elements in the window + scalingVec = ones(size(P_EST(3,:)))*N; + + else + + if (strcmpi(CEC.InterpWindow,'Centred')||strcmpi(CEC.InterpWindow,'Centered')) + if (~mod(CEC.FreqWindow,2)||~mod(CEC.TimeWindow,2)) + error('lte:error','Window size must be odd in time and frequency for centred/centered window type'); + end + end + + % Define number of subcarriers + K = size(H_EST,1); + + % Define an empty resource grid with the DC offset subcarrier + % inserted + grid = zeros(size(H_EST,1)+1,size(H_EST,2)); + + % Account for DC offset + P_EST(1,(P_EST(1,:)>K/2)) = P_EST(1,(P_EST(1,:)>K/2))+1; + + % Place the pilot symbols back into the received grid with the + % DC offset subcarrier in place + grid(sub2ind(size(grid),P_EST(1,:),P_EST(2,:))) = P_EST(3,:); + + % Define convolution window + kernel = ones(CEC.FreqWindow,CEC.TimeWindow); + + % Perform convolution + grid2=grid; + grid = conv2(grid,kernel,'same'); + + % Extract only pilot symbol location values and set the rest of + % the grid to zero + tempGrid = zeros(size(grid)); + tempGrid(sub2ind(size(grid),P_EST(1,:),P_EST(2,:))) = grid(sub2ind(size(grid),P_EST(1,:),P_EST(2,:))); + grid = tempGrid; + + % Normalize pilot symbol values after convolution + [grid, scalingVec] = normalisePilotAverage(CEC,P_EST,grid); + + % Place averaged values back into pilot symbol matrix + P_EST(3,:) = grid(sub2ind(size(tempGrid),P_EST(1,:),P_EST(2,:))); + + % Remove DC offset + P_EST(1,(P_EST(1,:)>K/2)) = P_EST(1,(P_EST(1,:)>K/2))-1; + + end + + + case 'testevm' + % As defined in TS36.141 Annex F.3.4 pilots are averaged in + % time across all pilot carrying subcarriers. The resulting + % vector is averaged in frequency direction with a moving + % averaging window of 19. + + % Get phase (theta) and magnitude(radius) of complex estimates + [theta,radius] = cart2pol(real(P_EST(3,:)),imag(P_EST(3,:))); + + % Declare vector to temporarily store averaged phase and + % magnitude + phasemagVec = zeros(size(H_EST,1),2); + + % Perform averaging in time and frequency for phase and + % magnitude separately and recombine at end + for phaseOrMag = 1:2 + % Define temporary P_EST vector and set estimates row to + % phase then magnitude + P_EST_temp = P_EST; + if phaseOrMag==1 + P_EST_temp(3,:) = theta; + elseif phaseOrMag==2 + P_EST_temp (3,:) = radius; + end + + % Declare averaging vector + avgVec = zeros(size(H_EST,1),1); + vec = []; + + % Average across pilot symbol carrying subcarriers in time, + % unwrapping performed on a subcarrier basis + for i = 1:size(H_EST,1) + if ~isempty(P_EST(3,(P_EST(1,:)==i))) + if phaseOrMag==1 + avgVec(i) = mean(unwrap(P_EST_temp(3,(P_EST(1,:)==i)))); + + elseif phaseOrMag==2 + avgVec(i) = mean(P_EST_temp(3,(P_EST(1,:)==i))); + end + vec = [vec i]; %#ok + end + end + + % Remove subcarriers which contained no pilot symbols, then + % perform frequency averaging + avgVec = avgVec(vec); + + % Perform averaging in frequency direction using a moving + % window of size N, N = 19 in TS36.141 Annex F.3.4 + % Performs a moving average of window size N. At the edges, where less than + % N samples are span the window size is reduce to span 1, 3, 5, 7 ... + % samples. + N = 19; + + % avgVec must be a column vector + if (size(avgVec,2) > 1) + error('lte:error','Input to moving average function must be a column vector.') + end + + % N max window size 19 as defined in Annex F.3.4 + if ~mod(N,2) + error('lte:error','Window size N must be odd'); + end + + if (length(avgVec) < N) + % sprintf ('Input signal must have at least %d elements',N); + error('lte:error','Input signal must have at least %d elements',N); + end + + % Use filter to perform part of the averaging (not normalized) + zeroPad = zeros(N-1,1); + data = [avgVec; zeroPad]; + weights = ones(N,1); + freqAvg = filter(weights,1,data); + + % Remove unwanted elements + removeIdx = [ 2:2:N ( length(data)-(2:2:N)+1 )]; + freqAvg(removeIdx) = []; + + % Normalization factor + normFactor = [1:2:N-2 N*ones(1,length(freqAvg)-(N-1)) N-2:-2:1]'; + freqAvg = freqAvg./normFactor; + + % Place frequency averaged pilots into temporary storage + % vector + phasemagVec(vec,phaseOrMag) = freqAvg; + + end + % Convert averaged symbols back to Cartesian coordinates + [X,Y] = pol2cart(phasemagVec(:,1),phasemagVec(:,2)); + P_EST = complex(X,Y); + scalingVec = []; + + otherwise + error('lte:error','Channel estimation structure field PilotAverage must be one of the following: "UserDefined" or "TestEVM"'); + + end +end + +function [avgGrid, scalingVec] = normalisePilotAverage(CEC,p_est,grid) + + % Determine total number of pilots within half a subframe + nPilots = length(p_est); + + avgGrid = zeros(size(grid)); + scalingVec = zeros(1,size(p_est,2)); + + for n = 1:nPilots + % Determine in which subcarrier and OFDM symbol pilot is located + sc = p_est(1,n); + sym = p_est(2,n); + % Determine number of REs to look at either side of pilot + % symbol in time and frequency + half_freq_window = floor(CEC.FreqWindow/2); + half_time_window = floor(CEC.TimeWindow/2); + % Set the location of the window at the back of the pilot + % to be averaged in frequency direction + upperSC = sc-half_freq_window; + % If this location is outwith the grid dimensions set it to + % lowest subcarrier value + if upperSC<1 + upperSC = 1; + end + % Set the location of the window in front of the pilot to + % be averaged in frequency direction + lowerSC = sc+half_freq_window; + % If this location is outwith the grid dimensions set it to + % highest subcarrier value + if lowerSC>size(grid,1) + lowerSC = size(grid,1); + end + % Set the location of the window in front of the pilot to + % be averaged in time direction + leftSYM = sym-half_time_window; + % If this location is outwith the grid dimensions set it to + % lowest OFDM symbol value + if leftSYM<1 + leftSYM = 1; + end + % Set the location of the window at the back of the pilot + % to be averaged in time direction + rightSYM = sym+half_time_window; + % If this location is outwith the grid dimensions set it to + % highest OFDM symbol value + if rightSYM>size(grid,2) + rightSYM = size(grid,2); + end + + % Define the window to average using the determined + % subcarrier and OFDM symbol values + avgVec = grid(upperSC:lowerSC,leftSYM:rightSYM); + + % Remove the zero values from the window so that the average is + % calculated using only valid pilots + avgVec = avgVec(avgVec~=0); + % Average the desired pilot using all pilots within + % averaging window + if isempty(avgVec) + avgVec = 0; + end + + avgGrid(sc,sym) = grid(sc,sym)/numel(avgVec); + scalingVec(n) = numel(avgVec); + end +end + +%createVirtualPilots Create virtual pilots. +% [VP] = createVirtualPilots(..) returns a 3xNvp matrix which contains +% the time/frequency and least squares estimate of virtual pilots. +% Nvp is the number of virtual pilots. +% +% The inputs are: +% +% ENB structure which must contain the following fields: +% NDLRB - Number of downlink resource blocks +% CyclicPrefix - Optional. Cyclic prefix length {'Normal'(default),'Extended'} +% +% P_EST - A matrix containing the pilots within the grid. The matrix must +% have three rows. Each column contains the information for each pilot: +% [subcarrier_indices symbol_indices pilot_channel_estimate] +% +% ANTENNA - The current antenna port +% CREATE_TOP,CREATE_BOTTOM,CREATE_FRONT,CREATE_END - a value {0 or 1} +% indicating whether to return virtual pilots in each position. A 1 +% indicates virtual pilots should be returned, a 0 indicates they should +% not be returned. +% +% Example: +% +% Create virtual pilots on the top of a grid of pilots (estimates) +% VP_TOP = createVirtualPilots(ENB,estimates,0,1,0,0,0); +% +% Add virtual pilots to current pilot estimates for interpolation +% estimates = [estimates VP_TOP]; + +% Copyright 2009-2010 The MathWorks, Inc. +function VP = createVirtualPilots(ENB,P_EST,CREATE_TOP,CREATE_BOTTOM,CREATE_FRONT,CREATE_END) + + % Initialize outputs + VP = zeros(3,0); + + % Get resource grid dimensions + dims = lteDLResourceGridSize(ENB); + K = dims(1); + + % Virtual pilot cut-off values - indices in time or frequency at which + % the virtual pilot shall be discarded. As the symbol indices of the + % pilots passed in may be positive and negative the start and end of + % the effective grid varies therefore affecting the front and end pilot + % cut-off values + coTop = -6; % + coBottom = K+6; % + coFront = min(P_EST(2,:))-7; % floor(min(P_EST(2,:))/L)*L-5; % + coEnd = max(P_EST(2,:))+7; % ceil(max(P_EST(2,:))/L)*L+5; + + % Number of pilots in a subcarrier is dependent upon antenna +% noPilotinSC = (2-floor(ANTENNA/2))*nSF; + noPilotinSC = numel(P_EST(P_EST(1,:)==P_EST(1,1))); + + % Number of pilots in a symbol + noPilotinSym = K/6; + + % Sort by subcarrier + [~, indices] = sortrows(P_EST(1,:).'); + p_est_sorted_sc = P_EST(:,indices).'; + + if CREATE_TOP + % Repeat first and second SC containing pilots + p_est_rep_above = p_est_sorted_sc(1+noPilotinSC:noPilotinSC*2,:); + p_est_rep_above = [p_est_rep_above.' p_est_sorted_sc(1:noPilotinSC,:).'].'; + + % Subtract six to normalize subcarrier indices + p_est_rep_above(:,1) = p_est_rep_above(:,1) - 6; + + % Remove virtual pilots which are too far + p_est_rep_above(p_est_rep_above(:,1) < coTop,:) = []; + + % Initialize vp vector + VP_TOP = zeros(size(p_est_rep_above)); + + % Create virtual pilots on top + for p = 1:size(p_est_rep_above) + VP_TOP(p,:) = calculatePilot(p_est_rep_above(p,:), p_est_sorted_sc); + end + + % Transpose to make suitable for output + VP = [VP VP_TOP.']; + end + + if CREATE_BOTTOM + % Repeat last and second last SC containing pilots + p_est_rep_below = p_est_sorted_sc(end-noPilotinSC*2+1:end-noPilotinSC,:); + p_est_rep_below = [p_est_rep_below.' p_est_sorted_sc(end-noPilotinSC+1:end,:).'].'; + + % Add six to normalize subcarrier indices + p_est_rep_below(:,1) = p_est_rep_below(:,1) + 6; + + % Remove virtual pilots which are too far + p_est_rep_below(p_est_rep_below(:,1) > coBottom,:) = []; + + % Initialize vp vector + VP_BOTTOM = zeros(size(p_est_rep_below)); + + % Create virtual pilots on bottom + for p = 1:size(p_est_rep_below) + VP_BOTTOM(p,:) = calculatePilot(p_est_rep_below(p,:), p_est_sorted_sc); + end + + % Transpose to make suitable for output + VP = [VP VP_BOTTOM.']; + end + + % Sort by symbol + [~, indices] = sortrows(P_EST(2,:).'); + p_est_sorted_sym = P_EST(:,indices).'; + + if CREATE_FRONT + % Repeat first symbol containing pilots + p_est_rep_front = p_est_sorted_sym(noPilotinSym+1:(noPilotinSym*2),:); + + % Add subcarrier above and below + p_est_rep_front = addSCs(p_est_rep_front,noPilotinSym); + + % Repeat second symbol containing pilots + p_est_rep_front_second = p_est_sorted_sym(1:noPilotinSym,:); + + % Add subcarrier above and below + p_est_rep_front_second = addSCs(p_est_rep_front_second,noPilotinSym); + + % Concatenate to create virtual pilots + p_est_rep_front = [p_est_rep_front.' p_est_rep_front_second.'].'; + + % Subtract seven to normalize symbol indices + p_est_rep_front(:,2) = p_est_rep_front(:,2) - 7; + + % Remove virtual pilots which are too far + p_est_rep_front = removeTooFar(p_est_rep_front,p_est_sorted_sym,coTop,coBottom,coFront,coEnd); + + % Initialize vp vector + VP_FRONT = zeros(size(p_est_rep_front)); + + % Create virtual pilots on front + for p = 1:length(p_est_rep_front) + VP_FRONT(p,:) = calculatePilot(p_est_rep_front(p,:), p_est_sorted_sc); + end + + % Transpose for output + VP = [VP VP_FRONT.']; + end + + if CREATE_END + % Repeat last symbol containing pilots + p_est_rep_end = p_est_sorted_sym(end-(noPilotinSym*2)+1:end-noPilotinSym,:); + + % Add subcarrier above and below + p_est_rep_end = addSCs(p_est_rep_end,noPilotinSym); + + % Repeat second symbol containing pilots + p_est_rep_end_second = p_est_sorted_sym(end-noPilotinSym+1:end,:); + + % Add subcarrier above and below + p_est_rep_end_second = addSCs(p_est_rep_end_second,noPilotinSym); + + % Concatenate to create virtual pilots + p_est_rep_end = [p_est_rep_end.' p_est_rep_end_second.'].'; + + % Add seven to normalize symbol indices, mod controls case when + % only 2 pilots per subcarrier + p_est_rep_end(:,2) = p_est_rep_end(:,2) + 7; + + % Remove virtual pilots which are too far + p_est_rep_end = removeTooFar(p_est_rep_end,p_est_sorted_sym,coTop,coBottom,coFront,coEnd); + + % Initialize vp vector + VP_END = zeros(size(p_est_rep_end)); + + % Create virtual pilots on front + for p = 1:length(p_est_rep_end) + VP_END(p,:) = calculatePilot(p_est_rep_end(p,:), p_est_sorted_sc); + end + + % Transpose for output + VP = [VP VP_END.']; + end +end + +function vp = calculatePilot(p_est_rep, p_est_sorted) + % Calculate Euclidean distance between virtual pilot and other + % pilots + ind = sqrt((p_est_rep(1)-p_est_sorted(:,1)).^2+ (p_est_rep(2)-p_est_sorted(:,2)).^2); + % ind = (abs(p_est_rep(1)-p_est_sorted(:,1))+ abs(p_est_rep(2)-p_est_sorted(:,2))); + + % Sort from shortest to longest distance + [~, ind] = sortrows(ind); + + % Take three closest pilots + ind = ind(1:10); + pilots_use = p_est_sorted(ind,:); + + % If first 3 pilots used are in the same subcarrier or symbol then use 4th instead of 3rd + while ((pilots_use(3,1) == pilots_use(2,1)) && (pilots_use(2,1) == pilots_use(1,1)) || (pilots_use(3,2) == pilots_use(2,2)) && (pilots_use(2,2) == pilots_use(1,2))) + pilots_use(3,:) = []; + end + + % Calculate virtual pilot value + vp = calculateVirtualValue(pilots_use(3,:),pilots_use(1,:),pilots_use(2,:),p_est_rep(1),p_est_rep(2)); +end + +function new = calculateVirtualValue(a,b,c,xnew,ynew) + % Calculate vectors + AB = b-a; + AC = c-a; + + % Perform cross product + cro = cross(AB,AC); + + % Break out X, Y and Z plane coeffs + x = cro(1); + y = cro(2); + z = cro(3); + + % Calculate normal in equation + c = sum(cro.*a); + + % Calculate new z value + znew = (c - xnew*x - ynew*y)/z; + + % Return new point + new = [xnew ynew znew]; +end + +function pilots = addSCs(pilots, noPilotinSym) + % Add on extra subcarrier above + pilots(end+1,:) = pilots(1,:); + pilots(end,1) = pilots(end,1) - 6; + + % Add on extra subcarrier below + pilots(end+1,:) = pilots(noPilotinSym,:); + pilots(end,1) = pilots(end,1) + 6; +end + +function pilots = removeTooFar(pilots,p_est_sorted_sym,coTop,coBottom,coFront,coEnd) + pilots((pilots(:,2) < coFront) | (pilots(:,1) < coTop) | (pilots(:,2)> coEnd) | (pilots(:,1) > coBottom),:) = []; + + % Remove virtual pilots which fall within resource grid + removeInd = []; + for i = 1:size(pilots) + if (find(p_est_sorted_sym(:,2) == pilots(i,2))) + removeInd = [removeInd i]; %#ok + end + end + pilots(removeInd,:) = []; +end + +%CreateEdgeVirtualPilots +% Creates virtual pilots on the edges of the resource grid to improve +% interpolation results. Also adds virtual pilots around the current pilot +% estimates. +function vps = createEdgeVirtualPilots(enb,p_est) + + % Determine dimensions of current subframe + dims=lteDLResourceGridSize(enb); + K=dims(1); + L=dims(2); + + % Calculate virtual pilots beyond upper and lower bandwidth edge based + % on pilots present on subcarriers. Also create virtual pilots 1/2 RB + % above and below extent of current pilots. + vps=createDimensionVirtualPilots(p_est,2,unique([-6 min(p_est(1,:)-6) max(p_est(1,:)+6) K+6])); + + % Combine these frequency-direction VPs with original pilots. + temp = [p_est vps]; + + % Calculate virtual pilots beyond start and end of subframe based on + % pilots and frequency-direction VPs present in symbols. Also create + % virtual 1 OFDM symbol before and after extent of current pilots. + vps = [vps createDimensionVirtualPilots(temp,1,unique([-1 min(p_est(2,:)-1) max(p_est(2,:)+1) L+1]))]; + +end + +% dim=1 adds VPs in time +% dim=2 adds VPs in frequency +function vps = createDimensionVirtualPilots(p_est,dim,points) + + vps=[]; + + pilots = unique(p_est(dim,:)); + for i = pilots + x=find(p_est(dim,:)==i); + rep=1+double(length(x)==1); + adj=(1:rep)-rep; + pil = interp1(p_est(3-dim,x)+adj,repmat(p_est(3,x),1,rep),points,'linear','extrap'); + for j=1:length(points) + vps = [ vps [i;points(j);pil(j)] ]; %#ok + end + end + + if (dim==2) + vps = vps([2 1 3],:); + end + +end + diff --git a/matlab/tests/lteDLChannelEstimate3.m b/matlab/tests/lteDLChannelEstimate3.m new file mode 100644 index 000000000..7b0618002 --- /dev/null +++ b/matlab/tests/lteDLChannelEstimate3.m @@ -0,0 +1,1204 @@ +%lteDLChannelEstimate Downlink channel estimation +% [HEST NOISEEST] = lteDLChannelEstimate(...) returns HEST, the estimated +% channel between each transmit and receive antenna and NOISEEST, an +% estimate of the noise power spectral density on the reference signal +% subcarriers. +% +% HEST is an M-by-N-by-NRxAnts-by-CellRefP (optionally +% M-by-N-by-NRxAnts-by-NLayers for UE-specific beamforming transmission +% schemes) array where M is the number of subcarriers, N is the number of +% OFDM symbols, NRxAnts is the number of receive antennas, CellRefP is +% the number of cell-specific reference signal antenna ports and NLayers +% is the number of transmission layers. Using the reference signals, an +% estimate of the power spectral density of the noise present on the +% estimated channel response coefficients is returned. +% +% [HEST NOISEEST] = lteDLChannelEstimate(ENB, RXGRID) returns the +% estimated channel coefficients using the method described in +% TS36.104/TS36.141 Annex E/F for the purposes of transmitter EVM +% testing. +% +% ENB is a structure and must contain the following fields: +% NDLRB - Number of downlink resource blocks +% CellRefP - Number of cell-specific reference signal antenna ports +% (1,2,4) +% NCellID - Physical layer cell identity +% NSubframe - Subframe number +% CyclicPrefix - Optional. Cyclic prefix length +% ('Normal'(default),'Extended') +% DuplexMode - Optional. Duplex mode ('FDD'(default),'TDD') +% Only required for 'TDD' duplex mode: +% TDDConfig - Optional. Uplink/Downlink Configuration (0...6) +% (default 0) +% SSC - Optional. Special Subframe Configuration (0...9) +% (default 0) +% Only required for CEC.Reference='CSIRS' below: +% CSIRefP - Number of CSI-RS antenna ports (1,2,4,8) +% CSIRSConfig - CSI-RS configuration index (TS 36.211 Table +% 6.10.5.2-1) +% CSIRSPeriod - Optional. CSI-RS subframe configuration: +% ('On'(default),'Off',Icsi-rs,[Tcsi-rs Dcsi-rs]) +% +% RXGRID is a 3-dimensional M-by-N-by-NRxAnts array of resource elements. +% The second dimension of RXGRID can contain any whole number of +% subframes worth of OFDM symbols i.e. for normal cyclic prefix each +% subframe contains 14 OFDM symbols, therefore N is a multiple of 14. +% Note: to adhere to the estimation method defined in TS36.104/TS36.141, +% RXGRID must contain 10 subframes. +% +% [HEST NOISEEST] = lteDLChannelEstimate(ENB, CEC, RXGRID) returns the +% estimated channel using the method and parameters defined by the +% configuration structure CEC. +% +% CEC is a structure which defines the type of channel estimation +% performed. CEC must contain a set of the following fields: +% PilotAverage - Type of pilot averaging ('TestEVM', 'UserDefined') +% FreqWindow - Size of window in resource elements used to average +% over frequency. +% TimeWindow - Size of window in resource elements used to average +% over time +% InterpType - Type of 2D interpolation used(see griddata for types) +% InterpWindow - Interpolation window type: ('Causal','Non-causal', +% 'Centred','Centered'). Note: 'Centred' and 'Centered' +% are equivalent +% InterpWinSize - Interpolation window size: +% 'Causal','Non-causal' - any number >=1 +% 'Centred','Centered' - odd numbers >=1 +% +% The 'TestEVM' pilot averaging will ignore other structure fields in +% CEC, and the method follows that described in TS36.104/TS36.141 Annex +% E/F for the purposes of transmitter EVM testing. +% +% The 'UserDefined' pilot averaging uses a rectangular kernel of size +% CEC.FreqWindow-by-CEC.TimeWindow and performs a 2D filtering operation +% upon the pilots. Note that pilots near the edge of the resource grid +% will be averaged less as they have no neighbors outside of the grid, +% or a limited number of neighbors outside of the grid obtained by the +% creation of virtual pilots. +% +% [HEST NOISEEST] = lteDLChannelEstimate(ENB, CHS, CEC, RXGRID) returns +% the estimated channel given cell-wide settings structure ENB, PDSCH +% transmission configuration CHS and channel estimator configuration CEC. +% CHS must include the following fields: +% TxScheme - Transmission scheme, one of: +% 'Port0' - Single-antenna port, Port 0 +% 'TxDiversity' - Transmit diversity scheme +% 'CDD' - Large delay CDD scheme +% 'SpatialMux' - Closed-loop spatial multiplexing scheme +% 'MultiUser' - Multi-user MIMO scheme +% 'Port5' - Single-antenna port, Port 5 +% 'Port7-8' - Single-antenna port, port 7 (when +% NLayers=1); Dual layer transmission, port 7 +% and 8 (when NLayers=2) +% 'Port8' - Single-antenna port, Port 8 +% 'Port7-14' - Up to 8 layer transmission, ports 7-14 +% PRBSet - A 1- or 2-column matrix, containing the 0-based Physical +% Resource Block indices (PRBs) corresponding to the resource +% allocations for this PDSCH. +% RNTI - Radio Network Temporary Identifier (16-bit) +% +% For the 'Port5', 'Port7-8', 'Port8' and 'Port7-14' transmission +% schemes, the channel estimator configuration structure CEC contains the +% following additional field: +% Reference - Optional. Specifies point of reference (signals to +% internally generate) for channel estimation. +% ('DMRS'(default), 'CSIRS') +% +% For the 'Port5', 'Port7-8', 'Port8' and 'Port7-14' transmission +% schemes, with Reference='DMRS', the channel estimation is performed +% using UE-specific reference signals and the returned channel estimate +% will be of size M-by-N-by-NRxAnts-by-NLayers. Alternatively, with +% Reference='CSIRS' the channel estimation is performed using the CSI +% reference signals (CSI) and the returned channel estimate will be of +% size M-by-N-by-NRxAnts-by-CSIRefP. For other transmission schemes the +% channel estimation is performed using cell-specific reference signals +% and the channel estimate will be of size M-by-N-by-NRxAnts-by-CellRefP. +% Note that CSI-RS based channel estimation and hence Reference='CSIRS' +% is strictly only valid within the standard for the 'Port7-14' +% transmission scheme. The optional CSIRSPeriod parameter controls the +% downlink subframes in which CSI-RS will be present, either always 'On' +% or 'Off', or defined by the scalar subframe configuration index Icsi-rs +% (0...154) or the explicit subframe periodicity and offset pair [Tcsi-rs +% Dcsi-rs] (TS 36.211 Section 6.10.5.3). +% +% For the 'Port7-8' and 'Port7-14' transmission schemes with +% 'UserDefined' pilot averaging, if CEC.TimeWindow = 2 or 4 and +% CEC.FreqWindow=1 the estimator will enter a special case where an +% averaging window of 2 or 4 pilots in time will be used to average the +% pilot estimates; the averaging is always applied across 2 or 4 pilots, +% regardless of their separation in OFDM symbols. This operation ensures +% that averaging is always done on 2 or 4 pilots. This provides the +% appropriate "despreading" operation required for the case of UE RS +% ports / CSI-RS ports which occupy the same time/frequency locations but +% use different orthogonal covers to allow them to be differentiated at +% the receiver. For the CSI-RS and any number of configured CSI-RS ports +% (given by ENB.CSIRefP), the pilot REs occur in pairs, one pair per +% subframe, that require averaging with CEC.TimeWindow=2 and will result +% in a single estimate per subframe. For the UE RS with between 1 and 4 +% layers (given by CHS.NLayers), the pilot REs occur in pairs, repeated +% in each slot, that require averaging with CEC.TimeWindow=2 and will +% result in two estimates per subframe, one for each slot; for between 5 +% and 8 layers, the pairs are distinct between the slots of the subframe +% and the required averaging is CEC.TimeWindow=4, resulting in one +% estimate per subframe. +% +% Example: +% Transmit RMC R.12 (4-antenna transmit diversity), model the propagation +% channel by combining all transmit antennas onto one receive antenna, +% OFDM demodulate and finally channel estimate. +% +% enb = lteRMCDL('R.12'); +% cec = struct('FreqWindow',1,'TimeWindow',1,'InterpType','cubic',... +% 'PilotAverage','UserDefined','InterpWinSize',3,... +% 'InterpWindow','Causal'); +% txWaveform = lteRMCDLTool(enb,[1;0;0;1]); +% rxWaveform = sum(txWaveform,2); +% rxGrid = lteOFDMDemodulate(enb,rxWaveform); +% hest = lteDLChannelEstimate(enb,cec,rxGrid); +% +% See also lteOFDMDemodulate, lteEqualizeMMSE, lteEqualizeZF, +% lteDLPerfectChannelEstimate, griddata. + +% Copyright 2009-2013 The MathWorks, Inc. + +function [H_EST, NoisePowerEst, AvgEstimates, Estimates] = lteDLChannelEstimate2(varargin) + + if(isstruct(varargin{2})) + if(isstruct(varargin{3})) + PDSCH = varargin{2}; + CEC = varargin{3}; + RXGRID = varargin{4}; + else + PDSCH = []; + CEC = varargin{2}; + RXGRID = varargin{3}; + avgrefs = varargin{4}; + end + else + % If no configuration structure then use the TestEVM method + PDSCH = []; + CEC.PilotAverage = 'TestEVM'; + RXGRID = varargin{2}; + end + + ENB = varargin{1}; + % Get dimensions of resource grid + Dims = lteDLResourceGridSize(ENB); + K = Dims(1); + Lsf = Dims(2); + + % Determine number of subframes + nsfs = size(RXGRID,2)/Dims(2); + + % Determine number of Tx- and RxAntennas + NRxAnts = size(RXGRID,3); + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port5' 'Port7-8' 'Port8' 'Port7-14'}))) + if (~isfield(CEC,'Reference')) + CEC.Reference='DMRS'; + end + if (strcmpi(CEC.Reference,'DMRS')==1) + NTx = PDSCH.NLayers; + elseif (strcmpi(CEC.Reference,'CSIRS')==1) + NTx = ENB.CSIRefP; + else + error('lte:error','Reference must be "DMRS" or "CSIRS", see help for details.'); + end + else + if (~isfield(CEC,'Reference')) + CEC.Reference='CellRS'; + end + NTx = ENB.CellRefP; + end + + % Initialize size of estimated channel grid + H_EST = zeros([size(RXGRID,1) size(RXGRID,2) size(RXGRID,3) NTx]); + + % Preallocate noise power estimate vector for speed + noiseVec = zeros(size(NRxAnts,NTx)); + + if (NTx == 4) + nn=3; + else + nn=NTx; + end + RefEstimates=zeros(nn,2*2*2*ENB.NDLRB*10); + + for rxANT = 1:NRxAnts + for nTxAntLayer = 1:NTx + + % Extract pilot symbols from received grid and calculate least + % squares estimate + [ls_estimates,specialvec,DwPTS]= GetPilotEstimates(ENB,nTxAntLayer-1,RXGRID(:,:,rxANT),PDSCH,CEC); + + Estimates(nTxAntLayer,1:length(ls_estimates(3,:))) = ls_estimates(3,:); + + % Average the pilots as defined by CEC.PilotAverage, this can + % be 'UserDefined' or 'TestEVM'. + % Note: Setting the window size to 1x1 is equivalent to no + % averaging. It is recommended that no averaging should be + % carried out when a high SNR is present, as this would have an + % adverse effect on the least squares estimates + if strcmpi(CEC.PilotAverage,'UserDefined')&&(CEC.FreqWindow==1)&&(CEC.TimeWindow==1) + P_EST = ls_estimates; + ScalingVec = 0; + else + [P_EST, ScalingVec]= PilotAverage(PDSCH,CEC,H_EST,ls_estimates); + end + + %AvgEstimates(nTxAntLayer,1:length(P_EST(3,:))) = P_EST(3,:); + + P_EST(3,:)=avgrefs; + + if strcmpi(CEC.PilotAverage,'TestEVM') + % Channel estimation as defined in TS36.141 Annex F.3.4 is + % performed, the pilots are averaged in time and frequency + % and the resulting vector is linearly interpolated. + %---------------------- + % Interpolate eq coefficients extrapolating to account for + % missing pilots at the edges + interpEqCoeff = interp1(find(P_EST~=0),P_EST(P_EST~=0),(1:length(P_EST)).','linear','extrap'); + % --------------------- + % The DC carrier needs to be accounted for during + % interpolation + % Get pilot symbol index values + ind = lteCellRSIndices(ENB,0); + % Using initial index value calculate position of pilot + % symbols within interpolated coefficients + if ind(1)-3>0 + vec = (ind(1)-3):3:K; + else + vec = ind(1):3:K; + end + % Determine the index values of the pilots either side of + % the DC carrier + preDC = vec(length(vec)/2); + postDC = vec(1+length(vec)/2); + % Interpolate the section using the DC carrier and insert + % this into the interpolated vector containing the + % coefficients + interpEqCoeffTemp = interp1(1:4:5,[interpEqCoeff(preDC) interpEqCoeff(postDC)],1:5,'linear'); + interpEqCoeff(preDC:K/2) = interpEqCoeffTemp(1:length(interpEqCoeff(preDC:K/2))); + interpEqCoeff(1+K/2:postDC) = interpEqCoeffTemp(2+length(interpEqCoeff(preDC:K/2)):end); + %---------------------- + % Generate grid by replicating estimated equalizer channel + % coefficients + H_EST(:,:,rxANT,nTxAntLayer) = repmat(interpEqCoeff,1,Lsf*nsfs); + % The value of the noise averaged pilot symbols are placed + % into P_EST(3,:) matrix and these are used to determine + % the noise power present on the pilot symbol estimates of + % the channel. + P_EST = ls_estimates; + tempGrid = squeeze(H_EST(:,:,rxANT,nTxAntLayer)); + P_EST(3,:) = tempGrid(sub2ind(size(H_EST),P_EST(1,:),P_EST(2,:))); + + else + % Channel estimation is performed by analyzing up to + % 'CEC.InterpWinSize' subframes together and estimating + % virtual pilots outwith the bounds of the subframe. + % Averaging window can be set to Causal, Non-causal or + % Centred, Centered. + + % Using the desired subframeAveraging type the initial + % position of the window is determined + if(CEC.InterpWinSize>=1) + if strcmpi(CEC.InterpWindow,'Centred')||strcmpi(CEC.InterpWindow,'Centered') + % For Centered averaging the window size must be odd + if ~mod(CEC.InterpWinSize,2) + error('lte:error','Window size must be odd for centered window type'); + end + % For Centered windowing both current,future and past + % data are used to estimate current channel coefficients + x = floor(CEC.InterpWinSize/2); + y = floor(CEC.InterpWinSize/2); + elseif strcmpi(CEC.InterpWindow,'Causal') + % For causal windowing current and past data are used + % to estimate current subframe channel coefficients + x = 0; + y = CEC.InterpWinSize-1; + elseif strcmpi(CEC.InterpWindow,'Non-causal') + % For non-causal windowing current and future data are + % used to estimate current subframe channel + % coefficients + x = CEC.InterpWinSize-1; + y = 0; + else + error('lte:error','Channel estimation structure field InterpWindow must be one of the following: Causal,Non-causal,Centred,Centered'); + end + else + error('lte:error','InterpWinSize cannot be less than 1'); + end + + % Interpolate using averaged pilot estimates and defined + % interpolation settings + for sf = 0:nsfs-1 + if specialvec(sf+1)=='U' + H_EST(:,sf*Lsf+1:(sf+1)*Lsf,rxANT,nTxAntLayer) = NaN; + else + % Extract the pilots from required subframes + p_use = P_EST; + p_use(:,p_use(4,:)>(sf+1)+x)=[]; + p_use(:,p_use(4,:)<(sf+1)-y)=[]; + + % Account for DC offset + p_use(1,(p_use(1,:)>K/2)) = p_use(1,(p_use(1,:)>K/2))+1; + + if strcmpi(CEC.InterpType,'Cubic')||strcmpi(CEC.InterpType,'Linear') + % VPVEC is used to determine if virtual pilots are + % needed at the beginning or end of the + % interpolation window, if the current subframe,sf, + % is located at beginning or end of the window then + % virtual pilots are created accordingly + vpvec = unique(p_use(4,:)); + + if (strcmp(CEC.Reference,'CellRS')==1) + % Create virtual pilots and append to vector + % containing pilots using cell-specific RS + % methodology. + vps = createVirtualPilots(ENB,p_use(1:3,:),1,1,(sf+1==min(vpvec)),(sf+1==max(vpvec))); + p_use = [p_use(1:3,:) vps]; + else + % Create edge virtual pilots suitable for + % UE RS which can have partial bandwidth, + % or CSI-RS. + if (~isempty(p_use)) + vps = createEdgeVirtualPilots(ENB,p_use(1:3,:)); + p_use = [p_use(1:3,:) vps]; + end + end + end + % Perform 2D interpolation + % Interpolation is carried out on a (K+1)-by-L + % matrix to account for DC offset being added in + if (~isempty(p_use)) + Htemp = griddata(p_use(2,:)-Lsf*sf,p_use(1,:),p_use(3,:),1:Lsf,(1:K+1)',CEC.InterpType); %#ok + % Remove DC offset + Htemp(1+(K/2),:) = []; + + if specialvec(sf+1)=='S' + Htemp(:,DwPTS+1:Lsf) = NaN; + end + H_EST(:,sf*Lsf+1:(sf+1)*Lsf,rxANT,nTxAntLayer) = Htemp; + + if isnan(H_EST) + error('lte:error','H_EST NaN'); + end + end + end + end + end + if nTxAntLayer<3 || (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port7-8' 'Port7-14'}))) + % The noise level present can be determined using the noisy + % least squares estimates of the channel at pilot symbol + % locations and the noise averaged pilot symbol estimates of the + % channel + noise = ls_estimates(3,~isnan(ls_estimates(3,:))) - P_EST(3,~isnan(P_EST(3,:))); + if strcmpi(CEC.PilotAverage,'UserDefined') + noise = sqrt(ScalingVec./(ScalingVec+1)).*noise; + end + + % Additional averaging for noise estimation in LTE-A case, + % to suppress interference from orthogonal sequences on + % other antennas in same time-frequency locations. + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port7-8' 'Port7-14'})) &&... + (CEC.TimeWindow==2 || CEC.TimeWindow==4) && CEC.FreqWindow==1) + if (~isempty(ls_estimates)) + if (strcmpi(CEC.Reference,'DMRS')==1) + temp=[]; + pilotSC = unique(ls_estimates(1,:)); + for i = pilotSC + x=find(ls_estimates(1,:)==i); + x=[x(1) x(end)]; + temp=[temp mean(noise(x))*2]; %#ok + end + noise=temp; + end + end + end + + % Taking the variance of the noise present on the pilot symbols + % results in a value of the noise power for each transmit and + % receive antenna pair + if (isempty(noise)) + noiseVec(rxANT,nTxAntLayer)=NaN; + else + noiseVec(rxANT,nTxAntLayer) = mean(noise.*conj(noise)); + end + end + end + end + + % The mean of the noise power across all the transmit/receive antenna + % pairs is used as the estimate of the noise power + NoisePowerEst = mean(mean(noiseVec)); + +end + +% GetPilotEstimates Obtain the least squares estimates of the reference +% signals +% [ls_estimates] = GetPilotEstimates(ENB,nTxAntLayer,RXGRID) Extracts the +% reference signals for a specific transmit/receive antenna pair and +% calculates their least squares estimate. The results are placed in a +% 3xNp matrix containing the subcarrier and OFDM symbol location, row and +% column subscripts, and value. Np is the number of cell specific +% reference (pilot) symbols per resource grid +% +% RXGRID is an MxN array where M is the number of subcarriers, +% N is the number of OFDM symbols. Dimension M must be +% 12*NDLRB where NDLRB must be {6,15,25,50,75,100}. +% Dimension N must be a multiple of number of symbols in a subframe L, +% where L=14 for normal cyclic prefix and L=12 for extended cyclic +% prefix. +% +% nTxAntLayer defines which transmit antennas' or layers' pilot symbols to extract +% +% ENB is a structure and must contain the following fields: +% NDLRB - Number of downlink resource blocks +% NCellID - Physical layer cell identity +% CellRefP - Number of transmit antenna ports {1,2,4} +% CyclicPrefix - Optional. Cyclic prefix length{'Normal'(default),'Extended'} +% +% +% Example +% Return the least squares estimate of the pilots symbols in a 3xNp array +% where the first row is the subcarrier index,k the second row is the OFDM +% symbol,l and the third row defines the least squares estimate of the +% pilot located at that position. +% +% enb=struct('NDLRB',6,'CellRefP',1,'NCellID',0,'CyclicPrefix','Normal'); +% rxGrid=ones(lteDLResourceGridSize(enb)); +% nTxAntLayer=0; +% ls_estimates = GetPilotEstimates(enb,nTxAntLayer,rxGrid).' +% ans = +% 1.0000 1.0000 -0.7071 - 0.7071i +% 7.0000 1.0000 -0.7071 - 0.7071i +% 13.0000 1.0000 -0.7071 - 0.7071i +% . . . +% . . . +% . . . +% 58.0000 12.0000 -0.7071 + 0.7071i +% 64.0000 12.0000 0.7071 - 0.7071i +% 70.0000 12.0000 0.7071 - 0.7071i + +% Copyright 2009-2010 The MathWorks, Inc. +function [ls_estimates,specialvec,DwPTS] = GetPilotEstimates(ENB,nTxAntLayer,RXGRID,PDSCH,CEC) + +% Get dimensions of resource grid +Dims = lteDLResourceGridSize(ENB); +nsfs = size(RXGRID,2)/Dims(2); +K = Dims(1); +L = Dims(2); + +if (rem(nsfs,1)~=0) + error('lte:error','The received grid input must contain a whole number of subframes.'); +end + +if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port5' 'Port7-8' 'Port8' 'Port7-14'}))) + if (strcmp(CEC.Reference,'DMRS')==1) + % create UE RS indices + PDSCH.NTxAnts=0; + PDSCH.W=[]; + linearIndSF = lteDMRSIndices(ENB,PDSCH,'mat'); + else + % create CSI-RS indices and symbols and deal with zero removal + linearIndSF = lteCSIRSIndices(ENB,'mat'); + CsiRS = lteCSIRS(ENB,'mat'); + CsiRS = CsiRS(:,nTxAntLayer+1); + linearIndSF(CsiRS==0,:)=[]; + CsiRS(CsiRS==0)=[]; + end + % If more than one TxAntenna is used we need to adjust indices values so that + % they start in first antenna plane + linearIndSF = linearIndSF(:,nTxAntLayer+1) - (K*L*nTxAntLayer); +else + % Extract linear indices of reference signals for particular TxAntenna + dumENB = ENB; + dumENB.NSubframe = 0; + linearIndSF = lteCellRSIndices(dumENB,nTxAntLayer); + % If more than one TxAntenna is used we need to adjust indices values so that + % they start in first antenna plane + linearIndSF = linearIndSF - (K*L*nTxAntLayer); +end + +% Create grid for pilot symbols +linearIndGrid = zeros(size(linearIndSF,1),nsfs); +idealPilotGrid = zeros(size(linearIndSF,1),nsfs); +offset = double(ENB.NSubframe); +specialvec = char(zeros(1,0)); +DwPTS = 0; +for sf = ENB.NSubframe:ENB.NSubframe+nsfs-1 + ENB.NSubframe = mod(sf,10); + linearIndGrid(:,sf-offset+1) = (linearIndSF+(sf*K*L)); + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port5' 'Port7-8' 'Port8' 'Port7-14'}))) + if (strcmp(CEC.Reference,'DMRS')==1) + % create UE RS symbols + DMRS = lteDMRS(ENB,PDSCH,'mat'); + idealPilotGrid(1:length(DMRS),sf-offset+1) = DMRS(:,nTxAntLayer+1); + else + % use CSI-RS symbols created earlier + idealPilotGrid(1:length(CsiRS),sf-offset+1) = CsiRS; + end + else + CellRS = lteCellRS(ENB, nTxAntLayer); + idealPilotGrid(1:length(CellRS),sf-offset+1) = CellRS; + end + dupinfo = lteDuplexingInfo(ENB); + if strcmpi(dupinfo.SubframeType,'Special') + specialvec = [specialvec 'S']; %#ok + DwPTS = dupinfo.NSymbolsDL; + elseif strcmpi(dupinfo.SubframeType,'Uplink') + specialvec = [specialvec 'U']; %#ok + elseif strcmpi(dupinfo.SubframeType,'Downlink') + specialvec = [specialvec 'D']; %#ok + else + error('lte:error','Invalid SubframeType, in structure dupinfo'); + end +end + +if (~isempty(idealPilotGrid)) + ulvec = find(idealPilotGrid(size(idealPilotGrid,1),:)==0); + linearIndGrid(:,ulvec) = []; + idealPilotGrid(:,ulvec) = []; +else + ulvec = []; +end + +linearIndGrid = linearIndGrid - (offset*K*L); +% Extract the row and column subscripts of the pilot symbols for entire +% grid +[p_estSC, p_estSym] = ind2sub(size(RXGRID),linearIndGrid); + +% Calculate least squares channel estimates at pilot locations +p_est = RXGRID(linearIndGrid)./idealPilotGrid; + +% Create a vector - [k;l;p_est] +sfref = repmat((1:nsfs),length(linearIndSF),1); +sfref(:,ulvec) = []; +ls_estimates = [double(p_estSC(:).') ; double(p_estSym(:).') ; p_est(:).';double(sfref(:)).']; + +end + +%PilotAverage Average reference signals +% [P_EST] = PilotAverage(PDSCH,CEC,H_EST,LS_EST) performs a moving average of pilot +% symbols +% +% LS_EST is a 3xNp matrix containing the least square estimates of the +% pilots symbols and their column and row indices within the received +% grid. +% LS_EST = [k;l;p_est] +% H_EST is an MxN matrix and defines the size of the grid that the +% averaging will be performed on +% +% CEC is a structure which defines the type of channel estimation +% performed. CEC must contain a set of the following fields: +% PilotAverage - Type of pilot averaging {'TestEVM', 'UserDefined'} +% FreqWindow - Size of window used to average in frequency in +% resource elements. +% TimeWindow - Size of window used to average in time in resource +% elements +% InterpType - Type of 2D interpolation used(see griddata for types) +% +% The dimensions of the averaging window are defined in structure CEC. +% The window is defined in terms of Resource Elements, and depending on +% the size of the averaging window, averaging will be performed in either +% the time or frequency direction only, or a combination of both creating +% a square/rectangular window. The pilot to be averaged will always be +% placed at the center of the window, therefore the window size must be +% an odd number. +% +% Frequency Direction: 9x1 Time Direction: 1x9 +% +% x +% x +% x +% x +% P x x x x P x x x x +% x +% x +% x +% x +% +% +% Square window: 9x9 +% x x x x x x x x x +% P x x x x x x P x +% x x x x x x x x x +% x x x x x x x x x +% x x x x P x x x x +% x x x x x x x x x +% x x x x x x x x x +% P x x x x x x P x +% x x x x x x x x x +% +% Performing EVM compliance testing as per TS36.141 AnnexF.3.4, requires +% time averaging be done over 10 subframes(i.e. 1 Frame) across each +% pilot symbol carrying subcarrier creating a TotalNumberPilotsx1 vector. +% This is then frequency averaged using a moving window average with a window +% size of 19. This type of averaging can be performed by setting the +% PilotAverage type in structure CEC to 'TestEVM'. + +% Copyright 2009-2013 The MathWorks, Inc. + +function [P_EST, scalingVec] = PilotAverage(PDSCH,CEC,H_EST,P_EST) + + switch lower(CEC.PilotAverage) + + case 'userdefined' + if(CEC.FreqWindow<1)||(CEC.TimeWindow<1) + error('lte:error','Frequency and time averaging window size cannot be less than 1'); + end + + if (~isempty(PDSCH) && any(strcmpi(PDSCH.TxScheme,{'Port7-8' 'Port7-14'})) &&... + (CEC.TimeWindow==2 || CEC.TimeWindow==4) && CEC.FreqWindow==1) + + % Perform averaging in time direction using a moving window + % of size N = CEC.TimeWindow + N = CEC.TimeWindow; + + if (~isempty(P_EST)) + % Average only subcarriers which contain RS symbols. Extract RS + % symbols on a per subcarrier basis and use window to average. + for subcarrier = unique(P_EST(1,:)) + + % Store DRS symbols from relevant slot in vector for easy access + symbVec = (P_EST(3,P_EST(1,:)==subcarrier)).'; + + % Define temporary vector to store the averaged values + avgVec = zeros(size(symbVec)); + + % Check length of input at least equal to size of window; + % if not, reduce window size. Despreading will be incomplete + % here + if (length(symbVec) < N) + N = length(symbVec); + end + + % Determines position of window w.r.t element being averaged, + % and performs the averaging + for symbNo = 1:length(symbVec) + if (symbNo-(N/2+1)<=0) + avgVec(symbNo)= sum(symbVec(1:N))/N; + else + avgVec(symbNo) = sum(symbVec(end-(N-1):end))/N; + end + end + + % Update P_EST with averaged values + P_EST(3,P_EST(1,:)==subcarrier)= avgVec; + end + + end + + % This vector is used to scale the noise by the number of averaging + % elements in the window + scalingVec = ones(size(P_EST(3,:)))*N; + + else + + if (strcmpi(CEC.InterpWindow,'Centred')||strcmpi(CEC.InterpWindow,'Centered')) + if (~mod(CEC.FreqWindow,2)||~mod(CEC.TimeWindow,2)) + error('lte:error','Window size must be odd in time and frequency for centred/centered window type'); + end + end + + % Define number of subcarriers + K = size(H_EST,1); + + % Define an empty resource grid with the DC offset subcarrier + % inserted + grid = zeros(size(H_EST,1)+1,size(H_EST,2)); + + % Account for DC offset + P_EST(1,(P_EST(1,:)>K/2)) = P_EST(1,(P_EST(1,:)>K/2))+1; + + % Place the pilot symbols back into the received grid with the + % DC offset subcarrier in place + grid(sub2ind(size(grid),P_EST(1,:),P_EST(2,:))) = P_EST(3,:); + + % Define convolution window + kernel = ones(CEC.FreqWindow,CEC.TimeWindow); + + % Perform convolution + grid = conv2(grid,kernel,'same'); + + % Extract only pilot symbol location values and set the rest of + % the grid to zero + tempGrid = zeros(size(grid)); + tempGrid(sub2ind(size(grid),P_EST(1,:),P_EST(2,:))) = grid(sub2ind(size(grid),P_EST(1,:),P_EST(2,:))); + grid = tempGrid; + + % Normalize pilot symbol values after convolution + [grid, scalingVec] = normalisePilotAverage(CEC,P_EST,grid); + + % Place averaged values back into pilot symbol matrix + P_EST(3,:) = grid(sub2ind(size(tempGrid),P_EST(1,:),P_EST(2,:))); + + % Remove DC offset + P_EST(1,(P_EST(1,:)>K/2)) = P_EST(1,(P_EST(1,:)>K/2))-1; + + end + + + case 'testevm' + % As defined in TS36.141 Annex F.3.4 pilots are averaged in + % time across all pilot carrying subcarriers. The resulting + % vector is averaged in frequency direction with a moving + % averaging window of 19. + + % Get phase (theta) and magnitude(radius) of complex estimates + [theta,radius] = cart2pol(real(P_EST(3,:)),imag(P_EST(3,:))); + + % Declare vector to temporarily store averaged phase and + % magnitude + phasemagVec = zeros(size(H_EST,1),2); + + % Perform averaging in time and frequency for phase and + % magnitude separately and recombine at end + for phaseOrMag = 1:2 + % Define temporary P_EST vector and set estimates row to + % phase then magnitude + P_EST_temp = P_EST; + if phaseOrMag==1 + P_EST_temp(3,:) = theta; + elseif phaseOrMag==2 + P_EST_temp (3,:) = radius; + end + + % Declare averaging vector + avgVec = zeros(size(H_EST,1),1); + vec = []; + + % Average across pilot symbol carrying subcarriers in time, + % unwrapping performed on a subcarrier basis + for i = 1:size(H_EST,1) + if ~isempty(P_EST(3,(P_EST(1,:)==i))) + if phaseOrMag==1 + avgVec(i) = mean(unwrap(P_EST_temp(3,(P_EST(1,:)==i)))); + + elseif phaseOrMag==2 + avgVec(i) = mean(P_EST_temp(3,(P_EST(1,:)==i))); + end + vec = [vec i]; %#ok + end + end + + % Remove subcarriers which contained no pilot symbols, then + % perform frequency averaging + avgVec = avgVec(vec); + + % Perform averaging in frequency direction using a moving + % window of size N, N = 19 in TS36.141 Annex F.3.4 + % Performs a moving average of window size N. At the edges, where less than + % N samples are span the window size is reduce to span 1, 3, 5, 7 ... + % samples. + N = 19; + + % avgVec must be a column vector + if (size(avgVec,2) > 1) + error('lte:error','Input to moving average function must be a column vector.') + end + + % N max window size 19 as defined in Annex F.3.4 + if ~mod(N,2) + error('lte:error','Window size N must be odd'); + end + + if (length(avgVec) < N) + % sprintf ('Input signal must have at least %d elements',N); + error('lte:error','Input signal must have at least %d elements',N); + end + + % Use filter to perform part of the averaging (not normalized) + zeroPad = zeros(N-1,1); + data = [avgVec; zeroPad]; + weights = ones(N,1); + freqAvg = filter(weights,1,data); + + % Remove unwanted elements + removeIdx = [ 2:2:N ( length(data)-(2:2:N)+1 )]; + freqAvg(removeIdx) = []; + + % Normalization factor + normFactor = [1:2:N-2 N*ones(1,length(freqAvg)-(N-1)) N-2:-2:1]'; + freqAvg = freqAvg./normFactor; + + % Place frequency averaged pilots into temporary storage + % vector + phasemagVec(vec,phaseOrMag) = freqAvg; + + end + % Convert averaged symbols back to Cartesian coordinates + [X,Y] = pol2cart(phasemagVec(:,1),phasemagVec(:,2)); + P_EST = complex(X,Y); + scalingVec = []; + + otherwise + error('lte:error','Channel estimation structure field PilotAverage must be one of the following: "UserDefined" or "TestEVM"'); + + end +end + +function [avgGrid, scalingVec] = normalisePilotAverage(CEC,p_est,grid) + + % Determine total number of pilots within half a subframe + nPilots = length(p_est); + + avgGrid = zeros(size(grid)); + scalingVec = zeros(1,size(p_est,2)); + + for n = 1:nPilots + % Determine in which subcarrier and OFDM symbol pilot is located + sc = p_est(1,n); + sym = p_est(2,n); + % Determine number of REs to look at either side of pilot + % symbol in time and frequency + half_freq_window = floor(CEC.FreqWindow/2); + half_time_window = floor(CEC.TimeWindow/2); + % Set the location of the window at the back of the pilot + % to be averaged in frequency direction + upperSC = sc-half_freq_window; + % If this location is outwith the grid dimensions set it to + % lowest subcarrier value + if upperSC<1 + upperSC = 1; + end + % Set the location of the window in front of the pilot to + % be averaged in frequency direction + lowerSC = sc+half_freq_window; + % If this location is outwith the grid dimensions set it to + % highest subcarrier value + if lowerSC>size(grid,1) + lowerSC = size(grid,1); + end + % Set the location of the window in front of the pilot to + % be averaged in time direction + leftSYM = sym-half_time_window; + % If this location is outwith the grid dimensions set it to + % lowest OFDM symbol value + if leftSYM<1 + leftSYM = 1; + end + % Set the location of the window at the back of the pilot + % to be averaged in time direction + rightSYM = sym+half_time_window; + % If this location is outwith the grid dimensions set it to + % highest OFDM symbol value + if rightSYM>size(grid,2) + rightSYM = size(grid,2); + end + + % Define the window to average using the determined + % subcarrier and OFDM symbol values + avgVec = grid(upperSC:lowerSC,leftSYM:rightSYM); + + % Remove the zero values from the window so that the average is + % calculated using only valid pilots + avgVec = avgVec(avgVec~=0); + % Average the desired pilot using all pilots within + % averaging window + if isempty(avgVec) + avgVec = 0; + end + + avgGrid(sc,sym) = grid(sc,sym)/numel(avgVec); + scalingVec(n) = numel(avgVec); + end +end + +%createVirtualPilots Create virtual pilots. +% [VP] = createVirtualPilots(..) returns a 3xNvp matrix which contains +% the time/frequency and least squares estimate of virtual pilots. +% Nvp is the number of virtual pilots. +% +% The inputs are: +% +% ENB structure which must contain the following fields: +% NDLRB - Number of downlink resource blocks +% CyclicPrefix - Optional. Cyclic prefix length {'Normal'(default),'Extended'} +% +% P_EST - A matrix containing the pilots within the grid. The matrix must +% have three rows. Each column contains the information for each pilot: +% [subcarrier_indices symbol_indices pilot_channel_estimate] +% +% ANTENNA - The current antenna port +% CREATE_TOP,CREATE_BOTTOM,CREATE_FRONT,CREATE_END - a value {0 or 1} +% indicating whether to return virtual pilots in each position. A 1 +% indicates virtual pilots should be returned, a 0 indicates they should +% not be returned. +% +% Example: +% +% Create virtual pilots on the top of a grid of pilots (estimates) +% VP_TOP = createVirtualPilots(ENB,estimates,0,1,0,0,0); +% +% Add virtual pilots to current pilot estimates for interpolation +% estimates = [estimates VP_TOP]; + +% Copyright 2009-2010 The MathWorks, Inc. +function VP = createVirtualPilots(ENB,P_EST,CREATE_TOP,CREATE_BOTTOM,CREATE_FRONT,CREATE_END) + + % Initialize outputs + VP = zeros(3,0); + + % Get resource grid dimensions + dims = lteDLResourceGridSize(ENB); + K = dims(1); + + % Virtual pilot cut-off values - indices in time or frequency at which + % the virtual pilot shall be discarded. As the symbol indices of the + % pilots passed in may be positive and negative the start and end of + % the effective grid varies therefore affecting the front and end pilot + % cut-off values + coTop = -6; % + coBottom = K+6; % + coFront = min(P_EST(2,:))-7; % floor(min(P_EST(2,:))/L)*L-5; % + coEnd = max(P_EST(2,:))+7; % ceil(max(P_EST(2,:))/L)*L+5; + + % Number of pilots in a subcarrier is dependent upon antenna +% noPilotinSC = (2-floor(ANTENNA/2))*nSF; + noPilotinSC = numel(P_EST(P_EST(1,:)==P_EST(1,1))); + + % Number of pilots in a symbol + noPilotinSym = K/6; + + % Sort by subcarrier + [~, indices] = sortrows(P_EST(1,:).'); + p_est_sorted_sc = P_EST(:,indices).'; + + if CREATE_TOP + % Repeat first and second SC containing pilots + p_est_rep_above = p_est_sorted_sc(1+noPilotinSC:noPilotinSC*2,:); + p_est_rep_above = [p_est_rep_above.' p_est_sorted_sc(1:noPilotinSC,:).'].'; + + % Subtract six to normalize subcarrier indices + p_est_rep_above(:,1) = p_est_rep_above(:,1) - 6; + + % Remove virtual pilots which are too far + p_est_rep_above(p_est_rep_above(:,1) < coTop,:) = []; + + % Initialize vp vector + VP_TOP = zeros(size(p_est_rep_above)); + + % Create virtual pilots on top + for p = 1:size(p_est_rep_above) + VP_TOP(p,:) = calculatePilot(p_est_rep_above(p,:), p_est_sorted_sc); + end + + % Transpose to make suitable for output + VP = [VP VP_TOP.']; + end + + if CREATE_BOTTOM + % Repeat last and second last SC containing pilots + p_est_rep_below = p_est_sorted_sc(end-noPilotinSC*2+1:end-noPilotinSC,:); + p_est_rep_below = [p_est_rep_below.' p_est_sorted_sc(end-noPilotinSC+1:end,:).'].'; + + % Add six to normalize subcarrier indices + p_est_rep_below(:,1) = p_est_rep_below(:,1) + 6; + + % Remove virtual pilots which are too far + p_est_rep_below(p_est_rep_below(:,1) > coBottom,:) = []; + + % Initialize vp vector + VP_BOTTOM = zeros(size(p_est_rep_below)); + + % Create virtual pilots on bottom + for p = 1:size(p_est_rep_below) + VP_BOTTOM(p,:) = calculatePilot(p_est_rep_below(p,:), p_est_sorted_sc); + end + + % Transpose to make suitable for output + VP = [VP VP_BOTTOM.']; + end + + % Sort by symbol + [~, indices] = sortrows(P_EST(2,:).'); + p_est_sorted_sym = P_EST(:,indices).'; + + if CREATE_FRONT + % Repeat first symbol containing pilots + p_est_rep_front = p_est_sorted_sym(noPilotinSym+1:(noPilotinSym*2),:); + + % Add subcarrier above and below + p_est_rep_front = addSCs(p_est_rep_front,noPilotinSym); + + % Repeat second symbol containing pilots + p_est_rep_front_second = p_est_sorted_sym(1:noPilotinSym,:); + + % Add subcarrier above and below + p_est_rep_front_second = addSCs(p_est_rep_front_second,noPilotinSym); + + % Concatenate to create virtual pilots + p_est_rep_front = [p_est_rep_front.' p_est_rep_front_second.'].'; + + % Subtract seven to normalize symbol indices + p_est_rep_front(:,2) = p_est_rep_front(:,2) - 7; + + % Remove virtual pilots which are too far + p_est_rep_front = removeTooFar(p_est_rep_front,p_est_sorted_sym,coTop,coBottom,coFront,coEnd); + + % Initialize vp vector + VP_FRONT = zeros(size(p_est_rep_front)); + + % Create virtual pilots on front + for p = 1:length(p_est_rep_front) + VP_FRONT(p,:) = calculatePilot(p_est_rep_front(p,:), p_est_sorted_sc); + end + + % Transpose for output + VP = [VP VP_FRONT.']; + end + + if CREATE_END + % Repeat last symbol containing pilots + p_est_rep_end = p_est_sorted_sym(end-(noPilotinSym*2)+1:end-noPilotinSym,:); + + % Add subcarrier above and below + p_est_rep_end = addSCs(p_est_rep_end,noPilotinSym); + + % Repeat second symbol containing pilots + p_est_rep_end_second = p_est_sorted_sym(end-noPilotinSym+1:end,:); + + % Add subcarrier above and below + p_est_rep_end_second = addSCs(p_est_rep_end_second,noPilotinSym); + + % Concatenate to create virtual pilots + p_est_rep_end = [p_est_rep_end.' p_est_rep_end_second.'].'; + + % Add seven to normalize symbol indices, mod controls case when + % only 2 pilots per subcarrier + p_est_rep_end(:,2) = p_est_rep_end(:,2) + 7; + + % Remove virtual pilots which are too far + p_est_rep_end = removeTooFar(p_est_rep_end,p_est_sorted_sym,coTop,coBottom,coFront,coEnd); + + % Initialize vp vector + VP_END = zeros(size(p_est_rep_end)); + + % Create virtual pilots on front + for p = 1:length(p_est_rep_end) + VP_END(p,:) = calculatePilot(p_est_rep_end(p,:), p_est_sorted_sc); + end + + % Transpose for output + VP = [VP VP_END.']; + end +end + +function vp = calculatePilot(p_est_rep, p_est_sorted) + % Calculate Euclidean distance between virtual pilot and other + % pilots + ind = sqrt((p_est_rep(1)-p_est_sorted(:,1)).^2+ (p_est_rep(2)-p_est_sorted(:,2)).^2); + % ind = (abs(p_est_rep(1)-p_est_sorted(:,1))+ abs(p_est_rep(2)-p_est_sorted(:,2))); + + % Sort from shortest to longest distance + [~, ind] = sortrows(ind); + + % Take three closest pilots + ind = ind(1:10); + pilots_use = p_est_sorted(ind,:); + + % If first 3 pilots used are in the same subcarrier or symbol then use 4th instead of 3rd + while ((pilots_use(3,1) == pilots_use(2,1)) && (pilots_use(2,1) == pilots_use(1,1)) || (pilots_use(3,2) == pilots_use(2,2)) && (pilots_use(2,2) == pilots_use(1,2))) + pilots_use(3,:) = []; + end + + % Calculate virtual pilot value + vp = calculateVirtualValue(pilots_use(3,:),pilots_use(1,:),pilots_use(2,:),p_est_rep(1),p_est_rep(2)); +end + +function new = calculateVirtualValue(a,b,c,xnew,ynew) + % Calculate vectors + AB = b-a; + AC = c-a; + + % Perform cross product + cro = cross(AB,AC); + + % Break out X, Y and Z plane coeffs + x = cro(1); + y = cro(2); + z = cro(3); + + % Calculate normal in equation + c = sum(cro.*a); + + % Calculate new z value + znew = (c - xnew*x - ynew*y)/z; + + % Return new point + new = [xnew ynew znew]; +end + +function pilots = addSCs(pilots, noPilotinSym) + % Add on extra subcarrier above + pilots(end+1,:) = pilots(1,:); + pilots(end,1) = pilots(end,1) - 6; + + % Add on extra subcarrier below + pilots(end+1,:) = pilots(noPilotinSym,:); + pilots(end,1) = pilots(end,1) + 6; +end + +function pilots = removeTooFar(pilots,p_est_sorted_sym,coTop,coBottom,coFront,coEnd) + pilots((pilots(:,2) < coFront) | (pilots(:,1) < coTop) | (pilots(:,2)> coEnd) | (pilots(:,1) > coBottom),:) = []; + + % Remove virtual pilots which fall within resource grid + removeInd = []; + for i = 1:size(pilots) + if (find(p_est_sorted_sym(:,2) == pilots(i,2))) + removeInd = [removeInd i]; %#ok + end + end + pilots(removeInd,:) = []; +end + +%CreateEdgeVirtualPilots +% Creates virtual pilots on the edges of the resource grid to improve +% interpolation results. Also adds virtual pilots around the current pilot +% estimates. +function vps = createEdgeVirtualPilots(enb,p_est) + + % Determine dimensions of current subframe + dims=lteDLResourceGridSize(enb); + K=dims(1); + L=dims(2); + + % Calculate virtual pilots beyond upper and lower bandwidth edge based + % on pilots present on subcarriers. Also create virtual pilots 1/2 RB + % above and below extent of current pilots. + vps=createDimensionVirtualPilots(p_est,2,unique([-6 min(p_est(1,:)-6) max(p_est(1,:)+6) K+6])); + + % Combine these frequency-direction VPs with original pilots. + temp = [p_est vps]; + + % Calculate virtual pilots beyond start and end of subframe based on + % pilots and frequency-direction VPs present in symbols. Also create + % virtual 1 OFDM symbol before and after extent of current pilots. + vps = [vps createDimensionVirtualPilots(temp,1,unique([-1 min(p_est(2,:)-1) max(p_est(2,:)+1) L+1]))]; + +end + +% dim=1 adds VPs in time +% dim=2 adds VPs in frequency +function vps = createDimensionVirtualPilots(p_est,dim,points) + + vps=[]; + + pilots = unique(p_est(dim,:)); + for i = pilots + x=find(p_est(dim,:)==i); + rep=1+double(length(x)==1); + adj=(1:rep)-rep; + pil = interp1(p_est(3-dim,x)+adj,repmat(p_est(3,x),1,rep),points,'linear','extrap'); + for j=1:length(points) + vps = [ vps [i;points(j);pil(j)] ]; %#ok + end + end + + if (dim==2) + vps = vps([2 1 3],:); + end + +end +