Geneerinen proseduuri on yhteinen nimi joukolle eri proseduureja, joista valitaan suoritettavaksi yksi todellisten argumenttien tyyppien perusteella.
Eri proseduurien on erottava toisistaan argumenttien tyyppien suhteen, jotta todellisten argumenttien perusteella voidaan päättää, mikä proseduuri suoritetaan. Pelkkä argumenttien järjestyksen vaihtaminen ei riitä (koska voidaan kutsua avainsana-argumenttien avulla).
module div2
private
interface half
module procedure half1, half2
end interface
public :: half
contains
integer function half1(i)
integer, intent(in) :: i
half1 = i/2
end function
real function half2(x)
real, intent(in) :: x
half2 = x/2.0
end function
end module
Proseduurin nimi half viittaa kahteen
eri proseduuriin, joista suoritetaan jompikumpi.
Tässä half1 ja half2 ovat yksityisiä, joten
niitä ei voi kutsua suoraan.
program divtest
use div2
integer :: i=10
real :: x=1.0
i=half(i)
x=half(x)
write(*,*) i,x
end program
Esimerkiksi yhteenlasku: määritellään yhteenlaskun suorittava funktio erikseen kullekin mahdolliselle tyyppikombinaatiolle.
interface operator (+)
module procedure ratint_add, intrat_add, add
end interface
...
contains
...
function add (p, q)
! lasketaan rationaalilukujen p ja q summa
type (rational) :: add
type (rational), intent(in):: p, q
...
end function
function ratint_add (p, i)
! lasketaan kokonaisluvun i ja rationaaliluvun p summa
type (rational) :: ratint_add
type (rational), intent(in):: p
integer, intent(in) :: i
integer a,b,c,d
type (rational) :: r
a=p%nominator; b=p%denominator
c=i; d=1
r%nominator = a*d + b*c
r%denominator = b*d
ratint_add = simplify (r)
end function
function intrat_add (i, q)
! lasketaan kokonaisluvun i ja rationaaliluvun q summa
type (rational) :: intrat_add
integer, intent(in) :: i
type (rational), intent(in):: q
integer a,b,c,d
type (rational) :: r
a=i; b=1
c=q%nominator; d=q%denominator
r%nominator = a*d + b*c
r%denominator = b*d
intrat_add = simplify (r)
end function
Jos sijoitus edellyttää epästandardia tyyppimuunnosta, sijoitusoperaattorin toiminta on määriteltävä.
Esimerkiksi sijoitus "rationaaliluku = kokonaisluku" ei ole automaattisesti määritelty, vaan moduulissa on kerrottava, miten tyyppimuunnos tehdään.
interface assignment (=)
module procedure intrat_subst
end interface
...
subroutine intrat_subst (p, i)
! sijoitetaan kokonaisluvun i arvo rationaaliluvuksi p
integer, intent(in) :: i
type (rational), intent(out) :: p
p%nominator = i
p%denominator = 1
return
end subroutine
Vain pisteellinen notaatio (kuten .eq., .and. jne)
interface operator (.inv.)
module procedure inverse
end interface
...
function inverse (p)
! lasketaan rationaaliluvun p kaanteisluku
type (rational) inverse
type (rational), intent(in):: p
type (rational) r
r%denominator=p%nominator
r%nominator=p%denominator
inverse = simplify (r)
return
end function
Mieti, mitkä määrittelyt on syytä pitää yksityisinä. Käyttäjän ei tarvitse päästä suoraan käsiksi kaikkeen, koska seurauksena voi olla hankalasti jäljitettäviä virheitä. Erityisesti mutkikkaita tietorakenteita käsiteltäessä on parempi antaa käyttöön vain joukko testattuja ja luotettavia työkaluja.
Kannattaa harkita huolella, mitä omien muuttujatyyppien laskutoimituksia toteuttaa operaattoreina. Liiallinen operaattorien käyttö saattaa hämärtää käsitystä, mitä ohjelma oikeastaan tekee.
Laskutoimitusten tulisi olla luonnollisia (mitä puhelinluetteloiden kertolaskun pitäisi tehdä?)
real, pointer, dimension (:) :: p
p on osoitin, joka voi osoittaa yksiulotteiseen
taulukkoon. Taulukolle ei varata tilaa tässä vaiheessa.
Osoitin asetetaan osoittamaan muuttujaa operaattorilla =>.
Osoitin ei ole pelkkä muuttujan talletuspaikan osoite, vaan kuvaaja, joka sisältää muutakin tietoa. Osoitin voi viitata myös taulukkosektioon:
real, dimension(100,100), target:: x
real, pointer, dimension (:) :: p1
real, pointer, dimension (:,:) :: p2
p1 => x(1,:)
p1 => x(:,10)
p1 => x(1, 1:91:10)
p2 => x
p2 => x(2:3, 10:20)
p2 => x(1:11:2, 10:100:10)
Fortranissa osoitin on lähinnä synonyymi sen
osoittaman olion nimelle.
HUOM!
p=p+1
C:ssä osoitin siirtyy osoittamaan seuraava alkiota.
Fortranissa p:n osoittamaan olioon lisätään 1.
Osoittimen tila (association status) voi olla:
real, target :: x
real, pointer :: p
! p maarittelematon
p => x
! p on liitetty x:aan
nullify(p)
! p on tyhja
Jos p on liitetty johonkin muuttujaan,
associated(p) on tosi.
Osoittimien käyttö dynaamiseen tilanvaraukseen:
real, dimension(:), pointer :: p
...
allocate (p(1000))
...
deallocate(p)
Sopii hyvin linkitettyihin tietorakenteisiin.
Esimerkiksi binääripuu:
type node
character (len=100) :: name
type (node), pointer :: left, right
end type node
type (node), pointer :: root, l, r
allocate(root, l, r)
root%name='puun juuri'
root%left => l ! linkki vasempaan alipuuhun
root%right => r ! linkki oikeaan alipuuhun
Esimerkki: linkitetyn listan käsittely.
program linkedlist
! rakennetaan ja tulostetaan lista, joka sisaltaa
! syotetyt luvut kasvavassa suuruusjarjestyksessa
implicit none
type node
integer :: value
type (node), pointer :: next
end type node
integer :: num
type (node), pointer :: first, current, p, q
! first osoittaa listan alkuun
! luodaan aluksi tyhja lista
nullify(first)
! luetaan alkioita ja lisataan listan alkuun
! lopetetaan, kun annettu luku on 0
do
read (*,*) num
if (num==0) exit
allocate(current) ! luodaan uusi tietue
current%value=num ! asetetaan arvokentta
! liitetaan luotu tietue listaan oikealle paikalleen
if (.not.associated(first)) then
! lista on tyhja
current%next=>first
first=>current
else if (num <= first%value) then
! luettu luku on pienempi kuin listan ensimmainen luku
! uusi tietue tulee listan alkuun
current%next=>first
first=>current
else
! muuten ruvetaan etenemaan pitkin listaa
! uusi tietue p:n ja q:n osoittamien tietuiden valiin
q=>first; p=>first%next
do
if (.not.associated(p)) exit ! lista kayty loppuun
if (num <= p%value) exit ! loytyi oikea vali
q=>p; p=>p%next ! siirrytaan eteenpain
end do
! linkitetaan tietue p:n ja q:n valiin
current%next=>p
q%next=>current
end if
end do
! lista on valmis, tulostetaan
current=>first ! listan 1. tietue
! edetaan listaa, kunnes vastaan tulee tyhja linkki
do
if (.not.associated(current)) exit
write (*,*) current%value
current=>current%next ! seuraava tietue
end do
end program linkedlist
Osoittimien mahdollisia ongelmia:
Roikkumaan jäänyt viite (dangling reference)
real, pointer :: p, q
...
allocate (p)
p=1.0
q=>p
deallocate(p) ! q = ?
Saavuttamaton muistialue:
real, dimension(:), pointer :: p
...
allocate (p(1000))
...
nullify(p)
Taulukolle on varattu tila, mutta siihen ei voi
enää viitata millään keinolla. Tila vapautuu
vasta ohjelman päättyessä.
dimension: määrittelee taulukon koon
parameter: muuttuja on vakio, jota ei saa muuttaa
save: muuttuja varataan staattisesti ja sen arvo säilyy aliohjelmakutsujen välillä
allocatable: dynaamisesti varattava taulukko
pointer: osoitin
target: muuttujaan saatetaan viitata osoittimen avulla
intent: aliohjelman argumentin käyttötapa
optional: valinnainen argumentti; voi puuttua proseduurin kutsusta
external: muuttuja on proseduurin nimi
intrinsic: muuttuja on varusfunktion nimi
private: moduulin yksityinen muuttuja
public: moduulin julkinen muuttuja
program leivo
implicit none
integer :: jauhot=1, sokeri=1, munat=1
namelist /kakku/ jauhot, sokeri, munat
...
read(*, nml=kakku)
...
write(*, nml=kakku)
end program
Syöttötiedoston sisältö voisi olla
&kakku munat=3, jauhot=4/
Tulostus:
&KAKKU
JAUHOT=4, SOKERI=1, MUNAT=3
/
Sopii esimerkiksi tilanteeseen, jossa useimmilla muuttujilla
käyttökelpoiset alkuarvot. Syöttötietoina tarvitsee antaa
vain muutettavat arvot.
F(u)= \int-\infty\infty f(x) e-2\pi iux dx.
ja käänteismuunnos
f(x)=\int-\infty\infty F(u) e2\pi iux du.
Muitakin määritelmiä käytetään yleisesti. Integraalin edessä voi olla vakiokertoimia, kuten 1/2\pi tai 1/\sqrt(2\pi). Eksponentin tekijä 2\pi voi puuttua tai merkit voivat olla toisinpäin. Tarkista aina, mitä määritelmää milloinkin on käytetty!
(h * f)(x) = \int_-\infty\infty h(x-x') f(x') dx'.
Konvoluution Fourier'n muunnos on funktioiden muunnosten tulo:
F (h * f)(u) = (F h)(u) (F f)(u) = H(u)F(u).
Funktioiden f ja g korrelaatio Rfg on
Rfg = \int-\infty\infty f(x') g(x + x') dx'.
Funktion korrelaatio itsensä kanssa on autokorrelaatiofunktio:
Rff = \int-\infty\infty f(x') f(x+x') dx',
jonka Fourier'n muunnos on
F Rff = | F(u)|2.
Jos signaali sisältää taajuuksia, jotka ovat suurempia kuin kaksi kertaa näytteidenottotaajuus, signaalia ei voi enää rekonstruoida yksikäsitteisellä tavalla. Nämä korkeat taajuudet saattavat aiheuttaa alias-ilmiön (aliasing), eli rekonstruoituun signaaliin ilmaantuu virheellisen tulkinnan vuoksi taajuuksia, joita alkuperäisessä signaalissa ei ole.
Diskreetti Fourier'n muunnos määritellään kaavalla
Fk = \sumn=0N-1 fn exp (-2\pi nik / N).
Käänteismuunnos on
fn= (1 / N) \sumk=0N-1 Fk exp (2\pi nik / N).
Normitustekijä 1/N esiintyy vain jommassakummassa muunnoksessa. Joskus sitä käytetään suoran muunnoksen lausekkeessa, jolloin se ei esiinny käänteismuunnoksessa.
Diskreetissä tapauksessa funktion ja sen muunnoksen arvoja esittävät luvut fn ja Fk liittyvät itse funktioiden f ja F arvoihin kaavoilla
fn =
\Delta x f(n\Delta x),
Fk =
F(k \Delta u),
missä
\Delta u \Delta x = 1 / N.
Kun munnettavan vektorin alkioiden määrä N on kakkosen potenssi, muunnos voidaan laskea ajassa N log N.
Danielsonin ja Lanczosin lemma(1942):
Fk =
\sumj=0N-1
fj exp (-2\pi ijk / N)
= \sumj=0N / 2-1
f2j
exp (-2\pi i (2j)k / N)
+ \sumj=0N / 2-1
f2j+1
exp (-2\pi i (2j+1)k / N)
= \sumj=0N / 2-1
f2j
exp (-2\pi i jk / (N /2) )
+ Wk
\sumj=0N /2-1
f2j+1
exp (-2\pi i jk / (N /2))
= Fkparillinen +
Wk
Fkpariton,
missä
W = e -2\pi i / N
Lasketaan siis erikseen parillisten ja parittomien alkioiden muunnokset, ja niistä yhdistetään lopullinen muunnos.
Kumpikin osa voidaan edelleen käsitellä samalla tavalla. Jatketaan rekursiivisesti, kunnes muunnettavana on vain yksi alkio.
program fft
implicit none
real, parameter :: pi=3.141592654
integer, parameter :: m=2, n=2**m
complex a(0:n)
integer i
a(0)=(0.0, 0.0) ; a(1)=(1.0, 0.0)
a(2)=(2.0, 0.0) ; a(3)=(3.0, 0.0)
write (*,*) 'alkuperainen data'
do i=0,n-1
write(*,'(2F8.4)') a(i)
end do
! ensin suora muunnos
call fourier(a, m, n, 1)
write (*,*) 'suora muunnos'
do i=0,n-1
write(*,'(2F8.4)') a(i)
end do
! sitten kaanteismuunnos
call fourier(a, m, n, -1)
write (*,*) 'kaanteismuunnos '
do i=0,n-1
write(*,'(2F8.4)') a(i)
end do
contains
subroutine fourier (a, m, n, sig)
! kompleksimuuttujan Fourier'n muunnos
! alkuperainen data on vektorissa a(0:n-1)
! muunnettavan vektorin pituuden n on oltava
! kakkosen potenssi, n=2**m
! muunnos korvaa alkuperaisen vektorin a(0..n-1)
! jos sig < 0, lasketaan kaanteinen muunnos
implicit none
integer, intent(in) :: m, n, sig
complex, dimension(0:n), intent(inout) :: a
integer i,j,k,l, half_n, le,le1,ip
real angle
complex x,u,w,t
half_n = n/2
! jos kaanteismuunnos, konjugoidaan a
if (sig < 0) a = conjg(a)
\page
! jarjestellaan vektori a s.e. a(s)=a(j), missa
! s on j:n binaariesitys takaperin
j=1
do i=1, n-1
if (i < j) then
x = a(j-1) ; a(j-1) = a(i-1) ; a(i-1) = x ;
end if
k = half_n
do while (k < j)
j = j-k ; k = k/2
end do
j = j+k
end do
! muunnetaan palasia, joiden pituus on
! le=2,4,8,...
le=2
do l=1, m
le1 = le/2
u = cmplx(1,0)
angle = pi/le1
w = cmplx(cos(angle),sin(angle))
do j=1,le1
do i=j,n,le
ip=i+le1
x = a(ip-1)
t = x * u
a(ip-1) = a(i-1)-t
a(i-1) = a(i-1)+t
end do
x = u * w
u = x
end do
le = 2*le
end do
! jos kaanteismuunnos, konjugoidaan a takaisin
! ja normitetaan
if (sig < 0) a=conjg(a)/n
end subroutine
end program
Ohjelman tulostus:
alkuperainen data 0.0000 0.0000 1.0000 0.0000 2.0000 0.0000 3.0000 0.0000 suora muunnos 6.0000 0.0000 -2.0000 -2.0000 -2.0000 0.0000 -2.0000 2.0000 kaanteismuunnos 0.0000 0.0000 1.0000 0.0000 2.0000 0.0000 3.0000 0.0000
Oletetaan yksinkertaisuuden vuoksi, että kuvan leviäminen on samanlaista kaikkialla kuvan alueella. Havaittu kuva g on silloin konvoluutio alkuperäisestä kuvasta f ja kuvausjärjestelmän virhettä kuvaavasta funktiosta h, johon tulee lisäksi kohina n.
g (x, y ) = \int\int h (x-x', y-y') f(x',y') dx' dy' + n (x, y).
Integrointi ulotetaan tässä sellaisen yleensä pienen alueen yli, jonka ulkopuolella h ei enää poikkea merkittävästi nollasta.
Käytännössä kuva on diskreetti pistejoukko, jolloin integraali korvataan äärellisellä summalla
gij = \sumi' \sumj' hi-i', j-j' fi'j' + nij.
Joskus kuvaussysteemin aiheuttama häiriö eli konvoluution ydin h on tunnettu. Tätä funktiota kutsutaan myös nimellä point spread function eli PSF.
Läheskään aina PSF ei ole tiedossa. Mikäli mahdollista, havaintoihin kannattaisi sisällyttää mittaus, josta PSF voidaan selvittää. Sellainen voi olla esimerkiksi hyvin terävä pulssi, pistemäinen valonlähde tai sakara-aalto.
Dekonvoluutio suotimien avulla
Fourier'n muunnos muuttaa konvoluution kertolaskuksi:
G (u, v) = H (u, v) F (u, v) + N (u, v),
missä G, H, F ja N ovat funktioiden g, h, f ja n Fourier'n muunnokset. Tässä esiintyvää funktiota H kutsutaan usein systeemin siirtofunktioksi (transfer function).
Alkuperäinen kuva voidaan ainakin periaatteessa ratkaista:
F (u, v) = G (u, v) / H (u, v) - N (u, v) / H(u, v).
Tyypillisesti h on Gaussin funktio, ja sen Fourier'n muunnos on myös Gaussin funktio. Se poikkeaa oleellisesti nollasta vain hyvin pienellä alueella.
Hieman parempi ratkaisu on käyttää 1/H:n sijasta suodinta
M (u, v) =
1/H (u, v),
jos u2 + v2 < r2,
M (u, v) = 1 muuten,
missä r on niin pieni, että H on varmasti nollasta poikkeava säteen r sisäpuolella.
Pienimmän neliösumman mielessä parhaan tuloksen antaa Wienerin suodin:
M (u, v) = [1 / H (u, v) ] [ |H (u, v)|2 / ( |H (u, v)2 + Snn (u, v) / Sff (u, v) ) ],
missä Snn ja Sff ovat kohinan ja alkuperäisen kuvan spektraalitiheydet. Käytännössä termiä Snn/Sff voi olla vaikea laskea, ja usein se korvataan vakiolla:
M (u, v) = [ 1 / H (u, v)] [ |H (u, v)|2 / ( |H (u, v)|2 + A ) ] = H*(u, v) / ( H*(u, v) H (u, v) + A ).
Vakio A on regularisointiparametri, joka tarkoituksena estää nimittäjää menemästä liian lähelle nollaa. Sopiva arvo löytyy usein helpoimmin kokeilemalla. Hyvin pienet arvot johtavat kyllä terävään kuvaan, mutta samalla vahvistavat kohinaa. Suuret arvot puolestaan tuottavat siistin mutta epätarkan kuvan.