Digitale Signalverarbeitung auf FPGAs · Spektrum abgetasteter Signale (Wdh.) ......

2016

Dr. Christian Münker

Digitale Signalverarbeitungauf FPGAs

DFT: Diskrete Fouriertransformation (Grundlagen)

Teil 1 Die Fourier-Familie

Dr. Christian Münker 12.09.16 ► Digitale Signalverarbeitung auf FPGAs DFT-2 von 34

DFT: Überblick

■ Die Fourier-Familie

■ DFT: Komplex aber nicht kompliziert

■ Frequenzauflösung

■ Schnell, schneller, FFT


Zusätzliches Material zur Wiederholung

■ edx-Kurs „ELEC301x: Discrete Time Signals and Systems“https://courses.edx.org/courses/RiceX/ELEC301x/T1_2014/(Kostenloser Account erforderlich)

■ Videos von Jörg Lovisach:

19.03.1: Diskrete Fourier-Transformation, FFT Teil 1 (www.youtube.com/watch?v=ls5MfqLZVbo)

19.03.2 Diskrete Fourier-Transformation, FFT Teil 2 (www.youtube.com/watch?v=wd6hAPS4ILc)

https://courses.edx.org/courses/RiceX/ELEC301x/T1_2014/

https://www.youtube.com/watch?v=ls5MfqLZVbo

https://www.youtube.com/watch?v=wd6hAPS4ILc


Spektrum abgetasteter Signale (Wdh.)

• =

* =

sin(2π tT 1 )/2+ cos(10π t

T 1

+ φ 0)/4Periodisches Spektrum !

-1

0

1

0 0.2 0.4 0.6 0.8 1

x(t)

→

t / T1 →

Analog Signal

0

1

0 10 20 30 40 50

Ideal Sampler

t / TS →

-1

0

1

0 10 20 30 40 50

x(

nT

S )

≡

x[

n]

→

t / TS ≡ n →

Sampled Signal

0

0.1

0.2

0.3

0.4

0.5

0 2 4 6 f / f1 →

|Sx(f

)|

→

Spectrum Analog Signal

f1

5 f1

0

0.1

0.2

0.3

0.4

0.5

0 1 2 f / f

S →

|Sy(

e j2π

F )|

→ Spectrum Sampled Signal

f1

5 f1

0

1

0 1 2 f / f

S →

Spectrum Ideal Sampler

TS = T

1 / 50

fS = 50 f

1


Anwendungen der DFT

Abschätzung des Spektrums eines analogen oder zeitdiskreten Signals x(t) bzw. x[n] (→ Matlab / Python, digitale Messwerterfassung)

Entwurf zeitdiskreter Filter in der Frequenzebene

Komprimierung zeitdiskreter Signale in der Frequenzebene

Schnelle Faltung (zeitdiskret, { ӿ } ○╱-● { • } )

Robuste Breitband-Datenübertragung über parallele schmale Frequenzbänder (OFDM, PowerLine, LTE, ...)

Kenngrößen bei Fourier-Transformation:

➢ Anzahl N der Zeit- / Frequenzpunkte

➢ Bandbreite Δf je Frequenzband

➢ Minimale Mess- bzw. Einschwingzeit TE

Zeit-Bandbreite-Gesetz: T

E ~ 1 / Δf !

(„Unschärferelation“ der Nachrichtentechnik)


Fourier-Familie: Fourierreihe (CFT*)

Analoges Signal: Spektralschätzung nur mit analogen Methoden möglich (Filterbank, durchstimmbares Bandpassfilter)

X (k ) =1T 1

∫0

T 1

x (t )e− j2π

kT 1

t

dt , k ∈ ℤ

x(t) periodisch-zeitkontinuierlich

X(f) aperiodisch-diskret

* Continuous Fourier Transform


Animation: Fouriereihe

https://commons.wikimedia.org/wiki/File:Fourier_transform_time_and_frequency_domains.gif

neu

Lucas V. Barbosa:

https://commons.wikimedia.org/wiki/File:Fourier_transform_time_and_frequency_domains.gif


Fourier-Familie: Fourierintegral (CFT)x(t) aperiodisch-zeitkontinuierlich

Analoges Signal: Spektralschätzung nur mit analogen Methoden möglich (Filterbank, durchstimmbares Bandpassfilter)

X(f) aperiodisch-kontinuierlich

X (f ) = ∫−∞

∞

x (t )e− j 2π f t dt


Animation: Fourierintegral

https://commons.wikimedia.org/wiki/File%3AContinuous_Fourier_transform_of_rect_and_sinc_functions.gif

