python分析PI特性曲线（2）

python分析PI特性曲线（2）

# -*- coding: utf-8 -*-

import pywt
import numpy as np
import matplotlib.pyplot as plt

import pandas as pd
from PIL.ImageQt import rgb
from matplotlib.collections import LineCollection
from numpy import double


def color_map(data, cmap):
    """数值映射为颜色"""
    dmin, dmax = np.nanmin(data), np.nanmax(data)
    cmo = plt.cm.get_cmap(cmap)
    cs, k = list(), 256 / cmo.N
    for i in range(cmo.N):

        c = cmo(i)
        # print(c)
        for j in range(int(i * k), int((i + 1) * k)):
            cs.append(c)
    cs = np.array(cs)
    data = np.uint8(255 * (data - dmin) / (dmax - dmin))

    return cs[data]



def Dividing_Line_Segments(arr, color):
    segments = np.stack((arr[:-1], arr[1:]), axis=1)
    cmap = color  # jet, hsv等也是常用的颜色映射方案
    colors = color_map(arr[1:, :-1], cmap)
    # colors = colors[::-1] #可以进行颜色的反转，但是有必要吗？
    line_segments = LineCollection(segments, colors=colors, linewidths=2, linestyles='solid', cmap=cmap)
    return line_segments


def cross_point(line1, line2):  # 计算交点函数
    x1 = line1[0]  # 取四点坐标
    y1 = line1[1]
    x2 = line1[2]
    y2 = line1[3]

    x3 = line2[0]
    y3 = line2[1]
    x4 = line2[2]
    y4 = line2[3]

    k1 = (y2 - y1) * 1.0 / (x2 - x1)  # 计算k1,由于点均为整数，需要进行浮点数转化
    b1 = y1 * 1.0 - x1 * k1 * 1.0  # 整型转浮点型是关键
    if (x4 - x3) == 0:  # L2直线斜率不存在操作
        k2 = None
        b2 = 0
    else:
        k2 = (y4 - y3) * 1.0 / (x4 - x3)  # 斜率存在操作
        b2 = y3 * 1.0 - x3 * k2 * 1.0
    if k2 == None:
        x = x3
    else:
        x = (b2 - b1) * 1.0 / (k1 - k2)
    y = k1 * x * 1.0 + b1 * 1.0
    return [x, y]

def drawplot(seg1,title, xlabel, ylabel):
    # sns.set(style='darkgrid')


    ax = plt.axes()
    ax.set_xlim((0,lx+double(lx/20.0)))
    ax.set_ylim((0,ly+double(ly/20.0)))
    if xlabel:
        plt.xlabel(xlabel)
    if ylabel:
        plt.ylabel(ylabel)
    plt.title(title,fontsize=15,fontweight = 'heavy')
    ax.add_collection(seg1)
    plt.show()

cut=14
color = "jet"
# cblabel = "asd"
# linewidth = 2

data = pd.read_excel(r'C:\Users\asus\Desktop\data.xlsx')
val=data.values
print(val)
plt.figure(1)
# plt.clf()
# plt.ylim(np.max(val[:,2]))
# plt.xlim(np.max(val[:,1]))
x=val[:,1]
y=val[:,2]
lx=np.max(x)
ly=np.max(y)
arr=np.stack([x,y])
print(arr)
arr=arr.T
print(arr.shape)
# seg=Dividing_Line_Segments(arr,color)
# print(seg)
#
# drawplot(seg,'P-I', 'Current (I/A)', 'Output power (P/mW)')




plt.figure(figsize=(20,20))
ax = plt.axes()
ax.set_xlim((-lx,lx))
ax.set_ylim((-ly,ly))
# plt.rcParams['font.sans-serif'] = ['SimHei']  # 绘图中文
plt.rcParams['axes.unicode_minus'] = False  # 绘图负号
plt.subplot().spines['left'].set_position('center')
plt.subplot().spines['bottom'].set_position('center')
plt.subplot().spines['right'].set_color('none')
plt.subplot().spines['top'].set_color('none')

plt.subplot().spines['left'].set_linewidth(4)
plt.subplot().spines['bottom'].set_linewidth(4)
plt.subplot().spines['right'].set_linewidth(4)
plt.subplot().spines['top'].set_linewidth(4)


plt.plot(x,y,'-k',linewidth=4)
z = np.polyfit(x[cut:], y[cut:], 1)
p = np.poly1d(z)
# 695.4 x - 12.27
plt.minorticks_on()
plt.subplot().grid(True, which='both')
plt.tick_params(pad=5,labelsize=20)     #调整坐标轴数字大小
print(p)
plt.plot([2, -2], p([2, -2]),'--',linewidth=4)
plt.title('P-I',weight='bold',fontsize=15)



z2 = np.polyfit(x[:11], y[:11], 1)
p2 = np.poly1d(z2)
print(p2)
lin1=[2,p(2),3,p(3)]
lin2=[2,p2(2),3,p2(3)]
c_pos=cross_point(lin1,lin2)
print(c_pos)
plt.plot([2, -2], p2([2, -2]),'--',linewidth=4)


# plt.annotate('a decision node',(0.02,5), (0.02,5), va="center",  ha="center",
#              xycoords="axes fraction", textcoords="axes fraction",
#              bbox=dict(boxstyle="sawtooth", fc="0.8"), arrowprops=dict(arrowstyle="<-"))

plt.annotate('cross_point', xy=(c_pos[0]+0.0005, c_pos[1]-0.2), xytext=(c_pos[0]-0.0005,-3.3), arrowprops=dict( arrowstyle="->",linewidth=4),fontsize=25,weight='bold')
plt.text(c_pos[0]-0.0007,-3.9,'(0.0177, 0.0533)',fontsize=25,weight='bold')

plt.plot(double(12.27/695.4),0,'o',linewidth=4)
plt.text(0.0086,6.5,'    695.4x-12.27',fontsize=30,weight='bold',fontdict={'style':'italic'})
# plt.text(double(12.27/695.4),0.2,'    (0.0176,0)',fontsize=20,weight='bold')
plt.text(-0.005,0.3,'2.852x+0.002785',fontsize=30,weight='bold',fontdict={'style':'italic'})

plt.show()






# plt.plot(val[:,1].tolist(),val[:,2].tolist())
# plt.title('P-I',fontsize=15,fontweight = 'heavy')
# plt.xlabel('Current (I/A)')
# plt.ylabel('Output power P/mW')
# plt.show()





# fig1 = plt.figure(1, figsize=(10, 8))
# plt.clf()
# plt.subplot(2, 1, 1)
# plt.plot(timems, res1[:, channel], 'r')
# plt.plot(timems, res2[:, channel], 'b:')
# plt.title('Time Course of P300 Amplitude')
# plt.legend(['standard', 'oddball'])
# plt.subplot(2, 1, 2)
# plt.plot(timems, ressq[:, channel])
# plt.title('Time Course of P300 r2')
# plt.xlabel('Time After Stimulus (ms)')
# plt.ylabel('Prop r2')
# plt.draw()


# print(pd.read_excel(xlsx, 'Sheet1'))