DIGITAL SIGNALBEHANDLING

Komplement till
DIGITAL SIGNALBEHANDLING
Fredrik Gustafsson
2005
2
Dokumentet inneh˚
aller nyskriva delar till nya upplagan.
Fredrik
Inneh˚
all
8.2
2.1
3.1
Station¨
ara kalmanfiltret . . . . . . . . . . . . . . . .
8.2.1 H¨
arledning av station¨
ara kalmanfiltret . . . .
Mikromodul: stokastiska processer . . . . . . . . . .
2.1.1 Stokastiska variabler . . . . . . . . . . . . . .
2.1.2 Stokastiska vektorer och kovarians . . . . . .
2.1.3 Stokastiska processer och kovariansfunktionen
Hilberttransformen . . . . . . . . . . . . . . . . . . .
3.1.1 Inledning . . . . . . . . . . . . . . . . . . . .
3.1.2 Modulering . . . . . . . . . . . . . . . . . . .
3.1.3 Hilberttransformen . . . . . . . . . . . . . . .
3.1.4 Utv¨
ardering . . . . . . . . . . . . . . . . . . .
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
7
9
9
9
11
15
15
15
17
20
8.2 Station¨ara kalmanfiltret
5
8.2 Station¨
ara kalmanfiltret
Tillst˚
andsmodell f¨
or m¨
atningar t = 1, 2, . . . , N :
x(t + 1) =Ax(t) + Bw(t),
y(t) =Cx(t) + v(t),
Cov(w(t)) = Q,
Cov(v(t)) = R,
E(x(1|0)) =x0 ,
Cov(x(1|0)) = P0 .
Utan information ger (8.3.a) att initialskattningen och dess os¨
akerhet
ska propageras enligt
x
ˆ(t|t − 1) = Aˆ
x(t − 1|t − 1), P (t|t − 1) = AP (t − 1|t − 1)AT + BQB T .
(8.2a)
Vi ans¨
atter nu ett godtycklig linj¨
ar funktion av informationen i senaste
m¨
atningen f¨
or att korrigera skattningen
a
(8.3a)
Vet vi skattningen x
ˆ(t − 1|t − 1) vid tiden t − 1, blir prediktionen
(s¨
att w(t − 1) = 0 i signalmodellen):
x
ˆ(t|t − 1) = Aˆ
x(t − 1|t − 1)
Prediktionsfelet, eller innovationen som det ocks˚
a kallas i kalmanfiltret,
ges av
ε(t) = y(t) − C x
ˆ(t|t − 1)
Vi ans¨
atter nu ett linj¨
art filter
x
ˆ(t|t) = x
ˆ(t|t − 1) + Ly(t) + M.
F¨or att bevara ett korrekt v¨
antev¨
arde m˚
aste 0 = E(Ly(t) + M ) =
E(L(Cx(t) + v(t)) + M ) = LCx(t) + M vilket kan erh˚
allas genom att
v¨
alja konstanttermen till M = −C x
ˆ(t|t − 1). Vi f˚
ar d˚
a filtret
x
ˆ(t|t) = x
ˆ(t|t − 1) + L(y(t) − C x
ˆ(t|t − 1)).
Skattningsfelet definieras som
x
˜(t + 1|t) = x(t + 1) − x
ˆ(t + 1|t),
6
Inneh˚
all
och p˚
a samma s¨
att f¨
or x
˜(t|t) = x(t) − x
ˆ(t|t). Skattningsfelet ges av
rekursionen
x
˜(t + 1|t) = x(t + 1) − x
ˆ(t + 1|t)
(8.4a)
= Ax(t) + Bw(t) − A (ˆ
x(t|t − 1) + L(y(t) − C x
ˆ(t|t − 1)))
(8.4b)
= Ax(t) + Bw(t) − A (ˆ
x(t|t − 1) + L(Cx(t) + v(t) − C x
ˆ(t|t − 1)))
(8.4c)
= (A − ALC)˜
x(t|t − 1) + Bw(t) − ALv(t)
(8.4d)
Felrekursionen kan expanderas till formen
x
˜(t + 1|t) = (A − ALC)t+1 (x(0) − x0 ) +
t
X
(A − ALC)t−k (w(k) − ALv(k))
k=0
Vi noterar f¨
oljande:
• F¨
or att initalfelet (x(0) − x0 ) ska avta mot noll, s˚
a m˚
aste K
v¨
aljas s˚
a att egenv¨
ardena till A − ALC ¨
ar innanf¨
or enhetscirkeln
(⇒ kA − ALCkt → 0, t → ∞).
• Valet av L ¨
ar f¨
or ¨
ovrigt en kompromiss mellan att minimera
initialfelets inverkan och att undertrycka b˚
ade processbrus och
m¨
atbrus.
M˚
alet kan nu uttryckas i att v¨
alja L s˚
a att felets kovariansmatris
Cov˜
x(t + 1|t) =(A − ALC)t P0 ((A − ALC)T )t+1 +
t
X
k=0
(A − ALC)t−k Q − ALRLT AT ((A − ALC)T )t−k
minimeras. N˚
agra specialfall att beakta:
1. Observat¨
or. Finns inget brus ¨
ar det endast transienten av initialfelet som beaktas. Observat¨
orens poler, dvs egenv¨
ardena till
A − ALC kan v¨
aljas godtyckliga om tillst˚
andsmodellen ¨
ar observerbar.
2. Station¨
ara kalmanfiltret fokuserar p˚
a att balansera summan av
brustermerna, eftersom att station¨
art har initialtillst˚
andets inverkan f¨
orsvunnit.
8.2 Station¨ara kalmanfiltret
7
3. Tidsvariabla kalmanfiltret balanserar inverkan av alla tre termerna s˚
a att skattningsfelets kovariansmatris minimeras.
8.2.1 H¨
arledning av station¨
ara kalmanfiltret
Kovariansmatrisen f¨
or v¨
anster- respektiva h¨
ogerled av (8.4) ger en bra
utg˚
angspunkt f¨
or att h¨
arleda kalmanfiltret,
P (t + 1|t) = = (A − ALC)P (t|t − 1)(A − ALC)T + BQB T − ALRLT AT .
(8.5)
Kvadratkomplettering ger (med den kompaktare notationen P = P (t|t−
1))
P (t + 1|t) =AP AT − ALCP C T LT AT + ALCP AT
+ AP C T LT AT + BQB T − ALRLT AT
=A P + L(CP C T + R)LT + LCP + P C T LT AT + BQB T
=A P + L − P C T (CP C T + R)−1 (CP C T + R) ×
T T
L − P C T (CP C T + R)−1
A − P C T (CP C T + R)−1 CP + BQB T
≥A P − P C T (CP C T + R)−1 CP AT + BQB T .
Det L som minimerar kovariansmatrisen P (t + 1|t) och dess minimum
ges av
L(t) = P C T (CP (t|t − 1)C T + R)−1 ,
T
T
−1
P (t + 1|t) = A P (t|t − 1) − P C (CP (t|t − 1)C + R)
(8.6a)
CP AT + BQB T .
(8.6b)
8
Inneh˚
all
2.1 Mikromodul: stokastiska processer
9
2.1 Mikromodul: stokastiska processer
2.1.1 Stokastiska variabler
En stokastiskt variabel X beskrivs av dess t¨
athetsfunktion pX (x), vars
viktigaste egenskaper sammanfattas i dess f¨
orsta tv˚
a s.k. moment:
medelv¨
arde och varians
Z
E(X) = xpX (x)dx
Zx
2
2
Var(X) =
x − E(X) pX (x)dx = E(X 2 ) − E(X)
x
Normalf¨
ordelningen har
(X−µ)2
1
pX (x; µ, σ 2 ) = √
e− 2σ2 ,
2πσ 2
athetsfunktionen.
d¨
ar parametrarna µ = E(X) och σ 2 = Var(X) beskriver hela t¨
Ett konfidensintervall med konfidensgrad 95% ges f¨
or normalf¨
ordelingen
av µ ± 1.96σ, dvs. Prob(|X − µ| < 1.96σ) = 0.95.
I Matlab genereras ett normalf¨
ordelat slumptal av x=mu+sigma*randn(1,1).
En vektor eller matris skapas av oberoende slumptal skapas av randn(n,m).
Konfidensgraden f¨
or ett intervall kan verifieras numeriskt med
>> x=randn(10000,1);
>> length(find(x<1.96 & x>-1.96))/10000
ans =
0.9512
dvs. 95% av alla slumptal hamnar i intervallet [−1.96, 1.96]. Analytiskt
kan sannolikheten Prob(|X| > 1.96) r¨
aknas ut med
>> erfc(1.96)
ans =
0.0056
2.1.2 Stokastiska vektorer och kovarians
En vektor X (dimension n × 1) av stokastiska variabler beskrivs p˚
a
analogt s¨
att av en vektorv¨
ard t¨
athetsfunktion, t.ex. den multivariabla
10
Inneh˚
all
normalf¨
ordelningen
pX (x; µ, P ) = p
1
(2π)n
T P −1 (X−µ)
det(P )
e−0.5(X−µ)
.
V¨antev¨
arde µ (n × 1) och kovarians P (n × n) definieras av
Z
µ = E(X) = xpX (x)dx
x
P = Cov(X) = E (x − µX )(x − µX )T
Z
= (x − µX )(x − µX )T pX (x)dx.
x
Kovariansen anger hur starkt kopplade slumptalen ¨
ar till varandra. Om
element Pij = 0 betyder det att elementen xi och xj ¨
ar okorrelerade.
Omv¨
ant, om Pij > 0 s˚
a¨
ar xi och xj positivt korrelerade s˚
a ett om den
ena ¨
ar positiv s˚
a¨
ar den andra det sannolikt ocks˚
a.
Konfidensintervall kan enkelt generaliseras till vektorer. Ett omr˚
ade
med konfidensgrad 95% ges f¨
or normalf¨
ordelningen av {X : (X −
µ)T P −1 (X − µ) < 1.96}.
En linj¨
artransformation Y = AX ger en ny stokastisk vektor Y (som
kan vara l¨
angre eller kortare ¨
an X. Fr˚
an definitionen h¨
arleds enkelt
f¨
oljande r¨
akneregler f¨
or de f¨
orsta momenten:
Y = AX ⇒
µY = E(Y ) = AE(X) = AµX
PY = Cov(Y ) = ACov(X)AT = APX AT .
En normalf¨
ordelad vektor av slumptal kan genereras i Matlab med
[U,S,V]=svd(P);
Psqrt=U*sqrt(S);
x=mu+Psqrt*randn(n,1);
De f¨
orsta tv˚
a raderna r¨
aknar ut en matriskvadratrot av P som definieras
av P 1/2 P T /2 = P . SVD faktoriserar en godtycklig matris P = U SV T ,
d¨
ar U ¨
ar unit¨
ar (dvs. U T U = U U T = I) och S ¨
ar diagonal. I detta
fall n¨
ar P ¨
ar symmetrisk och positivt semidefinit g¨
aller att U = V
an dessa egenskaper f¨
oljer att kvadratrotsdefinitionen
och Sii ≥ 0. Fr˚
ar uppfylld med konstruktionen P 1/2 = U S 1/2 .
¨
2.1 Mikromodul: stokastiska processer
11
F¨
oljande exempel illustrerar principerna f¨
or linj¨
ara avbildningar.
200 slumptal X fr˚
an standardnormalf¨
ordelningen genereras f¨
orst, och
avbildas sedan linj¨
art p˚
a 200 stokastiska vektorer Y .
X=randn(2,200);
A=[1 0.5;0 1];
Y=A*X;
figure(1)
plot(X(1,:),X(2,:),’.’)
figure(2)
plot(Y(1,:),Y(2,:),’.’)
Kovariansmatrisen ges av
P = AIAT =
1.25 0.5
,
0.5
1
vilken illustreras av niv˚
akurvan i figur 2.1.b
H¨
arur f¨
oljer t.ex. att Var(y1 ) = P11 = 1.25 och E(y1 y2 ) = 0.5.
Elementen y1 och y2 ¨
ar allts˚
a positivt korrelerade, s˚
a om den ena ¨
ar
positiv s˚
a¨
ar den andra det ocks˚
a med stor sannolikhet.
Den ellips som svarar mot 95% av alla punkter ritas med nedanst˚
aende
ber¨
akningar, d¨
ar enhetscirkeln avbildas av samma transformation y =
Ax och sedan skalas till ¨
onskad konfidensgrad.
P=A*A’
phi=linspace(0,2*pi,100);
x=[cos(phi);sin(phi)];
y=A*x;
plot(1.96*y(1,:),1.96*y(2,:),’-’)
Notera att den positiva korrelationen mellan y1 och y2 kan ses p˚
a
ellipsens form.
2.1.3 Stokastiska processer och kovariansfunktionen
Om x[k] ¨
ar en sekvens av stokastiska variabler ¨
ar det en stokastisk
process. Den stora skillnaden p˚
a en stokastisk vektor och en stokastisk
process ¨
ar att den senare oftast anses vara o¨
andligdimensionell. I ¨
ovrigt
12
Inneh˚
all
3
2.5
2
2
1.5
1
1
2
0
y
x
2
0.5
0
−0.5
−1
−1
−1.5
−2
−2
−2.5
−3
−3
−2
−1
0
1
2
3
−3
−2
x1
(a)
−1
0
y1
1
2
3
4
(b)
Figur 2.1 Tv˚
adimensionella normalf¨
ordelade slumptal med v¨
antev¨
arde noll och
kovariansmatris P = I respektive P = AAT . Ett konfidensomr˚
ade svarande mot
95% av alla punkter ¨
ar ocks˚
a markerat.
definieras v¨
antev¨
ardesfunktionen och kovariansfunktionen (notera ordet funktion i dessa begrepp) analogt:
mx [k] = E(x[k])
Rxx [k, l] = E (x[k] − mx [k])(x[l] − mx [l])
Ser vi x som en vektor med element x[k], ges v¨
antev¨
ardet av en vektorn
med element µk = m[k] och kovariansen av matrisen P med element
Pkl = Rxx [k, l].
Vi fokuserar p˚
a processer med medelv¨
arde 0 (mx [k] = 0) som ¨
ar
station¨
ara. Station¨
aritet definieras av att Rxx [k, l] = Rxx [k − l], dvs.
endast tidsskillnaden och inte absolut tid inverkar p˚
a kovariansfunktionen, och skriver
Rxx [l] = E(x[k]x[k − l]).
Om x tolkas som en vektor av l¨
angd N blir kovariansmatrisen


Rxx [0]
Rxx [1]
. . . Rxx [N − 1]
 Rxx [−1]
Rxx [0]
. . . Rxx [N − 2]


Cov(x) = 
.
..
..
.
.


.
.
.
Rxx [−N + 1] Rxx [−N + 2] . . .
Rxx [0]
P˚
a grund av symmetrin i matrisen (kallas f¨
or Toeplitz-struktur), samt
symmetrin Rxx [−k] = Rxx [k], r¨
acker det med att specifiera f¨
orsta
2.1 Mikromodul: stokastiska processer
13
raden i denna matris. Notera att f¨
orsta raden sammanfaller med kovariansfunktionen f¨
or den station¨
ara processen.
Korskovariansfunktionen definieras analogt som
Rxy [l] = E(x[k]y[k − l]).
I motsats till kovariansfunktionen som alltid ¨
ar symmetrisk g¨
aller h¨
ar
att Ryx [−l] = Rxy [l] 6= Rxy [−l] = Ryx [l].
Ett typexempel p˚
a station¨
ar stokastisk process ¨
ar filtrerat vitt brus.
Vitt brus svarar mot en sekvens av oberoende slumptal x[k] med
mx [k] = 0 och Rxx [k] = δ[k].
Exempel: en sekvens y[k] genereras fr˚
an vitt brus x[k] med f¨
oljande
rekursion
y[k] = 0.8y[k − 1] + x[k].
Rekursionen kan utvecklas till
y[k] = x[k] + 0.8x[k − 1] + . . . 0.8l x[k − l] + . . .
Kausalitet ger att y[k] ¨
ar oberoende av framtida x[k + l], s˚
a att korskovariansfunktionen blir
(
0.8l , l ≥ 0,
Ryx [l] = E(y[k]x[k − l]) =
0, l < 0.
P˚
a motsvarande s¨
att blir kovariansfunktionen for y[k]
Ryy [l] = E(y[k]y[k − l])
= E (x[k] + 0.8x[k − 1] + . . . x[k − l] + 0.8x[k − l − 1] + . . . )
∞
X
1
25
|l|
= 0.8
(0.82 )i = 0.8|l|
= 0.8|l|
≈ 2.78 · 0.8|l|
1 − 0.82
9
i=1
I Matlab kan exemplet testas med
x=randn(10000,1);
y=filter([1 0],[1 -0.8],x);
Ryy=covf(y,20)
k=0:19;
plot(k,Ryy,’-’,k,25/9*0.8.^k,’--’)
vilket genererar Figur 2.2.
14
Inneh˚
all
3
Skattad
Sann
2.5
Ryy[k]
2
1.5
1
0.5
0
0
5
10
k
15
Figur 2.2 Skattad och sann kovariansfunktion Ryy [k].
20
3.1 Hilberttransformen
15
3.1 Hilberttransformen
3.1.1 Inledning
Att best¨
amma frekvensen f¨or en sinusformad signal vore enkelt med
en komplex signalrepresentation, d˚
a signalen skulle ges av
x[k] = ei2πf kT .
(3.7)
Genom att utnyttja argumentfunktionen f¨
or komplexa tal kan man
studera differensen mellan signalens argument och f˚
ar
arg x[k] − arg x[k − 1] = 2πf kT − 2πf (k − 1)T = 2πf T.
(3.8)
Rent allm¨
ant brukar man definiera en signals momentana frekvens som
f [k] =
arg x[k] − arg x[k − 1]
,
2πT
(3.9)
och den allm¨
anna algoritmen 3.8 kan allts˚
a alltid appliceras p˚
a komplexa signaler f¨
or att ber¨
akna momentan frekvens. Fr˚
agan ¨
ar nu hur
man g¨
or f¨
or reellv¨
arda signaler? Det finns tv˚
a trick f¨
or detta:
• L˚
atsas att man har observerat realdelen av en s.k. analytisk signal, och fr˚
an denna ˚
aterskapa imagin¨
ardelen.
• Modulera signalen p˚
a en komplex b¨
arv˚
ag som i telekommunikation.
B˚
ada metoderna ¨
amnar att stryka det ena sidbandet f¨
or signalen och
beh˚
alla det andra, och p˚
a detta s¨
att f˚
a en komplex signal. T.ex. har
en ren sinus x[k] = cos(i2πf0 kT ) Fouriertransformen X(f ) = δ(f −
f0 ) + δ(f + f0 ), dvs. den har en spegelfrekvens i −f0 . Om vi kan stryka
det negativa sidbandet f¨
orsvinner denna, och vi har transformparet:
i2πf
kT
0
x[k] = e
↔ X(f ) = δ(f − f0 ).
3.1.2 Modulering
L˚
at x[k] vara den observerade signalen, och bilda
xc [k] = x[k]ei2πfc kT ,
(3.10)
arv˚
agsfrekvensen. Frekvensinneh˚
allet f¨
or x[k] kommer att
d¨
ar fc a
¨r b¨
flyttas fr˚
an att ha varit centrerat kring f = 0 till centrering kring
16
Inneh˚
all
Impulssvar för diskret LP−filter
0.6
0.5
0.4
0.3
0.2
0.1
0
−0.1
0
2
4
6
8
10
12
14
16
18
20
Figur 3.3 Impulssvar f¨
or hilbertfilter
b¨arv˚
agen f = fc . Vi har nu en komplex representation av signalen
som kan anv¨
andas f¨
or frekvensskattning. N˚
agra extrasteg f¨
orenklar
implementeringen:
1. V¨
alj fc = fs /4 och modulera xc [k] = x[k]ei2πfc kT .
ansfrekvens fc s˚
a att endast
2. L˚
agpassfiltrera xcf = h∗xc [k] med gr¨
det negativa sidbandet f¨
or X(f ) beh˚
alls. LP-filtret ges i figur 3.4.
3. Frekvenserna blir nu spegelv¨
anda, s˚
a momentana frekvensen som
kommer ut fr˚
an (3.8) beh¨
over reflekteras med fs /4 − f och vi f˚
ar
fs arg xcf [k] − arg xcf [k − 1]
−
.
fˆ[k] =
4
2πT
(3.11)
I denna algoritm beh˚
alls det negativa sidbandet f¨
or signalen enligt
figur 3.3.
3.1 Hilberttransformen
17
Illustration
1.4
1.2
Original
Modulerat
LP−filter
1
H(f)
0.8
0.6
0.4
0.2
0
−fs/2
−fs/4
0
f
fs/4
fs/2
Figur 3.4 Impulssvar f¨
or LP-filtret
3.1.3 Hilberttransformen
Matematisk bakgrund
Cauchy-Riemanns relation: Om f (z) = u(x, y) + iv(x, y) ¨
ar deriverbar
i en punkt z0 = x0 + iy0 , s˚
a g¨
aller f¨
oljande samband:
ux (x0 , y0 ) = vy (x0 , y0 )
(3.12)
uy (x0 , y0 ) = −vx (x0 , y0 )
(3.13)
Cauchys integralsats s¨
ager att om en komplex funktion ¨
ar entydigt
best¨
amd av sitt v¨
arde p˚
a randen av en region d¨
ar den ¨
ar analytisk.
Dessa samband g¨
or att vi entydigt kan best¨
amma en analytisk funktion
enbart genom att k¨
anna dess real- eller imagin¨
ar-del p˚
a randen av dess
definitionsomr˚
ade. Detta kallas i matematiken f¨
or Poissons formel.
Om dessa b˚
ada samband ist¨
allet anv¨
ands i signalbehandling och
appliceras p˚
a den kontinuerliga fouriertransformen erh˚
alls den s.k.
hilberttransformen. En signal x(t) kallas analytisk om dess fouriertransform ¨
ar kausal, dvs X(iω) = 0 f¨
or ω < 0. F¨
or en diskret sekvens
18
Inneh˚
all
x[k] definierar vi en analytisk signal genom att kr¨
ava att DFT:n uppjωn
fyller X(e ) = 0 f¨
or −π ≤ ω < 0, dvs. den ¨
ar noll p˚
a undre delen av
enhetscirkeln.
Allts˚
a, vi m¨
ater xr [k], och vi vill skapa en analytisk signal x[k] =
a att
xr [k] + xi [k] s˚
X(ejω ) = T DF T [x[k]] = 0,
−π ≤ ω < 0.
Fr˚
an grundl¨
aggande transformsamband g¨
aller att en reell sekvens har
komplexsymmetrisk transform, och samma sak f¨
or inverstransformen.
Det g¨
aller d˚
a att
Xr (ejω ) = 0.5 X(ejω ) + X ∗ (ejω )
iXi (ejω ) = 0.5 X(ejω ) − X ∗ (ejω )
Fr˚
an dessa samband h¨
arleder vi
−iXr (ejω ), 0 ≤ ω < π
Xi (ejωn ) =
iXr (ejω ), π ≤ ω < 0
eller p˚
a filterform Xi (ejω ) = H(ejω )Xr (ejω ) d¨
ar
H(ejω ) =
−i, 0 ≤ ω < π
i, π ≤ ω < 0
Matlabs funktion hilbert opererar direkt i FFT:n X[n] genom att
nollst¨
alla alla spegelfrekvenser. F¨
or j¨
amna N , blir Xc [n] = X[n]H[n]
d¨
ar

n = 0, N/2
 1,
2, n = 1, 2, . . . , , N/2 − 1
H[n] =

0, n = N/2 + 1, . . . , N − 1
Detta svarar mot ett diskret filter h[n] enligt figur 3.5. Eftersom impulssvaret avtar som 1/n ¨
ar det inte helt enkelt att implementera, ¨
aven
om varannan koefficient ¨
ar noll.
Numeriskt exempel
Ett exempel p˚
a frekvensskattning d¨
ar f0 = 1/6 ges nedan:
3.1 Hilberttransformen
19
Impulssvar för diskreta hilbertfiltret
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
0
5
10
15
20
25
30
35
Figur 3.5 Impulssvar f¨
or hilbertfilter
>> x=sin(2*pi*[0:5]/6)
x =
0
0.8660
0.8660
0.0000
-0.8660
-0.8660
>> xc=hilbert(x)
xc =
0-1.0i 0.8660-0.5i 0.8660+0.5i 0.0+1.0i -0.8660+0.5i
>> fft(x)
ans =
0.0 -0.0-3.0i
0.0
-0.0
0.0 -0.0+3.0i
>> fft(x)
ans =
0.0 -0.0-3.0i 0.0 -0.0 0.0 -0.0+3.0i
>> diff(angle(xc))/2/pi
ans =
0.1667
0.1667
0.1667
0.1667
-0.8333
Notera att ¨
ovre frekvensindex blir noll i Hilberttransformen, och att
en moduloopereration som skiftar in frekvenserna i intervallet [0,0.5]
-0.8660-0.5i
20
Inneh˚
all
Frekvensskattning med hilberttransform, N=30
Modulationsfrekvensskattning, N=30
AR(2) frekvensskattning, N=30
0.7
0.5
0.45
Referens
SNR = 10
SNR = 100
0.6
0.5
Referens
SNR = 10
SNR = 100
0.45
0.35
0.3
0.25
0.2
0.15
Skattad frekvens (normaliserad)
0.4
Skattad frekvens (normaliserad)
Skattad frekvens (normaliserad)
0.4
Referens
SNR = 10
SNR = 100
0.5
0.4
0.3
0.2
0.1
0.35
0.3
0.25
0.2
0.15
0.1
0.1
0.05
0
0
0.05
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Sann frekvens (normaliserad)
0.4
0.45
0.5
(a)
0
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Sann frekvens (normaliserad)
0.4
0.45
0.5
0
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Sann frekvens (normaliserad)
0.4
0.45
(b)
(c)
Figur 3.6 J¨
amf¨
orelse av de grundmetoder f¨
or frekvensskattning: baserat p˚
a
frekvensmodulering och fasdifferenser, baserat p˚
a hilberttransformering och fasdifferenser, och modellbaserad AR(2)-metod.
beh¨
ovs.
3.1.4 Utv¨
ardering
Vi ska h¨
ar studera hur de tv˚
a icke-parametriska metoder f¨
orm˚
ar skatta
en ren sinus fr˚
an brusiga m¨
atningar.
Figur 3.6 sammanfattar de numeriska resultaten:
• Modulationsbaserat skattning sker rekursivt och med cirka 16
multiplikationer och additioner per sampel f¨
or sj¨
alva moduleringen. Det tillkommer operationer f¨
or arctan och division i argumentber¨
akningarna f¨
or varje sampel, samt medelv¨
ardesbildning.
Efter en transient om 20 sampel f˚
as en fasdifferensskattning av
frekvens utan bias och med l˚
ag varians.
• Hilberttransformen implementerad i frekvensdom¨
anen arbetar
segmentvis p˚
a data och kr¨
aver d˚
a log2 (N ) operationer per sampel
i snitt. Det tillkommer operationer f¨
or arctan och division i argumentber¨
akningarna f¨
or varje sampel, samt medelv¨
ardesbildning.
Skattningen f¨
orefaller ha en viss bias f¨
or l˚
aga frekvenser med f˚
a
data, men variansen ¨
ar liten.
• AR-skattningen byggs upp av normalekvationerna i minstakvadratmetoden, och dessa kr¨
aver cirka 5 operationer per sampel.
Skattningen bildas med cirka 6 operationer samt l¨
osning av andragradsekvation och argumentber¨
akning, men beh¨
over inte ske
0.5
3.1 Hilberttransformen
21
varje sampel. AR-modellen modellerar inte en sinus i brus p˚
a
korrekt s¨
att s˚
a d¨
arf¨
or erh˚
alls en bias, och variansen ¨
ar st¨
orre ¨
an
f¨
or de b¨
agga andra metoderna.