concert.py－20170816

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 14 11:01:48 2017

@author: vicky
"""

from numpy import *
import numpy as np
import pandas as pd
import math 
import matplotlib.pyplot as plt

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

#波峰识别算法
#------------------------平滑处理-------------------------
#f＝y的四个区间累计
f=np.zeros((n,p), dtype=np.int)
for i in range(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 (p):
    for i in range (n-1):
        df1[i+1,j]=f[i+1,j]-f[i,j]
    for i in range (n-2):
        df2[i+2,j]=f[i+2,j]-2*f[i+1,j]+f[i,j]

#-------------------------波峰判断--------------------------
#起点波峰fmax=f的极大值且df2的局部极小值
del fmax
fmax=np.zeros((n,p), dtype=np.int)
for j in range(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=ffmax[~np.all(ffmax == 0, axis=1)]
ttmax=ttmax[~np.all(ttmax == 0, axis=1)]

#--------------------------去除假波峰-------------------------
[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]

del tmp
rffmax[1,:]=0 #从第2个波峰起开始判断
#fffmax=去除rffmax大于2.223的假波峰，起点波峰之后的波峰,tttmax=对应时间点
fffmax=np.zeros((n,p), dtype=np.int)
tttmax=np.zeros((n,p), dtype=np.int)
for j in range(0,p):
    tmp=np.array(np.where(rffmax[:,j]>2.223))
    k=0
    if size(tmp)>0:
        for i in range(tmp[0,0],n1):
            fffmax[k,j]=ffmax[i,j]
            tttmax[k,j]=ttmax[i,j]
            k=k+1
    else:
        fffmax[:,j]=0
        tttmax[:,j]=0

fffmax=fffmax[~np.all(fffmax == 0, axis=1)]
tttmax=tttmax[~np.all(tttmax == 0, axis=1)]        
[n2,p2]=np.shape(fffmax)

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

tymax=zeros([1,p],dtype=t.dtype)
for j in range(0,p):
    tymax[0,j]=t[tt[0,j]]
#--------------------削去中间峰值，并加上随机抖动----------------
ff=f; #ff＝原始数据f去除峰值
for j in range(p):
    if tttmax[0,j]>0:
        for i in range(tttmax[0,j],n):
            if f[i,j]>fffmax[0,j]:
                ff[i,j]=fffmax[0,j]*(0.99+0.02*np.random.rand(1,1))

#起点峰值下降算法

#--------------------------波谷判断--------------------------
fmin=zeros(shape(f),dtype=int) #所有波谷fmin=f的局部极小值且df2的局部极大值
for j in range(p):
    if tttmax[0,j]>0:
        for i in range(n-2):
            if f[i,j]<=min(f[i-1,j],f[i+1,j]):
                fmin[i,j]=f[i,j]

ffmin=zeros(shape(fffmax),dtype=int)
ttmin=zeros(shape(tttmax),dtype=int)
for j in range(p):  #区间波谷ffmin=每个波峰之间只取一个波谷后的fmin，去0;ttmin=波谷时间点
    for k in range(n2-1):
        if tttmax[k,j]>0 & tttmax[k+1,j]>0:
            for i in range(tttmax(k,j),tttmax(k+1,j)):
                if fmin[i,j]!=0:
                    ffmin[k,j]=fmin[i,j]
                    ttmin[k,j]=i
                    break
        else:
            tttmax[k,j]>0 & tttmax[k+1,j]==0
            for i in range(tttmax[k,j],n):
                if fmin[i,j]!=0: #& fmin(i,j)>100
                    ffmin[k,j]=fmin[i,j]
                    ttmin[k,j]=i
                    break

#-------------------------降低异常值-------------------------
fd=f #fd=把第一个异常波调整成平均值
fd=fd.astype(np.float64)
for j in range(0,p):
    if tttmax[0,j]>0 and ttmin[0,j]>0:
        for i in range(2*tttmax[0,j]-ttmin[0,j],ttmin[0,j]):
            fd[i,j]=(f[i,j]+ffmin[0,j])/2

s1=0;k1=0; #r1=法1峰值下降率
for j in range(p):
    if fffmax[0,j]>0 and tttmax[0,j]>0 and ttmin[0,j]>0:
        s1=s1+(fffmax[0,j]-max(fd[2*tttmax[0,j]-ttmin[0,j]:ttmin[0,j],j]))/fffmax[0,j]
        k1=k1+1
r1=s1/k1

#---法2:波峰波谷中间时间点为区间端点
tm1=zeros(shape(tttmax),dtype=int)#tm1=区间起点时间，tm2＝区间终点时间
tm2=zeros(shape(tttmax),dtype=int)
for j in range(p):
    for i in range(n2):
        if tttmax[i,j]>0 and ttmin[i,j]>0:
            tm1[i,j]=floor((3*tttmax[i,j]-ttmin[i,j])/2)
            tm2[i,j]=floor((tttmax[i,j]+ttmin[i,j])/2)

fm=zeros(shape(tttmax))#fm=波峰前后区间中间端点的均值
for j in range(p):
    for i in range(n2):
        if tm1[i,j]>0 and tm2[i,j]>0:
            fm[i,j]=(f[tm1[i,j],j]+f[tm2[i,j],j])/2

fe=f; #fe=f和fm取均值
fe=fe.astype(np.float64)
for j in range(p):
    if tm1[0,j]>0 and tm2[0,j]>0:
        for i in range(tm1[0,j],tm2[0,j]):
            fe[i,j]=(f[i,j]+fm[0,j])/2

s2=0;k2=0; #峰值下降率
for j in range(p):
    if fffmax[0,j]>0 and tm1[0,j]>0 and tm2[0,j]>0:
       s2=s2+(fffmax[0,j]-max(fe[tm1[0,j]:tm2[0,j],j]))/fffmax[0,j]
       k2=k2+1
r2=s2/k2

#----------------------削去第二个峰值，且加上随机波动------------------
fdmax=zeros((1,p))#fdmax=新的起点波峰
for j in range(p):
    if ffmax[0,j]>0 and ffmin[0,j]>0:
        fdmax[0,j]=(fffmax[0,j]+ffmin[0,j])/2

ffd=fd;#ffd＝fd中大于fdmax的降低为fdmax加上一个误差小于1%的随机波动
for j in range(p):
    if tttmax[0,j]>0 and fdmax[0,j]>0:
        for i in range(tttmax[0,j],n):
            if ffd[i,j]>fdmax[0,j]:
                ffd[i,j]=fdmax[0,j]*(0.99+0.02*np.random.rand(1,1))

femax=zeros((1,p)) #femax=新的起点波峰
for j in range(p):
    if ffmax[0,j]>0 and fm[0,j]>0:
        femax[0,j]=(fffmax[0,j]+fm[0,j])/2

ffe=fe;#ffe＝fde中大于femax的降低为femax加上一个误差小于1%的随机波动
for j in range(p): 
    if tttmax[0,j]>0 and femax[0,j]>0:
        for i in range(tttmax[0,j],n):
            if ffe[i,j]>femax[0,j]:
                ffe[i,j]=femax[0,j]*(0.99+0.02*np.random.rand(1,1))
                
#---------------------------图例--------------------------------
td=range(0,288)
for j in range(p):
    fig = plt.plot(td,ffd[:,j],td,ffe[:,j],td,f[:,j])
    plt.savefig('/Users/vicky/Documents/地图/演唱会/pythonfig/'+np.str(j)+'.jpg')
    plt.show()