neu

Lucas V. Barbosa:

https://commons.wikimedia.org/wiki/File%3AContinuous_Fourier_transform_of_rect_and_sinc_functions.gif


Zeitdiskret, aber keine numerische Berechnung möglich, da ∞ viele Samples. Näherungsweise Darstellung mit DFT und Zero-Padding.

x[n] aperiodisch-zeitdiskret

X(f) periodisch-kontinuierlich

Fourier-Familie: Discrete-Time Fourier Transform

F

X (f ) = ∑n=−∞

∞

x [n]e− j 2π

ff S

n


F

X [k] periodisch-diskret

Fourier-Familie: Discrete Fourier Transform DFT

X [k ] =1N

∑n=0

N−1

x [n]e− j 2π

kN

n

, k = 0 … N−1

x[n] periodisch-zeitdiskret

Zeitdiskret & endliche Anzahl Samples N herausgeschnitten: Für periodisch (fortgesetzte) zeitdiskrete Signale ist numerische Berechnung möglich!


Numerische spektrale Analyse einer Messsequenz (1)

Abtastung des Signals mit fS = 1 / T

S → kein Informationsverlust

Quantisierung des Signals → (hoffentlich) tolerierbarer Informationsverlust

DTFT würde im Bereich 0 … fS identisches Spektrum wie CFT liefern, aber

nicht numerisch berechenbar (∞ viele Samples)

Fensterung = Beschränkung auf N Messwerte: Daten außerhalb der NSamples bzw. Ausschnitts mit Länge T

1 = N

T

S werden zu Null gesetzt

(immer noch DTFT mit ∞ vielen Samples, nicht numerisch berechenbar)

DFT des periodisch mit T1 = N

T

S fortgesetzten Signalausschnitts liefert

identische Spektrumswerte bei k f1 = k f

S / N wie DTFT !

N Datenpunkte x[n] und N Frequenzpunkte X[k] → Block- oder Frame-Transformation → DFT ist numerisch gut berechenbar: Computer!


X2(

f )

F

X3(

f )

Numerische spektrale Analyse einer Messsequenz (2)

CFT

DTFT

DFT

Fensterung ↔ Abtastung

X1(

k

)

X1(k) = X

2,3(kf

S / N) = X

4[k]Wähle N Samples

F

X4 [k]

2016



Kap. 3a Diskrete Fouriertransformation (Grundlagen)

Teil 2 DFT: Komplex aber nicht kompliziert


DFT: Komplex, aber nicht kompliziert

„e t

o th

e pi

tim

es i“

(h

ttps

://x

kcd.

com

/179

/)

von

Ran

dall

Mun

roe

(xkc

d@xk

cd.c

om)

unte

r C

C-B

Y-N

C-2

.5 L

izen

z

https://xkcd.com/179/

mailto:[email protected]

http://creativecommons.org/licenses/by-nc/2.5/


Beispiel für DFT einer reellwertigen Sequenz

N / 2 + 1 Sinusfunktionen (Imaginärteil) →

Aber: k = 0 und N / 2 liefern keinen Beitrag!

N / 2 + 1 Kosinusfunktionen (Realteil) →

Analyse

Synthese

N = 16 Samples x[n]

Σ: N = 16 Koeffizienten X[k] für Basisfunktionen („Frequenzpunkte“)


Analysegleichung der DFT

■ n = 0 … N - 1 ist Laufvariable für N Samples x [n] mit Abstand TS = 1 / f

S .

→ Messfenster mit Länge T1 = NT

S periodisch fortgesetzt, x[n + N ] = x[n ]

■ k = 0 … N - 1 ist Laufvariable für N Frequenzpunkte X [k] zwischen f = 0 … fS ,

periodisch fortgesetzt, X[k + N ] = X[k ] (Wiederholspektren)

➔ fk = f

S k / N [F

k = k / N]

und ∆f = f

S / N = f

1 = 1 / T

1

■ Symmetrie bei reellem x[n]: X[-k ] = X*[ k ] → X[0] und X [ N-N/2

] = X* [

N/2

] reell

➔ Bei reellem x[n] muss nur die Hälfte der Punkte, k ≤ N / 2, berechnet werden

■ X[k] ist i.A. komplex, daher sind 2 Datenspeicher pro Frequenzpunkt k nötig(Real- und Imaginärteile sind Gewichtungsfaktoren für Kosinus- und Sinusfunkt.)

X [k ] =1N

∑n=0

N −1

x [n]e− j2π

knN

