findcrest2.py－20170814

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Aug  2 14:28:09 2017

@author: vicky
"""

from numpy import *;
import numpy as np;
import pandas as pd;
import math 

#y=np.delete(y,range(152,163),axis=1)

[n,p]=np.shape(y); #输入y=原始数据，t=时间点

#------------------------平滑处理-------------------------
#f＝y的四个区间累计
f=np.zeros((n,p), dtype=np.int)
for i in range(0,n-3): 
    for j in range (0,p):
        f[i+3,j]=y[i,j]+y[i+1,j]+y[i+2,j]+y[i+3,j]
#-----------------------取f的二阶差分----------------------

#f的一,二，三阶差分df1,df2
df1=np.zeros((n,p), dtype=np.int)
df2=np.zeros((n,p), dtype=np.int)
#df3=np.zeros((n,p), dtype=np.int)
for j in range (0,p):
    for i in range (0,n-1):
        df1[i+1,j]=f[i+1,j]-f[i,j]
    for i in range (0,n-2):
        df2[i+2,j]=f[i+2,j]-2*f[i+1,j]+f[i,j]
#    for i in range(0,n-3):
#        df3[i+3,j]=f[i+3,j]-3*f[i+2,j]+3*f[i+1,j]-f[i,j]

#-------------------------波峰判断--------------------------
#起点波峰fmax=f的极大值且df2的局部极小值
del fmax
fmax=np.zeros((n,p), dtype=np.int)
for j in range(0,p):
    for i in range(2,n-2):
        if (f[i,j]>=max(f[0:i,j]) and 
            f[i,j]>=max(f[i+1:i+3,j]) and 
            (df2[i,j]<=min(df2[i-1,j],df2[i+1,j])  or 
             df2[i+1,j]<=min(df2[i,j],df2[i+2,j]) or 
             df2[i-1,j]<=min(df2[i-2,j],df2[i,j]))):
            fmax[i,j]=f[i,j]

#-------------------------时间段处理-------------------------
#ffmax=fmax去0，且8点以后,且>100;ttmax=ffmax对应的时间点
del ffmax
del ttmax
ffmax=np.zeros((n,p), dtype=np.int)
ttmax=np.zeros((n,p), dtype=np.int)
for j in range(0,p):  
    k=0
    for i in range(97,n):
        if (fmax[i,j]!=0 and fmax[i,j]>100):
            ffmax[k,j]=fmax[i,j]
            ttmax[k,j]=i
            k=k+1
            
#ffmax=fmax去0，且8点以后
#--------------------------去除假波峰-------------------------
[n1,p1]=np.shape(ffmax)
rffmax=np.zeros(ffmax.shape) #rffmax＝波峰的增长率
for j in range(0,p1):
    for i in range(1,n1):
        if ffmax[i,j]!=0:
            rffmax[i,j]=np.abs(ffmax[i,j]-ffmax[i-1,j])/ffmax[i-1,j]


#找阈值的下界c：取每列第1行到最大值之间的最大值b，再取所有列计算结果的最大值max(b)
del b
b=np.zeros((1,p))
for j in range(0,p1): 
    tm=np.array(np.where(rffmax[:,j]!=0))
    ma=np.array(np.where(rffmax[:,j]==max(rffmax[:,j])))[0,0]
    if size(tm)>2:
        b[0,j]=max(rffmax[0:ma,j])
    else:
        b[0,j]=0

c=np.max(b)

#去除rffmax大于阈值c的假波峰          
fffmax=np.zeros((n,p), dtype=np.int)
tmp=np.zeros((1,p))
for j in range(0,p):  
     tmp = np.array(np.where(rffmax[:,j]>c))
     if size(tmp)!=0:
         fffmax[0,j] = ffmax[tmp[0,0],j]
     else:
         fffmax[0,j] = 0
         
tttmax=np.zeros((1,p))         
for j in range(0,p): #波峰fffmax在f中的位置tttmax（对应时间）
    np.where(f[:,j]==fffmax[j])
    tttmax[j]=min(np.where(f[:,j]==fffmax[j]))

#-------------------------返回对应到y------------------------
yymax=np.zeros((1,p))
for j in range(0,p): #yymax＝在y中取tttmax往前4个内的最大值
    if tttmax[j]>4:
        yymax[j]=max(y(tttmax[j]-4:tttmax[j],j))

tt=np.zeros((1,p))
for j in range(0,p): #tt=yymax对应的时间点
    if yymax[j]!=0:
        np.where(y[:,j]==yymax[j]
    tt[j]=t[min(np.where(y[:,j]==yymax[j])),1]