Home

Analyse des Orthogonal Frequency Division Multiplexing (OFDM) bei Digital Radio Mondiale (DRM)

Auf der Nacht der Wissenschaften habe ich in Erlangen im Fraunhofer Institut für Integrierte Schaltungen am Lehrstuhl für Informationstechnik mit dem Schwerpunkt Kommunikationselektronik (LIKE) eine Vorführung des digitalen Rundfunks über Kurzwelle mit DRM gesehen. Dort wurde ein Kurzwellensignal durch einen einfachen Mischer auf eine Zwischenfrequenz von 12kHz umgesetzt und direkt an den Soundeingang eines PC angeschlossen. Die Demodulation und Dekodierung lief als Software Defined Radio (SDR) auf dem PC.

Das OFDM wird in vielen aktuellen Kommunikationssystemen wie DVB-T, ADSL und WLAN (siehe OFDM bei Wikipedia) eingesetzt. Um dieses Verfahren besser zu verstehen, wollte ich es mir selbst näher anschauen. Die meisten Signale sind aber viel zu breitbandig, um sie mit "Hausmitteln" zu analysieren. Daher habe ich beschlossen, auf meinem PC die DRM-Signale zu analysieren.

Kurzwellen-Empfänger

Das DRM-Starterkit ist mit etwa 15€ ein wirklich sehr günstiger Empfänger. Diesen habe ich mir als Bausatz gekauft und zusammengelötet. Hier sieht man ihn angeschlossen an meinem Notebook auf einem kleinen Tisch im Garten:

Wenn man genau hinschaut, sieht man im rechten Bild den Draht, den ich als Antenne zum Gartenhaus gespannt habe.

Dieser Empfänger wurde mit zwei Quarzen geliefert. Der in der Beschreibung erwähnte Quarz mit 6107kHz legt den 6095kHz Träger von RTL Luxemburg auf die Zwischenfrequenz von 12kHz. In Erlangen habe ich mit dem ebenfalls beiliegenden 6097kHz Quarz bessere Ergebnisse erzielt, da dann der 6085kHz Sender des Bayerischen Rundfunks (Ismaning 50kW) gut zu empfangen ist.

Außerdem habe ich noch Versuche mit dem kleinen Empfänger simpleRX für das digitale Campusradio bit eXpress der Erlanger Universität gemacht. Dieser sendet mit 100W auf 15896kHz in Tennenlohe, das ist nur wenige Kilometer von mir entfernt.

Software

Im Internet gibt es z.B. das Software Radio Dream, was man aber selbst kompilieren muss (wohl wegen der Lizenz). Ich habe einiges aus den Quellen des Matlab-basierten Empfängers Diorama der Universität Kaiserslautern gelernt. Allerdings werden auch dort kompilerte Module genutzt, die ich nicht übersetzt habe und so die Software bei mir nie gelaufen ist.

Um die einzelnen Schritte der Signalverarbeitung wirklich zu verstehen, habe ich sie in Matlab selbst implementiert. Im folgenden zeige ich alle Schritte der Signalanalyse.

DRM Systemspezifikation

Die offizielle Systemspezifikation für DRM gibt es beim European Telecommunications Standards Institute ETSI in der Publications Download Area nach einer kostenlosen Registrierung. Dort liefert eine Suche nach "DRM" im Titel die aktuelle Spezifikation "ETSI ES 201 980 V3.1.1 (2009-08)". Auf diese beziehe ich mich bei der Signalanalyse als "Spec".

Aufzeichnung des Zwischenfrequenzsignals

Der Empfänger wird an den Audio-Eingang des PC angeschlossen und mit einer Batterie versorgt. Es hat bei mir einige Versuche benötigt, bis ich auswertbare Signale bekommen habe. Als große Störquelle hat sich das Netzteil des Notebook erwiesen, so dass ich dann alle Aufzeichnungen im Akkubetrieb ohne Netzteil durchgeführt habe. Außerdem bin ich in den Garten gegangen und habe eine Antenne aus einigen Metern Draht aufgespannt. Auch die Eingangskonfiguration der Audio-Hardware muss so eingestellt sein, dass eine gute Aussteuerung ohne Übersteuerung erzielt wird. Bei mir war die Einstellung "Mikrofon" aber ohne "Boost" gut.

