CN114514525A - 载波频率、初始相位、相位噪声的估计方法和相关设备 - Google Patents
载波频率、初始相位、相位噪声的估计方法和相关设备 Download PDFInfo
- Publication number
- CN114514525A CN114514525A CN202180005311.5A CN202180005311A CN114514525A CN 114514525 A CN114514525 A CN 114514525A CN 202180005311 A CN202180005311 A CN 202180005311A CN 114514525 A CN114514525 A CN 114514525A
- Authority
- CN
- China
- Prior art keywords
- time
- carrier
- covariance matrix
- phase
- nth
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/10—Noise analysis or noise optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computing Systems (AREA)
- Noise Elimination (AREA)
- Digital Transmission Methods That Use Modulated Carrier Waves (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本申请提供了一种载波频率、初始相位、相位噪声的估计方法和相关设备，其中方法包括：获取k+1个接收信号的幅度和相位，k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，k为大于或等于零的整数；根据k+1个接收信号的幅度确定第一协方差矩阵；根据第二向量和第二协方差矩阵确定第k时刻的载波频率和/或载波初始相位，第二向量为由k+1个接收信号的相位构成的向量，第二协方差矩阵为第一协方差矩阵与第三协方差矩阵的和。采用本申请，能够降低载波频率和/或载波初始相位的估计的计算量。
Description
技术领域
本申请涉及通信技术领域，具体涉及一种载波频率、初始相位、相位噪声的估计方法和相关设备。
背景技术
在通信、生物医学工程、雷达/声纳应用以及其他信号处理领域，如电网功率质量监测等，载波频率和相位估计都是一个经典而重要的问题。所有这些应用都可以借助于单一正弦载波的信号模型进行算法设计和性能分析。基于离散时间信号模型，在仅考虑加性高斯噪声(AWGN)时，多种频率、相位估计算法被提出并得到验证。目前最常用的方法是将频率和相位建模为未知的非随机参数，并应用极大似然(ML)估计理论。一个众所周知的解决方案是基于傅里叶变换的频域方法，即通过一维搜索来找到周期图的峰值，但这种方案的计算复杂度很大，计算量很大。
发明内容
本申请实施例提供了一种载波频率、载波初始相位的估计方法及相关设备，能够降低载波频率和/或载波初始相位的估计的计算量。
第一方面，本申请实施例提供了一种载波频率和/或载波初始相位的估计方法，包括：获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数；根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)；根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
在本申请实施例中，第一协方差矩阵∑∈是第一向量∈服从的高斯分布对应的协方差矩阵，第一向量∈是k+1个接收信号对应的加性观测相位噪声构成的向量，加性观测相位噪声等效于加性高斯噪声，第一协方差矩阵∑∈仅与接收信号的幅度有关；第二协方差矩阵∑r为第一协方差矩阵∑∈与第三协方差矩阵∑θ的和；第三协方差矩阵∑θ是由k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从的高斯分布对应的协方差矩阵，第三协方差矩阵∑θ是已知的；因此，在获取到k+1个接收信号的幅度和相位后，可以根据k+1个接收信号的幅度确定第一协方差矩阵∑∈，进而与第三协方差矩阵∑θ相加得到第二协方差矩阵∑r，再根据由k+1个接收信号的相位构成的第二向量∠r和第二协方差矩阵∑r即可确定出第k时刻的载波频率和/或载波初始相位，从而获取到k+1个接收信号的幅度和相位即可估计到载波频率和/或载波初始相位，降低了计算复杂度，也即降低了计算量。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波频率通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波频率，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述方法还包括：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波频率通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波初始相位通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波初始相位，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述方法还包括：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波初始相位通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
在一种可能的实现方式中，所述方法还包括：根据所述第二协方差矩阵∑r、所述第三协方差矩阵∑θ、所述第二向量∠r、所述第k时刻的载波频率和所述第k时刻的载波初始相位确定所述第k时刻的第二载波相位噪声。
在一种可能的实现方式中，k＝N-1，第N-1时刻的第二载波相位噪声通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的第二载波相位噪声，
表示第N-1时刻的载波频率，
表示第N-1时刻的载波初始相位，∑θ表示第N-1时刻的第三协方差矩阵，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，N＝[0，1，…，N-1]T，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述第N时刻的第二载波相位噪声通过以下公式确定：
其中：
其中：
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的第二载波相位噪声，
表示第N时刻的载波频率，
表示第N时刻的载波初始相位，∑θ (N)表示第N时刻的第三协方差矩阵，∑r (N)表示第N时刻的第二协方差矩阵，∠r(N)表示第N时刻的第二向量，N(N)＝[0，1，…，N-1，N]T，1(N)＝[1，1，…，1]T，∑r (N-1)表示第N-1时刻的第二协方差矩阵，N(N-1)＝[0，1，…，N-1]T，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(0)|表示第0时刻对应的接收信号的幅度，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
第二方面，本申请实施例提供了一种载波频率和/或载波初始相位的估计装置，包括：获取单元，用于获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数；确定单元，用于根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)；以及根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波频率通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波频率，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述获取单元还用于：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波频率通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波初始相位通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波初始相位，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述获取单元还用于：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波初始相位通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
在一种可能的实现方式中，所述确定单元还用于：根据所述第二协方差矩阵∑r、所述第三协方差矩阵∑θ、所述第二向量∠r、所述第k时刻的载波频率和所述第k时刻的载波初始相位确定所述第k时刻的第二载波相位噪声。
在一种可能的实现方式中，k＝N-1，第N-1时刻的第二载波相位噪声通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的第二载波相位噪声，
表示第N-1时刻的载波频率，
表示第N-1时刻的载波初始相位，∑θ表示第N-1时刻的第三协方差矩阵，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，N＝[0，1，…，N-1]T，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述第N时刻的第二载波相位噪声通过以下公式确定：
其中：
其中：
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的第二载波相位噪声，
表示第N时刻的载波频率，
表示第N时刻的载波初始相位，∑θ(N)表示第N时刻的第三协方差矩阵，∑r(N)表示第N时刻的第二协方差矩阵，∠r(N)表示第N时刻的第二向量，N(N)＝[0，1，…，N-1，N]T，1(N)＝[1，1，…，1]T，∑r (N-1)表示第N-1时刻的第二协方差矩阵，N(N-1)＝[0，1，…，N-1]T，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(0)|表示第0时刻对应的接收信号的幅度，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
本申请实施例第三方面提供了一种电子设备，包括处理器、存储器、通信接口，以及一个或多个程序，上述一个或多个程序被存储在上述存储器中，并且被配置由上述处理器执行，上述程序包括用于执行如上述第一方面中任一项所述的方法中的步骤的指令。
本申请实施例第四方面提供了一种芯片，其特征在于，包括：处理器，用于从存储器中调用并运行计算机程序，使得安装有上述芯片的设备执行如上述第一方面中任一项上述的方法。
本申请实施例第五方面提供了一种计算机可读存储介质，其存储用于电子数据交换的计算机程序，其中，上述计算机程序使得计算机执行如上述第一方面中任一项所述的方法。
本申请实施例第六方面提供了一种计算机程序产品，上述计算机程序产品使得计算机执行如上述第一方面中任一项所述的方法。
本申请的这些方面或其他方面在以下实施例的描述中会更加简明易懂。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本申请实施例提供的一种载波频率和/或载波初始相位的估计方法的流程示意图；
图2是本申请实施例提供的一几何向量的示意图；
图3是本申请实施例提供的一种载波相位噪声的估计方法的流程示意图；
图4是本申请实施例提供的一种载波频率和载波初始相位的迭代估计的方法流程示意图；
图5是本申请实施例提供的一种载波相位噪声的迭代估计的方法流程示意图；
图6是本申请实施例提供的一种载波频率和载波初始相位的估计的仿真性能图；
图7是本申请实施例提供的一种载波相位噪声的估计的仿真性能图；
图8是本申请实施例提供的另一种载波相位噪声的估计的仿真性能图；
图9是本申请实施例提供的一种载波频率和/或载波初始相位的估计装置的结构示意图；
图10是本申请实施例提供的一种电子设备的结构示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
本申请的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其他步骤或单元。
在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。
如前所述，现有方案采用频域估计方法，而现有的频率相位联合估计技术仅考虑脉冲噪声、高斯噪声，并没有考虑振荡器相位噪声的不利影响。本申请采用时域估计方法，将接收信号的相位作为观测数据样本输入估计器，利用极大似然估计理论来估计载波频率和相位。
众所周知，所有的天然和人造振荡器(无论是光学的、电子的、声学的、原子的或任何其他的)表现出相位和频率的不稳定性，统称为相位噪声。振荡器相位噪声是现代通信、雷达、光谱和计量系统的主要缺陷之一。此外，为了获取更高的频谱效率，超高速通信如今正在向更高的频段发展。使用的载波频率越高，所遇到的载波相位噪声就越大，系统损耗就越严重。例如，在实时测距和成像的雷达应用中，冗余相位噪声会引起时变偏移、杂散旁瓣、脉冲响应加宽以及雷达信号的低频相位调制，严重降低了雷达的检测和跟踪性能。因此，相位噪声是在实际应用中进行参数估计和性能分析时必须考虑的一个不可避免的因素。
考虑到无线、光学和雷达应用中实际振荡器的抖动行为，相位噪声通常是时变的。然而，为了简化理论分析，大多数算法都是基于准静态相位的假设设计的，没有考虑载波相位的时变特性，这导致估计性能显著下降。因此，针对时变相位噪声下的载波频率和相位估计问题，还没有系统的方法来实现频率和相位的联合估计，也没有确切的方法来估计这种情况下的相位噪声。
振荡器相位噪声是损害现代通信、雷达、光谱和计量等系统性能的主要因素之一。离散维纳相位噪声模型普遍适用于这些物理应用中的载波相位噪声，其中包括经常使用的半导体激光振荡器和无线射频振荡器。因此，本申请将解决在AWGN和维纳载波相位噪声的共同作用下，如何联合估计单一正弦信号的角度参数的问题。角度参数具体包含载波频率、载波相位和相位噪声三部分。我们将设计基于ML算法的未知载波频率和初始载波相位的联合时域估计方法，以及同时基于最大后验概率(MAP)算法的相位噪声的时域估计方法。并进一步给出具有低复杂度、逐样本迭代处理结构的估计方案，以确保在实际应用中可以实时迭代处理。
请参阅图1，图1是本申请实施例提供的一种载波频率和/或载波初始相位的估计方法的流程示意图，所述方法应用于电子设备，所述方法包括但不限于如下步骤：
步骤101、获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数。
其中，由于k＝0，1，2，3，......；故有第0时刻对应第1个接收信号，第1时刻对应第2个接收信号，第2时刻对应第3个接收信号，以此类推。
步骤102、根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)。
步骤103、根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
在本申请实施例中，第一协方差矩阵∑∈是第一向量∈服从的高斯分布对应的协方差矩阵，第一向量∈是k+1个接收信号对应的加性观测相位噪声构成的向量，加性观测相位噪声等效于加性高斯噪声，第一协方差矩阵∑∈仅与接收信号的幅度有关；第二协方差矩阵∑r为第一协方差矩阵∑∈与第三协方差矩阵∑θ的和；第三协方差矩阵∑θ是由k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从的高斯分布对应的协方差矩阵，第三协方差矩阵∑θ是已知的；因此，在获取到k+1个接收信号的幅度和相位后，可以根据k+1个接收信号的幅度确定第一协方差矩阵∑∈，进而与第三协方差矩阵∑θ相加得到第二协方差矩阵∑r，再根据由k+1个接收信号的相位构成的第二向量∠r和第二协方差矩阵∑r即可确定出第k时刻的载波频率和/或载波初始相位，从而获取到k+1个接收信号的幅度和相位即可估计到载波频率和/或载波初始相位，降低了计算复杂度，也即降低了计算量。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波频率通过以下公式确定：
在公式(1)中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波频率，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
作为一种示例，公式(1)可以通过以下方式得到：
在AWGN和维纳载波相位噪声的共同作用下，在第k时刻的单一正弦的接收信号时域表达式通常如公式(2)所示：
在公式(2)中，A是已知的发送信号幅度；
ω是未知的载波频率；
是未知的载波初始相位；所有的天然和人造振荡器(不论是光学、电子、声学、原子或其他)显示相位和频率的不稳定性统称为相位噪声，即θ(k)，θ(k)是在[-π，π)区间内随机游走的载波相位噪声，服从维纳过程：θ(k)＝θ(k-1)+Aθ(k)，其中，θ(0)＝0，{Δθ(k)}是一个独立、同分布、零均值、方差为
的实高斯随机变量序列；k为大于或者等于零的整数，例如，k＝0，1，2，3，......；{n(k)}则是一个离散时间、圆对称、零均值、协方差函数为E[n(k)n*(k-l)]＝N0δ(l)的复高斯随机变量序列。
进一步地，可以将接收信号改写为幅度-相位的极坐标形式，如公式(3)所示：
在公式(3)中，|r(k)|和∠r(k)分别是接收信号的幅度和相位；∈(k)是等效于AWGN的加性观测相位噪声，已证明在高信噪比下服从均值为0，方差为
的高斯分布，即
从图2(接收信号r的几何向量表示)中可以看出，∈(k)就是由AWGN所引起的相位角度的改变。
其中，可以将k＝0至N-1对应时间接收到的信号的相位用矢量形式表示，也即将第0时刻至第N-1时刻对应的接收信号的相位用矢量形式表示，如公式(4)所示：
在公式(4)中，第二向量∠r＝[∠r(0)，∠r(1)，…，∠r(N-1)]T，N＝[0，1，…，N-1]T，1＝[1，1，…，1]T，θ＝[θ(0)，θ(1)，…，θ(N-1)]T，第一向量∈＝[∈(0)，∈(1)，…，∈(N-1)]T。可以看出，以上向量均是N维列向量，且上标T表示向量的转置。θ和∈互相独立，且θ～N(0，∑θ)，∈～N(0，∑∈)，其中，第三协方差矩阵∑θ和第二协方差矩阵∑∈分别是：
和
最优估计算法推导主要是基于表达式(4)和一个理论事实：利用将接收振幅信息融入到加性观测相位噪声的接收信号相位∠r(k)的ML/MAP估计，等价于利用接收信号|r(k)|的ML/MAP估计；该方法可称之为基于相位的时域估计方法。
下面本申请将利用所有观测数据|r(k)|和∠r(k)(k＝0，1，2，3，…，N-1)对ω和
进行联合ML时域估计，并对θ进行MAP时域估计；其中，相应的估计值可以通过最大化联合概率密度函数
来获得，如公式(7)所示：
在公式(7)中，
和
分别代表k＝N-1时刻的载波频率和载波初始相位的ML估计值，
代表在k＝N-1时刻在载波相位噪声的MAP估计值。通过公式(4)可以很容易得出，在给定ω和
条件下，向量∠r服从联合高斯分布，也即
其中，协方差矩阵∑r是∑θ和∑∈的和。本申请通过解析公式(7)即可得出公式(1)。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波初始相位通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波初始相位，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
其中，本申请通过解析公式(7)即可得出公式(8)。
在一种可能的实现方式中，所述方法还包括：根据所述第二协方差矩阵∑r、所述第三协方差矩阵∑θ、所述第二向量∠r、所述第k时刻的载波频率和所述第k时刻的载波初始相位确定所述第k时刻的第二载波相位噪声。
在一种可能的实现方式中，k＝N-1，第N-1时刻的第二载波相位噪声通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的第二载波相位噪声，
表示第N-1时刻的载波频率，
表示第N-1时刻的载波初始相位，∑θ表示第N-1时刻的第三协方差矩阵，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，N＝[0，1，…，N-1]T，1＝[1，1，…，1]T。
其中，本申请通过解析公式(7)即可得出公式(9)。
应理解，通过公式(1)和(8)分别明确给出了载波频率和载波初始相位的联合ML估计值，将其代入(9)中可得到确定的载波相位噪声MAP估计值，也即得到
此外，在一示例中，对于无载波相位噪声的情况，也即
∑r＝∑∈，ML估计器，也即公式(1)可简化为公式(10)：
公式(8)可简化为公式(11)：
其中，公式(10)和(11)中的|r(k)|和∠r(k)分别为接收信号的幅度和相位。
本申请实施例，载波频率、载波初始相位以及第二载波相位噪声的估计值均可以随着样本数量N的增加迭代计算得到。下面，本申请示例性的给出一个简单的迭代实现方案。
首先，在k＝N-1时刻，令A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)，则公式(1)简化成公式(12)；其中，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝『1，1，…，1]T：
同理，公式(8)简化成公式(13)：
如此，随着N的增加，我们可以分别更新A(N-1)，B(N-1)，C(N-1)，D(N-1)和E(N-1)，继而实现估计值
和
的迭代计算。也即，在k＝N时刻，A(N)、B(N)、C(N)、D(N)和E(N)可以用A(N -1)、B(N-1)、C(N-1)、D(N-1)和E(N-1)迭代计算得到；将A(N)、B(N)、C(N)、D(N)和E(N)替换A(N-1)、B(N-1)、C(N-1)、D(N-1)和E(N-1)，分别代入公式(12)和公式(13)即可得到第N时刻的载波频率
和载波初始相位
因此，本申请当每计算一次载波频率或载波初始相位后，例如在第N-1时刻计算完一次载波频率或载波初始相位，仅需要存储对应的A(N-1)，B(N-1)，C(N-1)，D(N-1)和E(N-1)，在第N时刻获取到第N时刻对应的接收信号时，仅需要通过存储的A(N-1)，B(N-1)，C(N-1)，D(N-1)和E(N-1)计算得到A(N)、B(N)、C(N)、D(N)和E(N)，即可快速计算得到第N时刻的载波频率或载波初始相位。
在一种可能的实现方式中，所述方法还包括：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波频率通过以下公式确定：
其中，公式(14)中的A(N)、B(N)、C(N)、D(N)和E(N)与A(N-1)、B(N-1)、C(N-1)、D(N-1)和E(N-1)之间的迭代关系如公式(15)所示：
其中，在k＝1时刻，A(1)、B(1)、C(1)、D(1)和E(1)可以通过公式(16)直接推导得到：
公式(14)或公式(15)或公式(16)中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
应理解，该迭代过程充分利用了k＝N时刻新接收到的样本{|r(N)|，∠r(N)}，实现了实时性。上述迭代估计过程同样适用于无载波相位噪声的情况，即
在一种可能的实现方式中，所述方法还包括：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波初始相位通过以下公式确定：
同理，在公式(17)中的A(N)、B(N)、C(N)、D(N)和E(N)与A(N-1)、B(N-1)、C(N-1)、D(N-1)和E(N -1)之间的迭代关系如公式(15)所示；此外，在k＝1时刻，A(1)、B(1)、C(1)、D(1)和E(1)可以通过公式(16)直接推导得到。
应理解，本申请通过上述公式(14)和公式(17)计算得到第N时刻的
和
后.可以进一步计算得到第N时刻的第二载波相位噪声的估计值
在一种可能的实现方式中，所述第N时刻的第二载波相位噪声通过以下公式确定：
其中，公式(18)中的(∑r -1)(N)与(∑r -1)(N-1)的迭代关系如公式(19)所示：
在公式(19)中，
其中，在k＝1时刻，(∑r -1)(1)可通过公式(20)计算得到：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的第二载波相位噪声，
表示第N时刻的载波频率，
表示第N时刻的载波初始相位，∑θ(N)表示第N时刻的第三协方差矩阵，∑r (N)表示第N时刻的第二协方差矩阵，∠r(N)表示第N时刻的第二向量，N(N)＝[0，1，…，N-1，N]T，1(N)＝[1，1，…，1]T，∑r (N-1)表示第N-1时刻的第二协方差矩阵，N(N-1)＝[0，1，…，N-1]T，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(0)|表示第0时刻对应的接收信号的幅度，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
如此，在计算第二载波相位噪声时，(∑r -1)(N)可由(∑r -1)(N-1)递归得到，避免了复杂的矩阵求逆运算，从而减少了计算量。
需要说明的是，在本申请中，矩阵有上标(N)表示对应第N时刻的N+1维矩阵，矩阵有上标(N-1)表示第N-1时刻的N维矩阵；同理，向量有上标(N)表示对应第N时刻的N+1维向量，向量有上标(N-1)表示第N-1时刻的N维向量。
请参阅图3，图3是本申请实施例提供的一种载波相位噪声的估计方法的流程示意图，该方法应用于电子设备，具体应用于ML/MAP估计器，在相位噪声环境中，ML/MAP估计器的具体实现包括但不限于如下步骤：
步骤301、获取接收信号r(k)，k＝0，1，2，3，......，N-1；
步骤302、提取接收信号的幅度|r(k)|和相位∠r(k)；
步骤303、已知系统参数A、N0和
的条件下，采用公式(1)计算k＝N-1时刻的载波频率估计值；
步骤304、已知系统参数A、N0和
的条件下，采用公式(8)计算k＝N-1时刻的载波初始相位估计值；
步骤305、基于步骤303得到的载波频率估计值和步骤304得到的载波初始相位估计值，采用公式(9)计算得到k＝N-1时刻的载波相位噪声估计值。
请参阅图4，图4是本申请实施例提供的一种载波频率和载波初始相位的迭代估计方法的流程示意图，该方法应用于电子设备，具体应用于ML/MAP估计器，ML/MAP估计器采用迭代处理，包括但不限于如下步骤：
步骤401、获取接收信号r(k)，k＝0，1；
步骤402、提取接收信号的幅度和相位信息{|r(k)|，∠r(k)}；
步骤403、采用公式(16)计算A(1)、B(1)、C(1)、D(1)和E(1)；
步骤404、获取接收信号r(k)，并提取接收信号的幅度|r(k)|和相位∠r(k)，k＝2，3，…；
步骤405、获取存储的A(k-1)、B(k-1)、C(k-1)、D(k-1)和E(k-1)；
步骤406、采用公式(15)计算A(k)、B(k)、C(k)、D(k)和E(k)；
步骤407、采用公式(14)计算第k时刻的载波频率；
步骤408、采用公式(17)计算第k时刻的载波初始相位。
请参阅图5，图5是本申请实施例提供的一种载波相位噪声的迭代估计方法的流程示意图，该方法应用于电子设备，具体应用于ML/MAP估计器，ML/MAP估计器采用迭代处理，包括但不限于如下步骤：
步骤501、获取接收信号r(k)，k＝0，1，2，3，…；
步骤502、提取接收信号的幅度和相位信息{|r(k)|，∠r(k)}；
步骤503、k＝1时刻，采用公式(20)计算∑r -1；
步骤504、k-1时刻，采用公式(19)计算∑r -1；
步骤505、获取k时刻的载波频率估计值和载波初始相位估计值；
步骤506、采用公式(18)计算k时刻的载波相位噪声估计值。
示例性的，若给定载波频率0.05和初始相位0.25π，在不同相位噪声环境下，用蒙特卡洛仿真计算ML/MAP估计的逆均方误差(IMSE)，验证估计精度。以载波频率估计为例，即是计算：
其中，在公式(21)中，q＝105以保证仿真精度。
图6为本申请实施例提供一种在不同信噪比下ML估计载波频率和载波初始相位的IMSE和逆克拉美罗下限(inverseCRLB，ICRLB)的仿真性能，其中，图6所示的ML估计ω和
的仿真性能中，N＝16；考虑了
和10-2rad2两种情况，可以看出在约0dB信噪比时，估计器(1)和(8)，也即公式(1)和(8)的均方误差(MSE)性能就逼近克拉美罗下限(CRLB)，有效验证了ML时域估计方法的高精度。
图7和图8为本申请实施例提供一种在不同信噪比下MAP估计载波相位噪声{θ(k)，k＝1∶15}的IMSE和逆贝叶斯-克拉美罗下限(BCRLB)的仿真性能，其中，N＝16，图7考虑了
的情况，图8考虑了
的情况；可以看出，在低相位噪声
下，在约0dB信噪比时，估计器(9)，也即公式(9)的MSE性能就可达到BCRLB；而在
条件下，在约5dB信噪比时才能达到可接受的估计精度，即IMSE与逆BCRLB有约1dB的偏差。
本申请解决了在时变相位噪声影响下，单一正弦的载波频率和载波初始相位以及载波相位噪声的最优的实时估计问题，可以在较低的信噪比下获得较高的估计性能。估计器都表示为接收信号相位的加权线性组合，在实践中易于迭代实现。
在约0dB信噪比情况下，估计器(1)和(8)可以达到作为ML估计精度基准的克拉美罗下限(CRLB)，估计器(9)可以达到作为MAP估计精度基准的贝叶斯-克拉美罗下限(BCRLB)。估计器(1)和(8)同样适用于非相位噪声(纯AWGN)环境。
本申请与其他现有估计方法相比考虑了相位噪声的影响，具有更好的估计性能；时域估计复杂度低，每一次计算复杂度为O(1)，经过N次迭代的计算复杂度为o(N)。
本申请适用于通信、生物医学工程、雷达/声纳应用以及其他信号处理领域，如电网功率质量监测等；本申请方法不涉及复杂算法，便于硬件实现。
上述详细阐述了本申请实施例的方法，下面提供了本申请实施例的装置。
请参见图9，图9是本申请实施例提供的一种载波频率和/或载波初始相位的估计装置900的结构示意图，应用于电子设备，该估计装置900可以包括获取单元901和确定单元902，其中，各个单元的详细描述如下：
获取单元901，用于获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数；
确定单元902，用于根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)；以及根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波频率通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波频率，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述获取单元901还用于：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波频率通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N -1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
在一种可能的实现方式中，k＝N-1，所述N为正整数，第N-1时刻的载波初始相位通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波初始相位，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述获取单元901还用于：获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；所述第N时刻的载波初始相位通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
在一种可能的实现方式中，所述确定单元902还用于：根据所述第二协方差矩阵∑r、所述第三协方差矩阵∑θ、所述第二向量∠r、所述第k时刻的载波频率和所述第k时刻的载波初始相位确定所述第k时刻的第二载波相位噪声。
在一种可能的实现方式中，k＝N-1，第N-1时刻的第二载波相位噪声通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的第二载波相位噪声，
表示第N-1时刻的载波频率，
表示第N-1时刻的载波初始相位，∑θ表示第N-1时刻的第三协方差矩阵，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，N＝[0，1，…，N-1]T，1＝[1，1，…，1]T。
在一种可能的实现方式中，所述第N时刻的第二载波相位噪声通过以下公式确定：
其中：
其中：
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的第二载波相位噪声，
表示第N时刻的载波频率，
表示第N时刻的载波初始相位，∑θ (N)表示第N时刻的第三协方差矩阵，∑r (N)表示第N时刻的第二协方差矩阵，∠r(N)表示第N时刻的第二向量，N(N)＝[0，1，…，N-1，N]T，1(N)＝[1，1，…，1]T，∑r (N-1)表示第N-1时刻的第二协方差矩阵，N(N-1)＝[0，1，…，N-1]T，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(0)|表示第0时刻对应的接收信号的幅度，|r(1)|表示第1时刻对应的接收信号的幅度；
表示一方差，第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
需要说明的是，各个单元的实现还可以对应参照图1至图8所示的实施例的相应描述。当然，本申请实施例提供的估计装置900包括但不限于上述单元模块，例如：该估计装置900还可以包括存储单元903，存储单元903可以用于存储该估计装置900的程序代码和数据。图9所描述的估计装置900带来的有益效果可参照前述实施例的描述，此处不在重复描述。
请参见图10，图10是本申请实施例提供的一种电子设备1010的结构示意图，该电子设备1010包括处理器1011、存储器1012和通信接口1013，上述处理器1011、存储器1012和通信接口1013通过总线1014相互连接。
存储器1012包括但不限于是随机存储记忆体(random access memory，RAM)、只读存储器(read-only memory，ROM)、可擦除可编程只读存储器(erasable programmableread only memory，EPROM)、或便携式只读存储器(compact disc read-only memory，CD-ROM)，该存储器1012用于相关计算机程序及数据。通信接口1013用于接收和发送数据。
处理器1011可以是一个或多个中央处理器(central processing unit，CPU)，在处理器1011是一个CPU的情况下，该CPU可以是单核CPU，也可以是多核CPU。
该电子设备1010中的处理器1011用于读取上述存储器1012中存储的计算机程序代码，执行以下操作：获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数；根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)；根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
需要说明的是，各个操作的实现还可以对应参照图1至图8所示的实施例的相应描述。图10所描述的电子设备1010带来的有益效果可参照前述实施例的描述，此处不在重复描述。
本申请实施例还提供一种芯片，上述芯片包括至少一个处理器，存储器和接口电路，上述存储器、上述收发器和上述至少一个处理器通过线路互联，上述至少一个存储器中存储有计算机程序；上述计算机程序被上述处理器执行时，图1或图3至图5所示的方法流程得以实现。
本申请实施例还提供一种计算机可读存储介质，上述计算机可读存储介质中存储有计算机程序，当其在电子设备上运行时，图1或图3至图5所示的方法流程得以实现。
本申请实施例还提供一种计算机程序产品，当上述计算机程序产品在电子设备上运行时，图1或图3至图5所示的方法流程得以实现。
应理解，本申请实施例中提及的处理器可以是中央处理单元(CentralProcessing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital SignalProcessor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
还应理解，本申请实施例中提及的存储器可以是易失性存储器或非易失性存储器，或可包括易失性和非易失性存储器两者。其中，非易失性存储器可以是只读存储器(Read-Only Memory，ROM)、可编程只读存储器(Programmable ROM，PROM)、可擦除可编程只读存储器(Erasable PROM，EPROM)、电可擦除可编程只读存储器(Electrically EPROM，EEPROM)或闪存。易失性存储器可以是随机存取存储器(Random Access Memory，RAM)，其用作外部高速缓存。通过示例性但不是限制性说明，许多形式的RAM可用，例如静态随机存取存储器(Static RAM，SRAM)、动态随机存取存储器(Dynamic RAM，DRAM)、同步动态随机存取存储器(Synchronous DRAM，SDRAM)、双倍数据速率同步动态随机存取存储器(Double DataRate SDRAM，DDR SDRAM)、增强型同步动态随机存取存储器(Enhanced SDRAM，ESDRAM)、同步连接动态随机存取存储器(Synchlink DRAM，SLDRAM)和直接内存总线随机存取存储器(Direct Rambus RAM，DR RAM)。
需要说明的是，当处理器为通用处理器、DSP、ASIC、FPGA或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件时，存储器(存储模块)集成在处理器中。
应注意，本文描述的存储器旨在包括但不限于这些和任意其它适合类型的存储器。
还应理解，本文中涉及的第一、第二、第三、第四以及各种数字编号仅为描述方便进行的区分，并不用来限制本申请的范围。
应理解，本文中术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
应理解，在本申请的各种实施例中，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实施过程构成任何限定。
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，上述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。
上述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本申请各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。
上述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本申请各个实施例所示方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(Read-Only Memory，ROM)、随机存取存储器(Random Access Memory，RAM)、磁碟或者光盘等各种可以存储程序代码的介质。
本申请实施例方法中的步骤可以根据实际需要进行顺序调整、合并和删减。
本申请实施例装置中的模块可以根据实际需要进行合并、划分和删减。
以上所述，以上实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。
Claims (11)
1.一种载波频率和/或载波初始相位的估计方法，其特征在于，包括：
获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数；
根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)；
根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
2.根据权利要求1所述的方法，其特征在于，k＝N-1，所述N为正整数，第N-1时刻的载波频率通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波频率，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
3.根据权利要求2所述的方法，其特征在于，所述方法还包括：
获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；
所述第N时刻的载波频率通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
4.根据权利要求1所述的方法，其特征在于，k＝N-1，所述N为正整数，第N-1时刻的载波初始相位通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的载波初始相位，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T。
5.根据权利要求4所述的方法，其特征在于，所述方法还包括：
获取第N时刻对应的接收信号的幅度和相位，所述第N时刻对应的接收信号为在第N时刻接收到的信号；
所述第N时刻的载波初始相位通过以下公式确定：
其中：
其中，A(N-1)＝(NT∑r -1∠r)(N-1)，B(N-1)＝(1T∑r -11)(N-1)，C(N-1)＝(1T∑r -1N)(N-1)，D(N-1)＝(1T∑r -1∠r)(N-1)，E(N-1)＝(NT∑r -1N)(N-1)；
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的载波频率，∠r(N)表示第N时刻对应的接收信号的相位，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(N)|表示第N时刻对应的接收信号的幅度，N＝[0，1，…，N-1]T，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，1＝[1，1，…，1]T，∠r(1)表示第1时刻对应的接收信号的相位，|r(1)|表示第1时刻对应的接收信号的幅度；第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
6.根据权利要求2-5任一项所述的方法，其特征在于，所述方法还包括：
根据所述第二协方差矩阵∑r、所述第三协方差矩阵∑θ、所述第二向量∠r、所述第k时刻的载波频率和所述第k时刻的载波初始相位确定所述第k时刻的第二载波相位噪声。
7.根据权利要求6所述的方法，其特征在于，k＝N-1，第N-1时刻的第二载波相位噪声通过以下公式确定：
其中，上标(N-1)表示第N-1时刻，
表示第N-1时刻的第二载波相位噪声，
表示第N-1时刻的载波频率，
表示第N-1时刻的载波初始相位，∑θ表示第N-1时刻的第三协方差矩阵，∑r表示第N-1时刻的第二协方差矩阵，∠r表示第N-1时刻的第二向量，N＝[0，1，…，N-1]T，1＝[1，1，…，1]T。
8.根据权利要求7所述的方法，其特征在于，所述第N时刻的第二载波相位噪声通过以下公式确定：
其中：
其中：
其中：
其中，上标(N)表示第N时刻，上标(N-1)表示第N-1时刻，上标(1)表示第1时刻，
表示第N时刻的第二载波相位噪声，
表示第N时刻的载波频率，
表示第N时刻的载波初始相位，∑θ (N)表示第N时刻的第三协方差矩阵，∑r (N)表示第N时刻的第二协方差矩阵，∠r(N)表示第N时刻的第二向量，N(N)＝[0，1，…，N-1，N]T，1(N)＝[1，1，…，1]T，∑r (N-1)表示第N-1时刻的第二协方差矩阵，N(N-1)＝[0，1，…，N-1]T，N0表示白噪声的单边带功率谱密度，A表示发送信号的幅度，|r(0)|表示第0时刻对应的接收信号的幅度，|r(1)|表示第1时刻对应的接收信号的幅度；第一载波相位噪声θ(k)＝θ(k-1)+Δθ(k)，θ(0)＝0，{Δθ(k)}是一个方差为
的高斯随机变量序列。
9.一种载波频率和/或载波初始相位的估计装置，其特征在于，包括：
获取单元，用于获取k+1个接收信号的幅度和相位，所述k+1个接收信号为在连续的k+1个时刻分别接收到的k+1个信号，所述k为大于或等于零的整数；
确定单元，用于根据所述k+1个接收信号的幅度确定第一协方差矩阵∑∈，其中，由所述k+1个接收信号对应的加性观测相位噪声构成的第一向量∈服从以下高斯分布：∈～N(0，∑∈)；
以及根据第二向量∠r和第二协方差矩阵∑r确定第k时刻的载波频率和/或载波初始相位，其中，所述第二向量∠r为由所述k+1个接收信号的相位构成的向量，所述第二协方差矩阵∑r为所述第一协方差矩阵∑∈与第三协方差矩阵∑θ的和，由所述k+1个接收信号对应的随机游走的第一载波相位噪声构成的第三向量θ服从以下高斯分布：θ～N(0，∑θ)。
10.一种电子设备，其特征在于，包括处理器、存储器、通信接口，以及一个或多个程序，所述一个或多个程序被存储在所述存储器中，并且被配置由所述处理器执行，所述程序包括用于执行权利要求1-8任一项所述的方法中的步骤的指令。
11.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行以实现权利要求1-8任一项所述的方法。
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| PCT/CN2021/099179 WO2022257046A1 (zh) | 2021-06-09 | 2021-06-09 | 载波频率、初始相位、相位噪声的估计方法和相关设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114514525A true CN114514525A (zh) | 2022-05-17 |
| CN114514525B CN114514525B (zh) | 2022-10-25 |
Family
ID=81550698
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202180005311.5A Active CN114514525B (zh) | 2021-06-09 | 2021-06-09 | 载波频率、初始相位、相位噪声的估计方法和相关设备 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN114514525B (zh) |
| WO (1) | WO2022257046A1 (zh) |
Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20020159539A1 (en) * | 2001-02-15 | 2002-10-31 | Alcock Robert David | Method of estimating the carrier frequency of a phase-modulated signal |
| US7302016B1 (en) * | 2003-08-13 | 2007-11-27 | L-3 Communicaations Corporation | Phase estimator with bias correction |
| CN102656824A (zh) * | 2009-12-18 | 2012-09-05 | 阿尔卡特朗讯 | 用于相干光系统中非线性损害监视和减轻的载波相位估计器 |
| CN103095615A (zh) * | 2011-10-31 | 2013-05-08 | 华为技术有限公司 | 一种相位噪声估计、补偿的方法及装置 |
| US20140086373A1 (en) * | 2012-09-27 | 2014-03-27 | Scott Kenealy | Method of Channel Estimation and a Channel Estimator |
| WO2014100699A1 (en) * | 2012-12-20 | 2014-06-26 | Qualcomm Incorporated | Systems and methods to mitigate phase noise |
| US9559417B1 (en) * | 2010-10-29 | 2017-01-31 | The Boeing Company | Signal processing |
| CN106646543A (zh) * | 2016-12-22 | 2017-05-10 | 成都正扬博创电子技术有限公司 | 基于主从式aukf算法的高动态卫星导航信号载波追踪方法 |
| CN106788653A (zh) * | 2016-11-25 | 2017-05-31 | 西安电子科技大学 | 一种基于协方差矩阵重构的自适应波束形成方法 |
| JP2017138260A (ja) * | 2016-02-05 | 2017-08-10 | 沖電気工業株式会社 | 整相器および整相処理方法 |
| CN110231638A (zh) * | 2019-06-28 | 2019-09-13 | 北京邮电大学 | 一种无线网络下的载波相位定位方法及装置 |
| CN111007557A (zh) * | 2019-12-10 | 2020-04-14 | 中国矿业大学 | 自适应运动学模型辅助的gnss载波相位与多普勒融合测速方法 |
| CN111537966A (zh) * | 2020-04-28 | 2020-08-14 | 东南大学 | 一种适用于毫米波车载雷达领域的阵列天线误差校正方法 |
| EP3723332A1 (en) * | 2019-04-08 | 2020-10-14 | Mitsubishi Electric R&D Centre Europe B.V. | Circular pilot sequences for joint channel and phase noise estimation |
-
2021
- 2021-06-09 CN CN202180005311.5A patent/CN114514525B/zh active Active
- 2021-06-09 WO PCT/CN2021/099179 patent/WO2022257046A1/zh active Application Filing
Patent Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20020159539A1 (en) * | 2001-02-15 | 2002-10-31 | Alcock Robert David | Method of estimating the carrier frequency of a phase-modulated signal |
| US7302016B1 (en) * | 2003-08-13 | 2007-11-27 | L-3 Communicaations Corporation | Phase estimator with bias correction |
| CN102656824A (zh) * | 2009-12-18 | 2012-09-05 | 阿尔卡特朗讯 | 用于相干光系统中非线性损害监视和减轻的载波相位估计器 |
| US9559417B1 (en) * | 2010-10-29 | 2017-01-31 | The Boeing Company | Signal processing |
| CN103095615A (zh) * | 2011-10-31 | 2013-05-08 | 华为技术有限公司 | 一种相位噪声估计、补偿的方法及装置 |
| US20140086373A1 (en) * | 2012-09-27 | 2014-03-27 | Scott Kenealy | Method of Channel Estimation and a Channel Estimator |
| WO2014100699A1 (en) * | 2012-12-20 | 2014-06-26 | Qualcomm Incorporated | Systems and methods to mitigate phase noise |
| JP2017138260A (ja) * | 2016-02-05 | 2017-08-10 | 沖電気工業株式会社 | 整相器および整相処理方法 |
| CN106788653A (zh) * | 2016-11-25 | 2017-05-31 | 西安电子科技大学 | 一种基于协方差矩阵重构的自适应波束形成方法 |
| CN106646543A (zh) * | 2016-12-22 | 2017-05-10 | 成都正扬博创电子技术有限公司 | 基于主从式aukf算法的高动态卫星导航信号载波追踪方法 |
| EP3723332A1 (en) * | 2019-04-08 | 2020-10-14 | Mitsubishi Electric R&D Centre Europe B.V. | Circular pilot sequences for joint channel and phase noise estimation |
| CN110231638A (zh) * | 2019-06-28 | 2019-09-13 | 北京邮电大学 | 一种无线网络下的载波相位定位方法及装置 |
| CN111007557A (zh) * | 2019-12-10 | 2020-04-14 | 中国矿业大学 | 自适应运动学模型辅助的gnss载波相位与多普勒融合测速方法 |
| CN111537966A (zh) * | 2020-04-28 | 2020-08-14 | 东南大学 | 一种适用于毫米波车载雷达领域的阵列天线误差校正方法 |
Non-Patent Citations (5)
| Title |
|---|
| 侯永宏: "数字电视地面传输中的关键技术研究", 《中国优秀博硕士学位论文全文数据库(博士)信息科技辑(月刊)》 * |
| 王晓峰等: "用于正弦波频率估计的修正Kay算法", 《电讯技术》 * |
| 王芳等: "一种快速高精度的改进Fitz频率估计算法", 《江西师范大学学报(自然科学版)》 * |
| 郭履翔: "基于DFT的实正弦信号频率估计算法的研究", 《中国优秀博硕士学位论文全文数据库(硕士)信息科技辑(月刊)》 * |
| 顾建峰等: "基于伪协方差矩阵的频率和角度联合估计算法", 《通信学报》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114514525B (zh) | 2022-10-25 |
| WO2022257046A1 (zh) | 2022-12-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Fu et al. | MAP/ML estimation of the frequency and phase of a single sinusoid in noise | |
| O'Shea | On refining polynomial phase signal parameter estimates | |
| KR102458095B1 (ko) | 위상 교정 방법 및 장치 | |
| CN104142425B (zh) | 一种正弦信号频率估计的相位匹配方法 | |
| Wang et al. | Joint ML/MAP estimation of the frequency and phase of a single sinusoid with Wiener carrier phase noise | |
| Omelchuk et al. | A closed-form ARMA-based ML-estimator of a single-tone frequency | |
| Franz et al. | Frequency estimation using warped discrete Fourier transform | |
| US20120253742A1 (en) | Method and System for Determining Parameters of Sinusoidal Signals | |
| CN114514525B (zh) | 载波频率、初始相位、相位噪声的估计方法和相关设备 | |
| Zhou et al. | Efficient and accurate frequency estimator under low SNR by phase unwrapping | |
| US11196609B2 (en) | Interpolation method for polar signals in radio frequency transmitters | |
| CN114755629B (zh) | 提高室内定位精度的方法及装置 | |
| CN113726416B (zh) | 一种卫星通信载波同步方法、装置及通信设备 | |
| WO2022127476A1 (zh) | 谐波消除方法及装置、存储介质、终端 | |
| Murakami et al. | An analytical solution to Jacobsen estimator for windowed signals | |
| Elasmi-Ksibi et al. | Closed-form real single-tone frequency estimator based on a normalized IIR notch filter | |
| Zhou et al. | Improved frequency estimation algorithm by least squares phase unwrapping | |
| Selva | Efficient maximum likelihood estimation of a 2-d complex sinusoidal based on barycentric interpolation | |
| Hussain et al. | A fast hybrid DSC-GS-MLE approach for multiple sinusoids estimation | |
| Maslakov | Choice of regularization parameter in adaptive filtering problems | |
| Nizampatnam et al. | Signal Frequency Estimation via Kalman Filter and Least Squares Approach for Non-uniform Signals | |
| Li et al. | IFR estimation of overlapped multicomponent signals from time frequency rate distribution | |
| Petrović | New procedure for estimation of amplitude and phase of analog multiharmonic signal based on the differential irregular samples | |
| Ou et al. | A refined PCPF algorithm for estimating the parameters of multicomponent polynomial-phase signals | |
| Boccato et al. | Improving the efficiency of natural computing algorithms in DOA estimation using a noise filtering approach |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |