vzorek signál (sinusovka s periodou 0.8 $\dot\ 2\pi \approx$ 5.) + šum (rovnoměrně rozdělená data)
%pylab inline
from numpy.random import uniform
import matplotlib
matplotlib.rcParams['figure.figsize'] = [10, 5]
x1=uniform(size=1024)
from numpy import arange,sin
b=arange(1024)
y1=sin(b/0.8)
plot(b,y1+x1,'.')
xlim(0,1024)
fourierova transformace vzorku a čistého šumu
from matplotlib.pyplot import psd,plot
from numpy import fft
f1=fft.fft(x1)
ft1=fft.fft(x1+y1)
gr1=plot(b[:512]/1024.,abs(f1[:512]),b[:512]/1024.,abs(ft1[:512]))
ax=gr1[0].axes
ax.set_xlim(0.1,0.4)
ax.grid(1)
legend([u"šum",u"signál + šum"])
xlabel("frekvence")
#a1=psd(x1)
from math import pi
pteor=0.8*2*pi
print("očekávaná frekvence %.5f Hz"%(1/pteor))
tvorba 1/f šumu
#b=arange(1024)
f1corr=f1.copy()
f1corr[1:]/=b[1:]/1024.
f2corr=f1.copy()
f2corr[1:]/=b[1:]/330.
#f1corr[0]=0
x1corr=fft.ifft(f1corr).real
x2corr=fft.ifft(f2corr).real
plot(b,x1corr,b,x2corr)
legend(['1024','330'])
x1corr.std()
# prevzorkovani sumu s 1/f spektrem
def refreq(samp,npt=1024):
if npt==0: npt=samp.size
f1=fft.fft(samp,npt)
f1[1:]/=arange(1,npt)/float(npt)
return fft.ifft(f1).real
t1=array([refreq(uniform(size=1024)).std() for i in range(32)])
t1.mean(),t1.std()
p=cumsum(uniform(-1,1,size=1024)) #nahodna prochazka
ffp=abs(fft.fft(p))
p2=cumsum(normal(0,1,size=2*1024))
ffp2=abs(fft.fft(p2))
loglog(ffp[:512])
loglog(ffp2[:512])
title("Fourierovsky obraz kumulativnich souctu (stredni hodnota 0)")
legend(["1024 kroku","2048 kroku"])
#prvni parametr odpovida linear. clenu v log_log skale (1/x zavislost)
polyfit(np.log(r_[1:512]),np.log(ffp[1:512]),1),polyfit(np.log(r_[1:512]),np.log(ffp2[1:512]),1)
jak narůstá rozptyl u náhodných procházek s rostoucí délkou
t2=array([[refreq(uniform(size=1024*s),npt=0).std() for i in range(32)] for s in range(1,8,2)])
t2.mean(1),t2.std(1)
t3=array([[cumsum(uniform(-1,1,size=1024*s)).std() for i in range(32)] for s in range(1,8,2)])
t3.mean(1),t3.std(1)
použití Lomb-Scargle algoritmu (periodogram) čili LSSA (least square spectral analysis)
původní článek v Astrophys.Journal a pěkný ilustrativní rozbor Jakea Vanderplase
na ArXivu
rychlá implementace v pythonu zde
klasický periodogram vychází z výpočtu
$$P=\left|\frac{1}{N}\sum_{i=1}^N{g_i e^{-i\ 2 \pi \omega t_i}} \right|^2 =\\ = \frac{1}{N} \left[ \left( \sum_{i=1}^N{g_i \cos ( 2 \pi \omega t_i)}\right)^2 + \left(\sum_{i=1}^N{g_i \sin (2 \pi \omega t_i)} \right)^2\right ]$$
u nerovnoměrného vzorkování není zajištěna ortogonalita členů s $\cos$ a $\sin$; Scargle navrhl zobecnění pomocí fázového posunu $\tau$ a normalizačních parametrů $A, B$ do tvaru $$P= \frac{A^2}{2} \left( \sum_{i=1}^N{g_i \cos (2 \pi \omega (t_i-\tau))}\right)^2 + \frac{B^2}{2} \left(\sum_{i=1}^N{g_i \sin (2 \pi \omega (t_i-\tau))} \right)^2$$
Při volbě $A^2=1/\sum_{i=1}^N{ \cos^2 ( 2 \pi \omega (t_i-\tau))}$, $B^2=1/\sum_{i=1}^N{ \sin^2 ( 2 \pi \omega (t_i-\tau))}$ a $$\tau = \frac{1}{4 \pi \omega} \arctan\left(\frac{\sum_{i=1}^N{ \sin ( 4 \pi \omega t_i)}}{\sum_{i=1}^N{ \cos ( 4 \pi \omega t_i)}} \right)$$ dostaneme periodogram, který není závislý na volbě počátku časové osy (vliv konkrétních vzorkovacích bodů $t_i$ odstraňuje právě volba $\tau$) a pro rovnoměrné vzorkování přechází do tvaru klasického periodogramu.
from scipy.signal import lombscargle
bsel=uniform(size=b.size)<0.7# náhodný výběr bodů
freqs=arange(1,1.5,0.005)
data=y1[bsel]+x1[bsel]
data-=data.mean()
data/=data.std()
pgram=lombscargle(b[bsel].astype('float'),data,freqs)
gr2=semilogy(freqs,pgram)[0]
gr2.axes.set_ylim(-20,350)
xlabel("frequence")
freqs=arange(1.2,1.3,0.001)
data=y1[bsel]+x1[bsel]
data-=data.mean()
data/=data.std()
pgram2=lombscargle(b[bsel].astype('float'),data,freqs)
gr2=semilogy(freqs,pgram2)[0]
gr2.axes.set_ylim(-20,350)
xlabel("frequence")
gr2=plot(freqs,pgram2)[0]
title("linearni skala")
xlabel("frequence")
rekonstrukce periody z maximální frekvence
1/freqs[pgram.argmax()]
[t-win,t+win]
def get_disp(peri,pos,data,ndiv=50,win=0.05,ret=1):
#defaultne vrati sumu ctvercu rezidui bodu a stredni hodnoty
nperi=(pos/peri).astype(int)
phas=pos/peri-nperi
#plot(phas,data,'.')
phax=r_[:1:(ndiv+1)*1j]
rep=[]
cnt=[]
for a in phax:
sel=abs(phas-a)<win
if a<win: sel[phas>1-win+a]=True
if a>1-win: sel[phas<win-1+a]=True
cnt.append(sum(sel))
rep.append(sum(data[sel]))
#plot(phax,array(rep)/array(cnt))
cnt=array(cnt)
cnt[cnt==0]=1
if ret==4: return array(rep)/cnt
if ret==3: return array(rep),cnt,phas
if ret==2: return array(rep),cnt
from scipy import interpolate
prof=interpolate.interp1d(phax,array(rep)/cnt)
chi2=sum((prof(phas)-data)**2)
return chi2
pgram=get_disp(pteor,b[bsel],data,ret=3)
ndiv=50
plot(pgram[2],data,'.')
plot(r_[:1:(ndiv+1)*1j],pgram[0]/pgram[1],lw=2)
title("periodogram a stredni funkce p_t=%.3f"%pteor)
xlabel("faze")
pgram=get_disp(pteor-0.01,b[bsel],data,ret=3)
ndiv=50
plot(pgram[2],data,'.')
plot(r_[:1:(ndiv+1)*1j],pgram[0]/pgram[1],lw=2)
title("periodogram a stredni funkce p_t=%.3f"%(pteor-0.01))
xlabel("faze")
aplikace na původní data (y1
a šum x1
)
peris=r_[pteor-0.8:pteor+0.8:41j]
gmax=[get_disp(p,b[bsel],data) for p in peris]
plot(peris/2/pi,gmax)
xlabel("perioda")
ylabel("rezidualni rozptyl")
title("hledani char. periody")
grid()
peris=r_[pteor-0.09:pteor+0.09:41j]
gmax=[get_disp(p,b[bsel],data) for p in peris]
plot(peris/2/pi,gmax)
xlabel("perioda")
ylabel("rezidualni rozptyl")
title("hledani char. periody")
grid()
zp2p=lambda a:a.max()-a.min()
#zp2p(get_disp(p,b[bsel],data))
gamp=[zp2p(get_disp(p,b[bsel],data,ret=4)) for p in peris]
plot(peris/2/pi,gamp)
xlabel("perioda")
grid()
title("amplituda prumerovaneho fazogramu (peak2peak)")
pteor2=24.8*2*pi
y2=sin(b/pteor2*2*pi)
y3=sin(b/pteor2*3*pi/2.)*0.1
data2=y2[bsel]+y3[bsel]+x1[bsel]
plot(b[bsel],data2,'+')
pgram=get_disp(pteor2,b[bsel],data2,ret=3)
ndiv=50
plot(pgram[2],data2,'.')
plot(r_[:1:(ndiv+1)*1j],pgram[0]/pgram[1],lw=2)
title("periodogram a stredni funkce")
xlabel("faze")
pgram=get_disp(pteor2+0.4*2*pi,b[bsel],data2,ret=3)
ndiv=50
plot(pgram[2],data2,'.')
plot(r_[:1:(ndiv+1)*1j],pgram[0]/pgram[1],lw=2)
title("periodogram a stredni funkce")
xlabel("faze")
peris=r_[pteor2-40:pteor2+40:41j]
gmax=[get_disp(p,b[bsel],data2,ret=1) for p in peris]
plot(peris/2/pi,gmax)
xlabel("perioda")
ylabel("rezidualni rozptyl")
title("hledani char. periody")
grid()
odhad nejistoty z profilu minima
(asymetrickou křivku prokládáme polynomem 4. stupně, z něj bereme 2. derivaci)
p0=np.argmin(gmax)
psel=peris[p0-5:p0+5]/2/pi
minprof=polyfit(psel,gmax[p0-5:p0+5],4)
minprof2=polyfit(psel,gmax[p0-5:p0+5],2)
mder=polyval(polyder(polyder(minprof)),psel).mean() #profil 2. derivace
#polyval(polyder(polyder(minprof)),psel)
aroot=roots(polyder(minprof)) # koreny 1. derivace
print("poloha minima a druhe derivace")
aroot.real[aroot.imag==0][0],mder,minprof2[0]
plot(psel,gmax[p0-5:p0+5],'+',ms=15)
plot(psel,polyval(minprof,psel))
plot(psel,polyval(minprof2,psel))
xlabel("perioda")
legend(['data','pol.4.stup.','pol.2.stup.'])
# srovnávací test s prostou sinusovkou (bez y3)
gmax=[get_disp(p,b[bsel],x1[bsel]+y2[bsel]) for p in peris]
min2prof=polyfit(psel,gmax[p0-5:p0+5],4)
min2prof2=polyfit(psel,gmax[p0-5:p0+5],2)
mder2=polyval(polyder(polyder(min2prof)),psel).mean()
aroot=roots(polyder(min2prof))
aroot.real[aroot.imag==0][0],mder,mder2,min2prof2[0]
def phas_test(yval):
bsel=uniform(size=b.size)<0.7
gmax=[get_disp(p,b[bsel],uniform(size=sum(bsel))+yval[bsel]) for p in peris]
aroot=roots(polyder(polyfit(psel,gmax[p0-5:p0+5],4)))
return aroot.real[aroot.imag==0][0]
btot=np.array([phas_test(y2+y3) for i in range(200)]).T
2*pi*btot.mean()
2*pi*btot.std()
zpět pro měření s rovnoměrným vzorkováním
subias=(x1+y2).mean()
ft3=fft.fft(x1+y2-subias)
gr3=plot(b[:51]/1024.,abs(ft3[:51]))
print("ocekavana hodnota a perioda z polohy maxima")
xlabel("frekvence")
pteor2, 1/(b[5:51][argmax(abs(ft3[5:51]))]/1024.)
def four_test(yval):
ft3=fft.fft(uniform(size=1024)+yval)
mpos=3+argmax(abs(ft3[3:51]))
if mpos<3:mpos=3
ids=np.polyfit(r_[mpos-3:mpos+4],abs(ft3[:51])[mpos-3:mpos+4],2)
return mpos,-ids[1]/2/ids[0]/1024.
ztot,atot=np.array([four_test(y2+y3) for i in range(200)]).T
1/atot.mean(),atot.std()/atot.mean()**2
1/atot.mean()/2/pi