Die Aufzeichnung kann direkt in Matlab mit

    Fs=48000;
    T=4;
    u=wavrecord(T*Fs,Fs,1);

durchgeführt werden. Die Abtastrate Fs muss 48kHz betragen, da das gut zu den DRM-Frequenzen passt. Die Aufzeichnungsdauer T sollte einige Sekunden sein, um genügend Frames der Übertragung zu erhalten. Bei zu großer Aufzeichnungsdauer werden die Datenmengen sehr groß und die zeitlichen Änderungen durch Drift und Änderungen in der Atmosphäre erschweren die Ausswertung.

Zum suchen der optimalen Empfangseinstellungen ist das kleine Skript drmrecord.m nützlich, das kontinuierlich 0.4 Sekunden Blöcke aufnimmt und im Zeit- und Frequenzbereich darstellt.

Das Zwischenfrequenzsignal im Zeitbereich

Im Zeitbereich sollte das Signal wie ein Rauschsignal aussehen:

Die Spitzen dürfen nicht abgeschnitten sein. Ansonsten kann man im Zeitbereich nicht viel erkennen. Evtl. lohnt es sich, einen Ausschnitt von 0.5 Sekunden auf periodische Störungen aus dem Netz (Wiederholrate meist 100Hz) zu untersuchen.

Das Zwischenfrequenzsignal im Frequenzbereich

Die Berechnung des Spektrums erfolgt mit

	N=length(u);
	nn=1:N; tt=(nn-1)*dT; ff=(mod(nn-1+N/2,N)-N/2)/T;
	U=fft(u)*2/N;
	semilogy(ff*1e-3,abs(U));

In diesem Spektrum kann man das DRM-Signal sehr gut erkennen:

Typisch ist die sehr gleichmäßig über die Kanalbandbreite verteilte Leistungsdichte, wobei hier die genutze Bandbreite 10 kHz ist (Spec Table 76 Spectrum occupancy = 3). In diesem Band ragen die drei Pilottöne heraus. Sie werden als Frequency Refrence Pilot Cells bezeichnet und liegen bei 750Hz, 2250Hz und 3000Hz bezogen auf die Mittenfrequenz (Spec Abschnitt 8.4.2.1). Diese haben eine gegenüber den Datenträgern um den Faktor sqrt(2) erhöhte Amplitude. Sie ragen im obigen Spektrum viel weiter heraus, da durch die lange Aufzeichnung eine sehr hohe sepktrale Auflösung entsteht. Somit verteilt sich die Energie der Datenträger durch die Modulation auf viele benachbarte Spektrallinien, während die schmalbdanigen Pilottöne auf einer einzigen Linie liegen.

Die Spitze in der Mitte des Nutzbandes ist wohl ein Rest des Trägers, sie ist nicht immer so dominierend und wird daher bei der folgenden Auswertung nicht berücksichtigt. Die Stufe am Rande des Nutzbandes zeigt den Rauschabstand: je höher diese Stufe, desto besser der Rauschabstand.

Die Abtastung führt dazu, dass man das Spektrum symmetrisch gespiegelt sieht. Je nach Lage der Mischerfrequenz bezogen auf die Trägerfrequenz ist das linke oder das rechte Spektrum in der richtigen Lage.