⏟Python/Matlab 'fft'

=1N

∑n=0

N−1

x [n ]wNkn mit wN = e

−j2π

N (Drehfaktor)

und k, n = 0 ... N – 1

FrequenzauflösungPhys. Frequenzen


Synthesegleichung der DFT (IDFT)

x [n ] = ∑k=0

N−1

X [k ]e+ j 2π

knN

⏞N⋅Python / Matlab ifft

= ∑n=0

N−1

X [k ]wN−kn

■ N Samples ↔ N Frequenzpunkte ═► Blocktransformation !

■ X[k] im Allgemeinen komplex → N x N komplexe Multiplikationen

■ Aber aufgrund Symmetrie X[-k ] = X*[ k ] für reelle x[n] genügt IDFT über die Hälfte des Spektrums

■ Periodizität der IDFT x[n + N ] = x[n ]: Zeitsignal kann periodisch fortgesetzt werden

■ Wahl des Skalierungsfaktors (hier: 1/N bei DFT) so:

➢ dass nacheinanderfolgende DFT und IDFT von x [n] wieder x [n] ergibt

➢ dass DFT – Werte die Amplituden der Spektralkomponenten korrekt wiedergeben (optional, ist z.B. bei Matlab nicht so)

und w N = e−

j2π

N

mit k, n = 0 ... N – 1


Eigenschaften der DTFT (CFT)

Eigenschaft Zeitbereich x[n] Frequenzbereich X(ejΩ)

Definition x[n] ≡ x(nTS)

Linearität ax[n] + by[n] aX(ejΩ) + bY(ejΩ)

Zeitverschiebung x[n-k], k ∊ ℤ X(ejΩ) e-jΩk

Frequenzverschiebung (Modulation)

x[n] ejnΩ1, Ω1 ∊ ℝ X(ej(Ω-Ω1))

Zeitumkehr x[-n] X(e-jΩ)

Faltung im Zeitbereich x[n] ӿ y[n] X(ejΩ) · Y(ejΩ)

Multiplikation im Zeitbereich x[n] · y[n] X(ejΩ) ӿ Y(ejΩ) / 2π

Inverse DTFT

X (e j Ω) = ∑

n=−∞

∞

x [n ]e- jnΩ

x [n ] = T S ∫−f S

f S

X (e j Ω)ejnΩdΩ


Eigenschaften der DFT (N Punkte)Eigenschaft Zeitbereich x[n] Frequenzbereich X [k]

Definition IDFT / DFT

Periodizität x [n] = x [n - N] X [k] = X [k - N]

Linearität ax [n] + by [n] aX [k] + bY [k]

Zyklische Zeitverschiebung x[n-m], m ∊ ℤ X [k] e-j2πkm / N

Modulation x[n] ej2πnm / N X [k-m]

Zeitumkehr x[-n] = x[N-n] X [-k] = X [N-k]

Faltung im Zeitbereich x[n] ӿ y[n] X [k] · Y [k]

Multiplikation im Zeitbereich x[n] · y[n] X [k] ӿ Y [k] / N

X [k ] =1N

∑n=0

N−1

x [n ]e− j 2π

knNx [n] = ∑

n=0

N−1

X [k ]e j2

knN


Symmetrie der FouriertransformationenBei reellwertigen Zeitsignalen x gelten für CFT / DTFT / DFT x ○-● X

die folgenden Symmetrieeigenschaften des Spektrums:

Reellwertiges x ○-● X*(ejΩ) = X(e-jΩ) (konjugiert-symmetrisches Spektrum)

Re{X} ist gerade

Im{X} ist ungerade

|X| ist gerade (Betragsgang)

∠ X ist ungerade (Phasengang)

gerader Teil von x ○-● Re{X}

ungerader Teil von x ○-● j Im{X}


0 0 0 000 0 0

DFT reellwertiger und komplexer SequenzenR

eell

wer

tig

reell: xR

[n] = x[n]

imaginär: xI [n] = 0

0 N-1

Zeitbereich x[n] mit n = 0 … N -

1 Frequenzbereich X[k]; k = 0 … N

-

1

Ko

mp

lex

0 0

N reellwertige Datenpunkte

N kompl. = 2N reellw. Datenpunkte

reell: xR

[n]

imaginär: xI [n]

N relevante Datenpunkte

reell: XR

[k]

imag.: XI [k]

XR

[N/2 - i ] = X

R [N/2 + i

]

XI [N/2 - i

] = -X

I [N/2 + i

]

2N relevante Datenpunkte0 N/2 N-1

reell: XR

