Frekvenční analýza

vzorek signál (sinusovka s periodou 0.8 $\dot\ 2\pi \approx$ 5.) + šum (rovnoměrně rozdělená data)

In [10]:
%pylab inline 
from numpy.random import uniform
x1=uniform(size=1024)
from numpy import arange,sin
b=arange(1024)
y1=sin(b/0.8)

Populating the interactive namespace from numpy and matplotlib

C:\Users\Admin\Anaconda2\lib\site-packages\IPython\core\magics\pylab.py:161: UserWarning: pylab import has clobbered these variables: ['pi']
`%matplotlib` prevents importing * from pylab and numpy
  "\n`%matplotlib` prevents importing * from pylab and numpy"
In [11]:
plot(b,y1+x1,'.')
xlim(0,1024)
Out[11]:
(0, 1024)

fourierova transformace vzorku a čistého šumu

In [12]:
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)
Out[12]:
<matplotlib.text.Text at 0xcf65cc0>
In [13]:
from math import pi
pteor=0.8*2*pi
print("očekávaná frekvence %.5f Hz"%(1/pteor))

očekávaná frekvence 0.19894 Hz

náhodná procházka (odbočka :)

tvorba 1/f šumu

In [14]:
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()
Out[14]:
5.9735154549550247
In [15]:
# prevzorkovani sumu s 1/f sumem
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()
Out[15]:
(7.7619438416653947, 2.2877718829224349)
In [16]:
p=cumsum(uniform(-1,1,size=1024))
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"])
polyfit(np.log(r_[1:512]),np.log(ffp[1:512]),1),polyfit(np.log(r_[1:512]),np.log(ffp2[1:512]),1)
Out[16]:
(array([-0.94309579,  7.5644454 ]), array([-0.95086342,  9.80787595]))

jak narůstá rozptyl u náhodných procházek s rostoucí délkou

In [17]:
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)
Out[17]:
(array([  7.90501853,  13.27062508,  17.98484879,  22.54067068]),
 array([ 2.2383581 ,  3.51301196,  4.17368994,  6.24110727]))
In [18]:
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)
Out[18]:
(array([  6.43137671,  12.31831753,  15.59801301,  17.88216474]),
 array([ 1.87163557,  5.29105632,  5.79917662,  6.79705782]))

Nerovnoměrné vzorkování

použití Lomb-Scargle algoritmu (periodogram) čili LSSA (least square spectral analysis)
původní článek v Astrophys.Journal a pěkný ilustrativní rozbor 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.

In [19]:
from scipy.signal import lombscargle
bsel=uniform(size=b.size)<0.7# náhodný výběr bodů
In [20]:
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")
Out[20]:
<matplotlib.text.Text at 0xdc76fd0>

rekonstrukce periody z maximální frekvence

In [21]:
1/freqs[pgram.argmax()]
Out[21]:
0.80000000000000338

vlastní konstrukce periodogramu

  • pro každý bod dopočítáme fázi
  • průběh fázové křivky zrekonstruujeme jako průměr přes interval [t-win,t+win]
  • minimalizovaným kritériem je suma čtverců odchylek jednotlivých měření od rekonstruované křivky
In [22]:
def get_disp(peri,pos,data,ndiv=50,win=0.05,ret=1):
    phas=(pos/peri).astype(int)
    phas=pos/peri-phas
    #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))
    if ret==3: return array(rep),array(cnt),phas
    if ret==2: return array(rep),array(cnt)
    from scipy import interpolate
    prof=interpolate.interp1d(phax,array(rep)/array(cnt))
    chi2=sum((prof(phas)-data)**2)
    return chi2

ukázka "správné" periody a odchýlené o 0.01 (1.3%)

In [23]:
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")
xlabel("faze")
Out[23]:
<matplotlib.text.Text at 0xdc69908>
In [24]:
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")
xlabel("faze")
Out[24]:
<matplotlib.text.Text at 0xddbfc88>

aplikace na původní data (y1 a šum x1)

In [27]:
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()

rekonstrukce dlouhoperiodické křivky

In [28]:
pteor2=24.8*2*pi
y2=sin(b/pteor2*2*pi)
y3=sin(b/pteor2*2*pi/2.)*0.1
data2=y2[bsel]+y3[bsel]+x1[bsel]
plot(b[bsel],data2,'+')
Out[28]:
[<matplotlib.lines.Line2D at 0xd6e1630>]
In [29]:
pgram=get_disp(pteor2*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")
Out[29]:
<matplotlib.text.Text at 0xcf5bf60>
In [30]:
peris=r_[pteor2-20:pteor2+20:101j]
gmax=[get_disp(p,b[bsel],data2) for p in peris]
p0=argmin(gmax)
plot(peris/2/pi,gmax)
print "nalezené minimum"
peris[argmin(gmax)]/2/pi

nalezené minimum
Out[30]:
24.800000000000001

odhad nejistoty z profilu minima
(asymetrickou křivku prokládáme polynomem 4. stupně, z něj bereme 2. derivaci)

In [39]:
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),minprof2[0]
Out[39]:
(array([ 198.93206732,  186.9831774 ,  176.47811771,  167.41688823,
         159.79948896,  153.62591992,  148.89618109,  145.61027248,
         143.76819409,  143.36994592]), 79.49938434057141)
In [40]:
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.'])
Out[40]:
<matplotlib.legend.Legend at 0x105505f8>
In [41]:
# 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))
print("poloha minima a druhe derivace")
aroot.real[aroot.imag==0][0],mder,mder2,min2prof2[0]

poloha minima a druhe derivace
Out[41]:
(24.785699963082493, 162.4880253119656, 162.4880253119656, 79.49938434057141)

srovnání s FFT

zpět pro měření s rovnoměrným vzorkováním

In [45]:
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")
pteor2, 1/(b[5:51][argmax(abs(ft3[5:51]))]/1024.)
xlabel("frekvence")

ocekavana hodnota a perioda z polohy maxima
Out[45]:
<matplotlib.text.Text at 0xd542358>