Für die weitere Verarbeitung ist die genaue Lage des Spektrums wichtig. Ich gebe zunächst nach Augenmaß die Frequenz f0 kurz links von dem ersten Pilotton ein. Folgende Befehle suchen dann die Pilottöne und bestimmen die exakte Mittenfrequenz:

	f0=-11.5;
	frange=4000;
	fpilotmin=f0*1000;
	fpilotmax=fpilotmin+frange;
	ipilotrange=find((ff>fpilotmin)&(ffpk*.5)>0);
	fpilots=ff(ipeaks)';
	fcenter=mean(fpilots'-[750 2250 3000]);

Mit der Frqeuenz fcenter kann nun die Umsetzung in das Basisband erfolgen.

Das Basisbandsignal im Frequenzbereich

Die Umsetzung ins Basisband erfolgt durch Multiplikation mit einem drehenden komplexen Zeiger:

	v48=u.*exp(-2*pi*sqrt(-1)*fcenter/48000*nn);

Für die Bandbreite von unter 12kHz genügt eine Abtastrate für das komplexe Basisbandsignal von 12kSample/Sekunde, d.h. es wird nur jeder vierte Wert des ursprünglichen 48kSample/Sekunde Datenstroms benötigt. Vor der Unterabtastung erfolgt noch eine Tiefpassfilterung zur Verbesserung des Rauschabstandes:

	v48f=filter(h,1,v48);
	v=v48f(length(h):4:end);

Hier die verwendete Impulsantwort h des Filters:

So sieht nun das Spektrum des komplexen Basisbandsignals v aus:

Erkennung des Robustness Mode

Der nächste Schritt der Signalanalyse dient der Zerlegung der Zeitfunktion in einzelne Symbole. Die Nutzinformation wird bei OFDM durch Modulation der Real- und Imaginärteile von einzelnen Trägersignalen übertragen. die über die "useful symbol duration Tu" (Spec 4.4.2.2) orthogonal sind. Die Spec spezifiziert in Table 2 vier verschiedene "Robustness Modes" mit unterschiedlichen Tu. Nach Tu wird das Signal über das "guard interval" Tg periodisch fortgesetzt, um robust gegen Signallaufzeiten zu werden. Die gesamte Dauer eines Symbols ist somit Ts=Tu+Tg.

Die Erkennung des verwendeten Robustness Mode und die zeitliche lokalisierung der Symbole erfolgt über die Autokorrelation. Folgende Schleife probiert alle vier Modes durch:

	for mode=1:4
		switch mode
			case 1, Nu=288;Ng=32;Ns=Nu+Ng;
			case 2, Nu=256;Ng=64;Ns=Nu+Ng;
			case 3, Nu=176;Ng=64;Ns=Nu+Ng;
			case 4, Nu=112;Ng=88;Ns=Nu+Ng;
		end;
		vcorr=v(1:(end-Ns)).*conj(v((1+Nu):(end-Ns+Nu)));
		Nvcorrsym=floor(Nv/Ns-1)*Ns;
		vcorravg{mode}=abs(sum(reshape(vcorr(1:Nvcorrsym),Ns,Nvcorrsym/Ns).'));
		vcorrfilt=filter(ones(1,Ng)/Ng,1,[vcorravg{mode}(end-Ng+1:end) vcorravg{mode} vcorravg{mode}(1:Ng)]);
		vsymsync{mode}=vcorrfilt(1+1.5*Ng:end-0.5*Ng).';
		vsymbase{mode}=min(vsymsync{mode});
		[symsyncval(mode) symsyncind(mode)]=max(vsymsync{mode}-vsymbase{mode});
	end;
	[egal mode]=max(symsyncval);
	m=modes(mode);

Das Signal vcorr ist dort, wo eine periodische Fortsetzung mit der Periode Nu vorliegt, positiv reel. In allen anderen Bereichen ist es eine komplexe Zufallszahl. Mit reshape() wird dieses Signal in Stücke der Länge Ns zerlegt und mit sum() über alle diese Stücke summiert (vcorravg). Eine weitere Glättung ergibt sich durch eine gleitende Mittelung entsprechend der Länge des Guard-Intervalls (vcorrfilt). Im folgenden Bild zeigen die vier glatten Kurven die Ergebnisse vcorrfilt für die entsprechend Modes.

Das Signal ist deutlich als Mode B zu erkennen. Die Struktur modes enthält diverse Konstanten aus der Spec, die zum aktuellen Mode gehörenden werden in m kopiert. Der Kreis zeigt die in symsyncind gespeicherte Lage des Maximums, welches die zeitliche Lage der Symbole angibt. Für den erkannten Mode B ist außerdem das ungefilterte vcorravg dargestellt.

Speicherung in Symbolstruktur

Nach Erkennung des Modes und Bestimmung der zeitlichen Lage kann das Signal v in Symbolstruktur umgespeichert werden:

	isym1=symsyncind(mode);
	nsym1=floor((length(v)+1-isym1)/m.Ns);
	w=reshape(v(isym1:isym1-1+nsym1*m.Ns),m.Ns,nsym1);

Jede Spalte der Matrix w ist ein empfangenes Symbol. Die ersten Nu Werte jeder Spalte werden nun verwendet, um für jedes empfangene Symbol die komplexen Amplituden der einzelnen Träger zu berechnen:

	W=fft(w(1:m.Nu,:)).'.*exp(sqrt(-1)*2*pi*0.5*m.Ng/m.Nu*ones(nsym1,1)*(0:m.Nu-1));

Dabei berücksichtigt der komplexe Zeiger die Phasendrehung von einem Symbolanfang zum nächsten entsprechend der mathematischen Beschreibung des Signals in der Spec in Abschnitt 8.1. Der Operator .' entspricht der Transposition, d.h. in der Matrix W entspricht jede Zeile einem Symbol. Die Summe über alle Symbole zeigt, dass die Frequency Pilots wie spezifiziert in allen Symbolen die gleiche Phase haben und sich somit zu großen Spitzen addieren, während die anderen Werte sich im wesentlichen aufheben:

Außerdem erkennt man, dass die Pilottöne bei Carrier Index 16, 48 und 64 liegen, was genau den Angaben in Table 87 der Spec für Robustness Mode B entspricht. Somit ist die Indizierung konsistent mit der Spec. Allerdings bezeichnet die Spec die oberen Indizes mit negativen Nummern, d.h. vom hier angezeigten Index muss Nu abgezogen werden, um den Wert in der Spec zu erhalten.

Erkennung der Frames

Nachdem nun die Symbole in Frequenz und Zeit lokalisiert sind, muss als nächstes der Beginn der Frames gesucht werden. Table 82 der Spec zeigt, dass bei Mode B ein Frame aus 15 Symbolen besteht (m.Nc). Table 89 gibt die Carrier an, auf denen jeweils im ersten Symbol eines Frames die Time Reference übertragen wird (m.trefi). Dies sind Cells mit fest vorgegebener Phase und um den Faktor sqrt(2) erhöhter Amplitude. Sie zeigen sich deutlich in der Korrelation der Cells:

	nf=floor(nsym1/m.Nc)*m.Nc;
	Wtcorr=W(1:end-m.Nc,1+m.trefi).*conj(W(1+m.Nc:end,1+m.trefi));

Daraus lässt sich die Symbol Nummer itref, bei der ein Frame beginnt, bestimmen:

	% nur Frequenzen nehmen, die keine fref sind
	ctrefi=m.trefi~=m.frefi(1) & m.trefi~=m.frefi(2) & m.trefi~=m.frefi(3);
	tcorr=sqrt(sum(Wtcorr(:,ctrefi).*conj(Wtcorr(:,ctrefi)),2));
	Ntcorr=floor(length(tcorr)/m.Nc)*m.Nc;
	[egal itref]=max(sum(reshape(tcorr(1:Ntcorr),m.Nc,Ntcorr/m.Nc)'));
	iitref=mod(itref-1,m.Nc)+1:m.Nc:nsym1;

Das folgende Bild zeigt die Phasenfehler der Frequenzpiloten (farbige Punkte) für alle Symbole und die Phasenfehler der Time Reference Cells für die Anfangssymbole der Frames (schwarze Kreise):

	plot(mod(angle(S(:,1-m.Kmin+m.frefi))-2*pi*ones(nsym1,1)*m.frefph+pi,2*pi)-pi,'.-');
	hold on;
	for i=1:length(m.trefi)
		plot(iitref,mod(angle(S(iitref,1-m.Kmin+m.trefi(i)))-2*pi*ones(length(iitref),1)*m.trefph(i)+pi,2*pi)'-pi,'ko-');
	end;

Hier erkennt man zum ersten Mal, dass die berechneten Winkel der Referenzzellen zu den vorgegebenen Referenzwinkeln passen. Bei Zuordnungs- oder Vorzeichenfehlern liegen die Werte nicht auf glatten Linien, sondern quasi zufällig verteilt. Die verbleibenden, relativ glatten Abweichungen werden im nächsten Schritt korrigiert.

Lokale Amplituden- und Phasenkorrektur mit den Gain Pilots

Auf dem Übertragungskanal gibt es Signalverfälschungen, die sich über Zeit und Frequenz verändern. Um diese korrigieren zu können, sind im gesamten Zeit-Frequenz-Bereich eines Frames Gain Reference Cells mit vorgegebener Amplitude und Phase verteilt (Spec 8.4.4). Die Position der Gain Reference Cells wird entsprechend Table 93 berechnet:

	s=0:m.Nc-1;
	Ms=s'*ones(1,kmax-kmin+1);
	Mn=mod(Ms,m.y);
	Mm=floor(Ms/m.y);
	Mp=(Mk-m.k0-m.x*Mn)/(m.x*m.y);
	Mgainflag=floor(Mp)==Mp;
	[gs, gk]=find(Mgainflag);
	gref=find(Mgainflag);

Die zugehörigen Amplituden sind in Spec 8.4.4.2 und die Phasen in Spec 8.4.4.3 wie folgt spezifiziert:

	Mth=mod(4*Mz256+Mp.*Mw1024+Mp.*Mp.*(1+Ms)*q1024,1024)/1024*2*pi;
	% time_ref anf freq_rfef cell phase definition has priority
	Mth(1,-kmin+1+m.trefi)=m.trefph*2*pi;
	Mth(:,-kmin+1+m.frefi)=ones(m.Nc,1)*m.frefph*2*pi;
	Mamp=[2*ones(m.Nc,1-m.kboost-m.Kmin) sqrt(2)*ones(m.Nc,2*m.kboost-1) 2*ones(m.Nc,1-m.kboost+m.Kmax)];

Für den ersten vollständigen Frame F1 (mit fno=0)

	fno=0
	sym1=itref+m.Nc*fno;
	F1=S(s+sym1,:);

ergibt sich die im folgenden Bild dargestellte Amplituden- und Phasenkorrektur:

	Gdth=(Mamp.*exp(sqrt(-1)*Mth))./F1;

Der Blick "von oben" zeigt die Lage der Pilottöne im Frame:

Durch Interpolation zwischen den Pilottönen ergibt sich die Kanalentzerrung für alle Zellen des Frames:

	Mdth=griddata(Ms(gref),Mk(gref),Gdth(gref),Ms,Mk);

Der entzerrte Frame

Mit obiger Kanalentzerrung wird der entzerrte Frame Z berechnet:

	% Entzerrte Werte
	Z=F1.*Mdth;

Folgendes Bild zeigt die so empfangenen Zellen in der komplexen Ebene:

Die gain reference cells liegen genau auf den Kreisen, da diese ja zur Kanalentzerrung verwendet wurden. Die time references und die frequency references liegen auch etwa auf dem erwarteten Kreis, wobei die time references bewusst unterschiedliche Phasen haben.

Die Nutzdaten des Fast Access Channel (FAC) nutzen laut Spec 8.5.2.2 QAM-4 mit Amplitude 1/sqrt(2) nach Spec figure 39, die roten Punkte passen genau dazu.

Die blauen Punkte zeigen die Nutzdaten, die mit QAM-64 nach Spec p.131 mit einer Quantisierung von 1/sqrt(42) übertragen werden (blaue Quadrate).

Zusammenfassung

Die Demodulation ist damit geschafft, die Konstellationen entsprechen sehr gut dem für QAM-4 und QAM-64 erwarteten Muster.

In einem DRM-Empfänger würden nun aus den Konstellationen die Datenbits gebildet und die Datenströme decodiert. Vielleicht werde ich mir diesen Teil der Signalverarbeitung auch mal anschauen.