[k]

imag.: XI [k]

DC fS / 2

N/2

fS / N○-●

○-●


Frequenzauflösung der DFTJe größer die Länge N der Datenfolge, desto besser die

Frequenzauflösung der DFT: ∆f = 1 / NTS

Anmerkung: Die höchste Frequenzkomponente der DFT wird natürlich durch die Abtastfrequenz bestimmt: f

k,max = f

N / 2 = f

S / 2

F = k / NΔF = 1 / N


DFT des Dirac-Stoßes (1)

DFT des Dirac-Stoßes: |X[k]| = ? φ[k] = ?

x [n] = δ[n ] ⇒ X [k ] =1N

∑n=0

N−1

δ[n]e− j2π

knN =

1N

x [n] = δ[n−m ] ⇒ X [k ] =1N

∑n=0

N−1

δ[n−m]e− j 2π

knN =

1N

e− j2 π

kmN

→ Alle Basisfunktionen gewichtet mit Faktor 1/N, rein reell.

→ Alle Basisfunktionen gewichtet mit Faktor 1/N,

lineare Phase - 2

π

k

m

/ N .


DFT des Dirac-Stoßes (2)

δ[n - 1]

Zeit Betrag Phase

δ[n]

2016



Kap. 3a Diskrete Fouriertransformation (Grundlagen)

Teil 4 Schnell, schneller ...


DFT als Signal Flow Graph für N = 8

Je N 2 komplexe Multiplikationen und Additionen:

4N 2 reelle Mult.+ 4N 2 reelle Add.

Σ: 4 N 2 MAC- Operationen(N 2 für reellwertige Signale)

N X [k ] = ∑n=0

N−1

x [n ]e− j 2π

knN = ∑

n=0

N−1

x [n]wNkn mit wN

kn= e

− j2πknN (Drehfaktor)

w8

0

w8

0

w8

0

w8

0

w8

0

w8

0

w8

0

w8

0

w8

0

w8

1

w8

2

w8

3

w8

4

w8

5

w8

6

w8

7

w8

0

w8

2

w8

4

w8

6

w8

8=

w8

0

w8

2

w8

4

w8

6

w8

0

w8

3

w8

6

w8

9=

w8

1

w8

4

w8

7

w8

2

w8

5

w8

0

w8

4

w8

0

w8

4

w8

0

w8

4

w8

0

w8

4

w8

0

w8

5

w8

2

w8

7

w8

4

w8

1

w8

6

w8

3

w8

0

w8

6

w8

4

w8

2

w8

0

w8

6

w8

4

w8

2

w8

0

w8

7

w8

6

w8

5

w8

4

w8

3

w8

2

w8

1

N X [1]

N X

[2]

N X

[3]

N X

[7]

N X

[0]

N X [4]

N X

[5]

N X

[6]

x [1]x

[0] x

[2] x

[3] x

[4] x

[5] x

[6] x

[7]

N X = W N x

Matrix-schreibweise:

n →

k →


DFT als FIR-Filterbank

DFT kann auch interpretiert werden als Filterbank aus N FIR-Filtern mit z.T. komplexwertigen Koeffizienten:

N X [k ] = ∑n=0

N−1

x [n ]e− j 2π

knN = ∑

n=0

N−1

x [n]wNkn mit wN

kn= e

− j2πknN

N X

[0]

x [1]x

[0] x

[2] x

[3] x

[4] x

[5] x

[6] x

[7]

MA-Tiefpass → DC

N X

[2]

−1

−1+1+1 −

j−

j+ j+ j

−1

N X [4]

−1

−1

−1

+1 +1 +1 +1

Bandpass → fS /4

MA-Hochpass → fS /2


Aufteilung in zwei N/2-Punkte DFTs

Aufteilung in N/2 – Punkte DFTs für ungerade und gerade Samples und geeignete Kombination der Teilergebnisse reduziert Gesamt-MACs

Umsortierung der Eingangssamples notwendig

Herleitung nicht trivial, siehe Literatur oder 9. Kapitel (vielleicht)

[„D

ecim

atio

n in

Tim

e“] N X

[1]

N X

[2]

N X

[3]

N X

[7]

N X

[0]

N X [4]

N X

[5]

N X

[6]

x [2]

x [4]

x [8]

x [7]

x [0]

x [1]

x [3]

x [5]

w 83

w 82

w 81

w 80

N/2 – PunkteDFT:

N²/4 MACs

N/2 – PunkteDFT:

N²/4 MACs


FFT: Sukzessive Aufteilung in Teil-DFTs

