Python实现rician莱斯衰落和rician莱斯信道

谁践踏了优雅 2022-11-07 04:14 174阅读 0赞

# 基本概念 #

（1）莱斯分布  
在概率论中，Rice分布或Rician分布（或更不常见的是Ricean分布）是圆对称双变量正态随机变量的大小（可能具有非零均值（非中心））的概率分布。它以斯蒂芬·赖斯（Stephen O. Rice）的名字命名。

（2）莱斯信道  
Rician (rice)衰落信道一般有两种实现方式

*  第一种:将高斯白噪声通过线性时变滤波器获得,滤波器的传递函数由Droppler功率谱密度函数S(f)的平方根
 *  第二种:通过利用不同相位、幅度和频率的余弦波叠加来实现。以下代码就是该方法实现  
    ![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkzNTY5Ng_size_16_color_FFFFFF_t_70]  
    参考资料
 *  [Rician fading][]
 *  [Rice衰落信道的仿真-百度文库][Rice_-]

# Example #

来自[Github-Rician-Fading-Python][]

## 2.1 Rician衰落类 ##

# -*- coding: utf-8 -*-
    """ Created on Wed Nov 4 18:55:46 2020 @author: Jonathan Browning """
    
    import numpy as np
    from scipy.stats import gaussian_kde as kdf
    from scipy import special as sp
    
    class Rice:
        # numSamples = 2*(10**6) # the number of samples used in the simulation
        # numSamples = 64 #表示要产生的随机数的长度
        r = np.linspace(0, 6, 6000) # theoretical envelope PDF x axes
        theta = np.linspace(-np.pi, np.pi, 6000)    # theoretical phase PDF x axes
        
        def __init__(self, K, r_hat_2, phi,numSamples):
            
            # # user input checks and assigns value
            self.K = self.input_Check(K, "K", 0, 50)
            self.r_hat_2 = self.input_Check(r_hat_2, "\hat{r}^2", 0.5, 2.5)  
            self.phi = self.input_Check(phi, "\phi", -np.pi, np.pi)
            self.numSamples = numSamples
            # user input checks and assigns value
            
            # simulating and theri densities 
            self.multipathFading = self.complex_Multipath_Fading()
            self.xdataEnv, self.ydataEnv = self.envelope_Density()
            self.xdataPh, self.ydataPh = self.phase_Density()
            
            # theoretical PDFs calculated
            self.envelopeProbability = self.envelope_PDF()
            self.phaseProbability = self.phase_PDF()
    	# 输入格式检查
        def input_Check(self, data, inputName, lower, upper):
            # input_Check checks the user inputs
            
            # has a value been entered
            if data == "":
                raise ValueError(" ".join((inputName, "must have a numeric value")))
            
            # incase of an non-numeric input 
            try:
                data = float(data)
            except:    
                raise ValueError(" ".join((inputName, "must have a numeric value")))
        
            # data must be within the range
            if data < lower or data > upper:
                raise ValueError(" ".join((inputName, f"must be in the range [{lower:.2f}, {upper:.2f}]")))
            
            return data
    	
        def calculate_Means(self):
            # calculate_means calculates the means of the complex Gaussians representing the
            # in-phase and quadrature components
            
            p = np.sqrt(self.K * self.r_hat_2 / (1+self.K)) * np.cos(self.phi)
            q = np.sqrt(self.K * self.r_hat_2 / (1+self.K)) * np.sin(self.phi)
            
            return p, q
        
        def scattered_Component(self):
            # scattered_Component calculates the power of the scattered signal component
            
            sigma = np.sqrt(self.r_hat_2 / ( 2 * (1+self.K) ) )
            
            return sigma
        
        def generate_Gaussians(self, mean, sigma):
            # generate_Gaussians generates the Gaussian random variables
            
            gaussians = np.random.default_rng().normal(mean, sigma, self.numSamples)
            
            return gaussians
        
        def complex_Multipath_Fading(self):
            # complex_Multipath_Fading generates the complex fading random variables
            
            p, q = self.calculate_Means()
            sigma = self.scattered_Component()
            multipathFading = self.generate_Gaussians(p, sigma) + (1j*self.generate_Gaussians(q, sigma))
            
            return multipathFading
        
        def envelope_PDF(self):
            # envelope_PDF calculates the theoretical envelope PDF
            
            PDF = 2 * (1+self.K) * self.r / self.r_hat_2 \
                * np.exp(- self.K - ((1+self.K)*self.r**2)/self.r_hat_2) \
                * np.i0(2 * self.r * np.sqrt(self.K*(1+self.K)/self.r_hat_2))
                
            return PDF
    
        def phase_PDF(self):
            # phase_PDF calculates the theoretical phase PDF
            
            def q_func(x):
            # Q-function
            
                return 0.5-0.5*sp.erf(x/np.sqrt(2))     
            
            PDF = (1/(2*np.pi))* np.exp(- self.K) * (1 + (np.sqrt(4*np.pi*self.K) \
                  * np.exp(self.K * (np.cos(self.theta-self.phi))**2) *np.cos(self.theta-self.phi)) \
                  * (1 - q_func(np.sqrt(2*self.K) * np.cos(self.theta-self.phi))))
            
            return PDF
            
        
        def envelope_Density(self):
            # envelope_Density finds the envelope PDF of the simulated random variables
            
            R = np.sqrt((np.real(self.multipathFading))**2 + (np.imag(self.multipathFading))**2)
            kde = kdf(R)
            x = np.linspace(R.min(), R.max(), 100)
            p = kde(x)
            
            return x, p
        
        def phase_Density(self):
            # phase_Density finds the phase PDF of the simulated random variables
            
            R = np.angle(self.multipathFading)
            kde = kdf(R)
            x = np.linspace(R.min(), R.max(), 100)
            p = kde(x)
            
            return x, p

写一个Demo测试程序

from rice import *
    K =1 # K因子
    hat = 1 # 振幅
    phi = 1 # 相位
    numSamples = 16	# 随机样本的长度
    s = Rice(1,1,1,16)
    # rician衰落随机数
    print(s.multipathFading)

## 2.2 Ricican信道 ##

# signal输入的信号
    # channelResponse rician衰落随机数,就是以上代码求的s.multipathFading
    # SNRdb,信噪比
    def rician_channel(signal, channelResponse, SNRdb):
    	# 卷积
        convolved = np.convolve(signal, channelResponse)
        signal_power = np.mean(abs(convolved**2))
        sigma2 = signal_power * 10**(-SNRdb / 10)
        # 噪声
        noise = np.sqrt(sigma2 / 2) * (np.random.randn(*convolved.shape) + 1j * np.random.randn(*convolved.shape))
        return convolved + noise

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkzNTY5Ng_size_16_color_FFFFFF_t_70]: /images/20221023/64ef9bd91436480c9cdf2943009068b9.png
[Rician fading]: http://www.wirelesscommunication.nl/reference/chaptr03/ricepdf/rice.htm
[Rice_-]: https://wenku.baidu.com/view/c588e63d580216fc700afd15.html
[Github-Rician-Fading-Python]: https://github.com/Jonathan-Browning/Rician-Fading-Python