Analizuojam signalus su python

Matlab DIE !!!

Taikyti python kasdieniams skaičiavimas pradėjau visai neseniai. Ankščiau mane labai tenkino ir octave ar calc ar maxima. Tačiau dabar, kai duomenų srauto atvaizdavimo poreikis kyla vis labiau, tenka ieškoti kitų galimų įrankių.

Pirma man atėjusi mintis ( kaip visad pirmosios mintys turi būti neteisingos ) buvo tiesiog pradėti rašyti savo programinį paketą. Tai užimtų amžius, bet galiausiai, pabaigoje turėčiau naudojamą įrankį.

Laiko stoka ir į pagalba tiesiog atėjo mano mylimas python! O su juo analizuoti signalus yra tikrai lengva panaudojus tiesiog vieną cmath modulį.

Ilgai nelaukiant, paimkime paprastą pavyzdį – neperiodinį, nesimetrinį signalą su faze. Tarkim mūsų signalas sklinda nuo 0 iki tau. Tuomet spektrinės analizės formulė atrodys

V(f) = A/(-i*2*pi*f)*(e**(-i*2*pi*f) -1)

Tikiuosi kas nors atpažįstate spektrinės analizės formą.

Kadangi dabar turim galutinę formulės formą, galima pradėti skaičiavimus! Paimkime amplitudę Atau = 2 V.

def result( f ):
  from cmath import pi, e, sqrt, polar
  A = 2. # 2V
  i = sqrt(-1)
  try:
    r, phi = polar ( A/(-i*2.*pi*f)*(e**(-i*2*pi*f) - 1) )
    return r, phi
  except ZeroDivisionError:
    return A, 0

Klaidos tikrinimą reikia daryti neatsitiktinai, o tikslingai, kadangi kai dažnis pasiekia nulį, mes turim dalinimą iš nulio, o tai sąlygoją klaidą. Kaip žinoma, funkcijos python gali gražinti keletą kintamųjų iškarto ir šiuo atvėju mums tai labai pagelbėja ( nors dirbančioje programoje turėjau išskaldyti r ir phi, dėl klaidos ).

Dabar teliko parašyti pagrindinę funkciją, kuri kreipiasi į mūsų parašytą ankstesnę funkciją.

def main():
  f = open("phase_shift_amplitude.dat", "w") # duomenu failas
  p = open("phase_shift_phase.dat", "w") # duomenu failas
  for fr in [x*0.001 for x in range(-3000,3001)]:
    r, phi = result(fr)
    print >>f, fr, r
    print >>p, fr, p

Galutinis, veikiantis variantas atrodo taip

#!/usr/bin/env python
def result ( f ):
  from cmath import pi, sqrt, polar
  A = 2. # 2V
  i = sqrt(-1)
  try:
    r, phi = polar ( A/(-i*2.*pi*f)*(e**(-i*2.*pi*f) - 1))
    return r
  except ZeroDivisionError:
    return A

def phiesult ( f ):
  from cmath import pi, e, sqrt, polar
  A = 2. # 2V
  i = sqrt(-1)
  try:
    r, phi = polar ( A/(-i*2.*pi*f)*(e**(-i*2.*pi*f) - 1))
    return phi
  except ZeroDivisionError:
    return 0

def main():
  import math
  f = open("phase_shift_amplitude.dat", "w") # duomenu failas
  p = open("phase_shift_phase.dat", "w") # duomenu failas
  for fr in [x*0.001 for x in range(-3000,3001)]:
    r = result(fr)
    phi = phiesult(fr)
    print >>f, fr, r
    print >>p, fr, phi

if __name__ == '__main__':
  main()

Dabar mes turime du failus su amplitudės ir fazės ( radianais ) skaitinėm vertėm. Teliko tik juos atvaizduoti. Tam galime pasitelkti gnuplot:

plot 'phase_shift_amplitude.dat' smooth frequency

plot 'phase_shift_phase.dat' smooth frequency

Arba LaTeX su pgfplots:

\begin{tikzpicture}
  \begin{axis} [
    width=500pt,
    height=230pt,
    grid=major,
    ylabel={$|V(f)|$},
    xlabel={$f, Hz$}
    ]
    \addplot[black] file {phase_shift_amplitude.dat};
  \end{axis}
\end{tikzpicture}\\
\textsl{Figure 1. Signal, with phase shift, amplitude spectrum.}\\

\begin{tikzpicture}
  \begin{axis} [
    width=500pt,
    height=230pt,
    grid=major,
    ylabel={$argV(f)$},
    xlabel={$f, Hz$}
    ]
    \addplot[black] file {phase_shift_phase.dat};
  \end{axis}
\end{tikzpicture}\\
\textsl{Figure 2. Signal phase, with phase shift.}\\

Asmeniškai aš daugiau esu linkęs prie Latex su pgfplots. Gnuplot labai geras įrankis greitam duomenų pateikimui. O jau atsakaitai geriausias įrankis yra LaTeX ( tikrų vyrų dokumentų tvarkyklė! ).

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s