python负指数信号的生成和波形变换

功能：

实现负指数信号的参数可调
高斯白噪声的均值和方差可调 椒盐噪声的生成概率和幅度可调
波形变换为 负指数->矩形->锯齿->梯形

复制 安装numpy，matplotlib，scipy安装包即可用

import numpy as np
import matplotlib.pyplot as plt
import random
import math
from scipy.integrate import simps

'''
    descent：负指数曲线的陡峭程度  
    length_half:负指数一个周期的半径
    resolution_rate：分辨率/采样间隔 
    A：信号幅度
    Duty_Ratio：矩形脉冲占空比
    Num_T：周期个数

    mean：高斯白噪声的均值
    bias：高斯白噪声的方差

    Sp_prob：椒盐噪声的生成概率
    Sp_A：椒盐噪声的尖刺幅度
'''

descent = 1
length_half = 5
resolution_rate = 0.01
A = 50
Duty_Ratio = 0.001
Num_T = 10

mean = 0
bias = 1

Sp_prob = 0.001
Sp_A = 200

def exp_comptue(x,A):
    signal1 = np.exp(-1*descent*x)*A
    return signal1

def Sp_Noise_Generate(Sp_A,Sp_prob):
    a = random.random()
    if(a<Sp_prob):
        return -1*Sp_A
    elif a>(1-Sp_prob):
        return Sp_A
    else :
        return 0

def SNR(array1,array2):
    SNR = 10*math.log(exp_comptue_energy(array1)/exp_comptue_energy(array2),10)
    return SNR

#生成负指数周期信号
x_generate = np.linspace(0,2*length_half*Num_T,((int((2*length_half)/resolution_rate))+1)*Num_T)
x_generate_temp = np.linspace(0,2*length_half,(int((2*length_half)/resolution_rate))+1)

y_generate_temp = exp_comptue(x_generate_temp,A)
y_generate = y_generate_temp
for i in range(1,Num_T):
    y_generate = np.hstack([y_generate,y_generate_temp])

#生成高斯白噪声
White_Gaussian_Noise_y = np.random.normal(mean,bias,((int((2*length_half)/resolution_rate))+1)*Num_T)

#生成椒盐噪声
Sp_Noise = [Sp_Noise_Generate(Sp_A,Sp_prob)]
for i in range(1,((int((2*length_half)/resolution_rate))+1)*Num_T):
    Sp_Noise.append(Sp_Noise_Generate(Sp_A,Sp_prob))
Sp_Noise=np.array(Sp_Noise)

#负指数周期信号变成矩形脉冲周期信号
def TOrectangle(array):
    array = np.where(array<A*Duty_Ratio,array,1*A)
    array = np.where(array>A*Duty_Ratio,array,0)
    return array

#矩形脉冲周期信号积分变成三角波
def TOtriangle(array):
    array = list(array)
    a = []
    b = []
    pad = 0
    for x in range(1,Num_T+1):
        for i in range(1,int((2*length_half)/resolution_rate+2)):
            for j in range(1,i+1):
                pad = resolution_rate*array[j-1]+pad
            #pad = array[i] + pad
            a.append(pad)
            pad = 0
        b = np.hstack([b,a])
        a = []
    ladder = b
    k = simps(array,x_generate)/Num_T
    #print(k)
    b = np.where(b<=k,b,0)

    parameters = {'ladder'  :ladder,
                  'triangle':b,
                  'simps'   :k
    }

    return parameters

#三角波信号生成梯形波信号
def TOladder(parameters):
    ladder = list(parameters.get('ladder'))
    k = parameters.get('simps')
    for i in range(1,Num_T):
        for j in range(((int((2*length_half)/resolution_rate))+1)*i,((int((2*length_half)/resolution_rate))+1)*Num_T):
            ladder[j] = k - ladder[j]
    ladder = np.array(ladder)
    return ladder

#计算信号能量
def exp_comptue_energy(array):
    array = list(np.square(array))
    energy = 0
    for i in range(1,int((2*length_half)/resolution_rate+1)*10+1):
        energy = resolution_rate*array[i-1]+energy
    return energy

B = TOladder(TOtriangle(TOrectangle(y_generate)))
print(SNR(B,White_Gaussian_Noise_y))
print(exp_comptue_energy(TOladder(TOtriangle(TOrectangle(y_generate)))))


#幅度从0到100  白噪声不变  添加高斯白噪声周期负指数信号的 信噪比  SNR_draw(x_generate,White_Gaussian_Noise_y)
def SNR_draw(array1,array2):
    SNR_draw = []
    for A in range(1,((int((2*length_half)/resolution_rate))+1)*Num_T+1):
        SNR_draw.append(SNR(exp_comptue(array1,A),array2))
    return SNR_draw


plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


plt.figure(1)
plt.subplot(2,2,1) 
plt.title('负指数信号')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate_temp,y_generate_temp)
plt.subplot(2,2,2) 
plt.title('负指数周期信号')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,y_generate)
plt.subplot(2,2,3) 
plt.title('负指数周期信号添加高斯噪声')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,y_generate+White_Gaussian_Noise_y)
plt.subplot(2,2,4) 
plt.title('负指数周期信号添加椒盐噪声')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,y_generate+Sp_Noise)
plt.show()

plt.figure(2)

plt.subplot(2,2,1) 
plt.title('脉冲信号')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOrectangle(y_generate))

plt.subplot(2,2,2)
plt.title('锯齿波')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOtriangle(TOrectangle(y_generate)).get('triangle'))

plt.subplot(2,2,3)
plt.title('梯形波')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOladder((TOtriangle(TOrectangle(y_generate)))))

plt.subplot(2,2,4)
plt.title('噪声信号')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,White_Gaussian_Noise_y)
plt.plot(x_generate,exp_comptue(x_generate,A)+White_Gaussian_Noise_y)

plt.show()

plt.figure(3)
plt.subplot(1,2,1)
plt.title('高斯白噪声')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,White_Gaussian_Noise_y)
plt.subplot(1,2,2)
plt.title('椒盐噪声')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,Sp_Noise)
plt.show()


plt.figure(4)
plt.subplot(1,2,1)
plt.title('添加高斯白噪声的负指数信号转化成矩形波')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOrectangle(y_generate)+White_Gaussian_Noise_y)
plt.subplot(1,2,2)
plt.title('添加椒盐噪声的负指数信号转化成矩形波')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOrectangle(y_generate+Sp_Noise))
plt.show()

plt.figure(5)
plt.subplot(1,2,1)
plt.title('添加高斯白噪声的负指数信号转化成三角波')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOtriangle(TOrectangle(y_generate)+White_Gaussian_Noise_y).get('triangle'))
plt.subplot(1,2,2)
plt.title('添加椒盐噪声的负指数信号转化成三角波')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,TOtriangle(TOrectangle(y_generate+Sp_Noise)).get('triangle'))
plt.show()

plt.figure(6)
plt.subplot(1,2,1)
plt.title('添加高斯白噪声的负指数信号 A = 50')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.plot(x_generate,y_generate+White_Gaussian_Noise_y)
plt.subplot(1,2,2)
plt.title('改变负指数信号的幅度计算信噪比')
plt.xlabel('Amplitude')
plt.ylabel('10*log(S/N)  (db)')
plt.plot(x_generate,SNR_draw(x_generate,White_Gaussian_Noise_y))
plt.show()