ca.

2*

N lo

g2

N r

eel

lw. M

ult

. (a

bh

äng

ig v

on

A

lgo

rith

mu

s),

vg

l. D

FT

: N

2 r

eellw

. Mu

lt.

N X [1]

N X

[2]

N X

[3]

N X

[7]

N X

[0]

N X [4]

N X

[5]

N X

[6]

x [4]

x [2]

x [6]

x [7]

x [0]

x [1]

x [5]

x [3]

w 80 w 8

2

w 80

w 83

w 82

w 80 w 8

1

w 80

w 80 w 8

2

w 80

w 80

log2 N

* „2N“, da kaum Einspareffekte für reellwertige Eingangssignale !

2-Pkte DFTs Kombination Teilergebnisse

N


Symmetrien und „Butterfly“ - Operation

w Nk+8/2

w Nk

w Nk

wNN= e

− j2πN

N = 1; w NN /2

= e− j

πNN =−1

wNN +k = w N

N w Nk = w N

k (Periodizität )

wNN−m = wN

−m = wNm ∗

wNN /4

=− j; wN3N /4

= j

w Nk+ N / 2

= w Nk w N

N / 2= −w N

k

wNk= e

−j 2k π

N

Nutze Symmetrien des komplexenDrehfaktors („twiddle factor“) zur Vereinfachung der Rechnung:

Damit Optimierung der 2-Punkte DFT („butterfly computation“)


„Optimale“ FFT ?

■ Für NFFT

= 2k, lässt sich die FFT so zerlegen, dass man die

minimale Anzahl Rechenoperationen erhält

■ Aber: Verschiedene sehr gute FFT-Algorithmen wenn eine Zerlegung von N

FFT in kleine Primfaktoren möglich ist

■ Für diese Fälle ist NMAC

= O(N log2 N)

■ Weitere Optimierungsziele:

■ Minimaler Speicherverbrauch

■ minimale Anzahl von Sortiervorgängen

■ Eingangssamples werden in korrekter Reihenfolge verarbeitet

■ ...

Siehe auch http://lighthouseinthesky.blogspot.de/2010/03/flops-and-fft.html

http://lighthouseinthesky.blogspot.de/2010/03/flops-and-fft.html


FFT in Python / MatlabFFT ist in Python / Matlab schnell aufgestellt:

Konstanten: NFFT = 150; fS = 1000

Signal: n = np.arange(NFFT); y = np.sin(2*pi*n/16)

FFT (skaliert): Y = fft(y)/NFFT

Frequenzen: f = fftfreq(NFFT, 1/fS) # Angabe von fs optional

f = np.arange(NFFT)* fS/NFFT # Alternativ

Zweiseit. Spektrum: Y2 = fftshift(Y)

f2 = fftshift(f)

Plotten: plot(f2, np.abs(Y2))

Plotten bis fS /2: plot(f[0:NFFT/2], np.abs(Y[0:NFFT/2])

Python: import fft, ... from numpy.fft

Diese Folien und die zugehörigen Videos sind unter Creative-Commons-Lizenz CC-BY-NC-SA 3.0 de veröffentlicht.

Bei Verwendung dieses Werks müssen Sie auf die entsprechende CC-Lizenzurkunde verweisen, in diesem Fall http://creativecommons.org/licenses/by-nc-sa/3.0/de/ .

Sie müssen ferner die folgenden Angaben machen ("BY“, attribution)

Author („Christian Münker“)

Titel („Digitale Signalverarbeitung auf FPGAs“)

URL zu Werk (https://github.com/chipmuenk/dsp_fpga) und / oder Author ( http://www.chipmuenk.de)

Außerdem ist die Verwendung auf folgende Weise eingeschränkt:

Diese Materialien dürfen nur nicht kommerziell genutzt werden („NC“, non-commercial).

Dieses Werk oder Teile daraus dürfen nur unter gleichen Lizenzbedingungen weiterverteilt werden („SA“, share alike).

Fragen, Anmerkungen, Anregungen, Bugs, Bierbons bitte an [email protected] wünsche viel Erfolg und Spaß (?!) mit den Materialien!

http://creativecommons.org/licenses/by-nc-sa/3.0/de/

https://github.com/chipmuenk/dsp_fpga

http://www.chipmuenk.de/

mailto:[email protected]

Digitale Signalverarbeitung auf FPGAs · Spektrum abgetasteter Signale (Wdh.) ......

Documents

Transcript of Digitale Signalverarbeitung auf FPGAs · Spektrum abgetasteter Signale (Wdh.) ......