note[1]：python Hilbert transform and calculate phase

系统输出两个相位不同的信号，要求用hilbert transform求相位差，参考

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html

"""
Created on Tue Feb 13 23:07:25 2018

@author: JZY
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert

duration = 1.0
fs = 1000.0
samples = int(fs*duration)
t = np.arange(samples) / fs
signal = np.cos(2.0*np.pi*50*t+(np.pi/180*75))

signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) )
analytic_signal = hilbert(signal)

amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) /
                           (2.0*np.pi)*fs )
ph = instantaneous_phase*180/np.pi
#ph = (instantaneous_phase-2.0*np.pi*50*t)*180/np.pi
fig = plt.figure()
ax0 = fig.add_subplot(311)
ax0.plot(t, signal, label='signal')
ax0.plot(t, amplitude_envelope, label='envelope')
ax0.set_xlabel("time in seconds")
ax0.legend()
ax1 = fig.add_subplot(312)
ax1.plot(t[1:], instantaneous_frequency)
ax1.set_xlabel("time in seconds")
#ax1.set_ylim(0.0, 120.0)
ax1 = fig.add_subplot(313)
print(ph[0])
ax1.plot(t, ph)
ax1.set_xlabel("time in seconds")

由于信号的频率是已知的，我只需要测信号的增益和相位就行，按照变换后分析信号的瞬时相位为2*pi*f*t+w，code中ph的第一个是t=0是的初相w。另外在初相那里开始时自己犯了一些常识错误！另外hilbert transform用来求取信号相位有一定的局限性，它需要的信号要有较高的SNR才能保证精度，在仿真时将输入信号经过一个Passband很窄的IIR bandpass Filter再求相位才稍稳定下来(目前还在调下位机，调python仿真，定量的精确数据我还无法给出)。当然用fft再求信号相位的方法还在测试，到时还要对比这两种方法用单片机来处理速度怎样。
在实际MCU对采样信号处理时这里并不会直接使用hilbert变换计算，因为变换需要FFT或者滤波器方式。实际上用到的方法反而类似IQ解调，就像下面给出的图片一样，这样MCU处理起来反而更方便：

上面图片是一路信号处理，两个信号做同样处理之后，再做简单的乘，和，差，除就就可以得到他们的相差和相对增益了，很方便！

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy
from scipy.signal import hilbert
from ctypes import *
fs = 48000
fo=5000
amp_ref=3000

t=np.arange(0,48*(1/fs),1/fs)

y=np.array(np.sin(2*np.pi*fo*t))
y1=np.array(np.cos(2*np.pi*fo*t))
r=np.array(10*np.cos(2*np.pi*fo*t+(np.pi/180*19)))
s=np.array(30*np.cos(2*np.pi*fo*t+(np.pi/180*90)))

hr=hilbert(r)
hs=hilbert(s)
o1=y1*s#cos(wt)*Samp
o2=y*s#sin(wt)*Samp

r1=y1*r#cos(wt)*Ref
r2=y*r#sin(wt)*Ref

rs=0
rc=0
ss=0
sc=0
for i in o1:#cos(wt)*Samp
    ss+=i#Acos(dif)
for i in o2:#sin(wt)*Samp
    sc+=i#Asin(dif)
for i in r1:#cos(wt)*Ref
    rs+=i#Acos(dif)
for i in r2:#sin(wt)*Ref
    rc+=i#Asin(dif)
rr=rc*rc+rs*rs
gmar=(sc * rc + ss * rs) / rr
gmai=(ss * rc - sc * rs) / rr
print('QD:%f' %(np.arctan(gmai/gmar)*180/np.pi) )
print('abs(gmai/gmar):%f' %np.sqrt(gmar**2+gmai**2))
#P=math.atan(np.angle(hr))
#pr = np.unwrap(P)
#= np.unwrap(np.angle(analytic_signal))
pr = np.unwrap(np.angle(hr)-2*np.pi*fo*t)*180/np.pi
ps = np.unwrap(np.angle(hs)-2*np.pi*fo*t)*180/np.pi
print(np.arctan(-rc/rs)*180/(math.pi))
print(np.arctan(-sc/ss)*180/(math.pi))
print('hilbert:ps-pr %f'%(ps[0]-pr[0]))

#plt.figure(1)
#plt.plot(gmar,color='r')
#plt.figure(2)
#plt.plot(t,y1)
#plt.figure(3)
#plt.plot(t,s)
#plt.figure(4)
#plt.plot(t,hr,t,hs)
plt.show()

一些注意点：

http://blog.csdn.net/wordwarwordwar/article/details/62078726
上面是关于fft分析初相时要注意采样点数，采样频率，输入信号频率的问题