Podstawy transmisji radiowej
Autorzy: Wiesław Ludwin, Jacek Wszołek, Tomasz Zieliński
Pliki: rozdzial_22.zip
Przykład 22.1
W przykładzie 22.1 (plik przyklad_22_1.m) zaprezentowano program w języku MATLAB realizujący pracę kodera modulatora dla modulacji M-PSK oraz 16QAM (tab. 22.1 i 22.2).%% Przykład 22.1 -- Koder modulatora dla M-PSK oraz 16 QAM
A=2; % Amplituda
M=4; % wartościowowść modulacji M-PSK
n=1:M;
I = A*cos((2*n)*pi/M); % wartość składowej synfazowej
Q = A*sin((2*n)*pi/M); % wartość składowej kwadraturowej
fi=0:0.01:2*pi; %
tytul = strcat(num2str(M), '-PSK'); figure(1);
plot(I,Q,'ro',A*cos(fi),A*sin(fi),'b-'); xlabel('I'); ylabel('Q');
title(tytul); grid on;
%% Koder modulatora dla 16QAM (w tablicy)
Iqam(1) = 0.4472; Qqam(1) =0.4472; Iqam(2) = 0.4472; Qqam(2) =1.3416;
Iqam(3) = 1.3416; Qqam(3) =0.4472; Iqam(4) = 1.3416; Qqam(4) =1.3416;
Iqam(5) = 0.4472; Qqam(5) =-0.4472; Iqam(6) = 0.4472; Qqam(6) =-1.3416;
Iqam(7) = 1.3416; Qqam(7) =-0.4472; Iqam(8) = 1.3416; Qqam(8) =-1.3416;
Iqam(9) = -0.4472; Qqam(9) =0.4472; Iqam(10) = -0.4472; Qqam(10) =1.3416;
Iqam(11) = -1.3416; Qqam(11) =0.4472; Iqam(12) = -1.3416; Qqam(12) =1.3416;
Iqam(13) = -0.4472; Qqam(13) =-0.4472; Iqam(14) = -0.4472; Qqam(14) =-1.3416;
Iqam(15) = -1.3416; Qqam(15) =-0.4472; Iqam(16) = -1.3416; Qqam(16) =-1.3416;
Aqam = unique(abs(Iqam + 1i*Qqam)); % wartości amplitud
figure(2); plot2 = plot(Iqam,Qqam,'ro');
xlabel('I'); ylabel('Q'); title(tytul); grid on; title ('16QAM');
hold on
for k=1:length(Aqam)
plot(Aqam(k)*cos(fi),Aqam(k)*sin(fi),'b-')
end
hold off
Przykład 22.2
W przykładzie 22.2(1) (plik przyklad_22_2.m) przedstawiono implementację uproszczonego modelu transmisji radiowej przez kanał AWGN w paśmie radiowym. Program został napisany z wykorzystaniem dokumentacji i przykładów prezentowanych na stronie internetowej www.mathworks.com/.%% Przykład 22.2 - #1. Transmisja radiowa
%% Nadajnik
% Koder modulatora
M = 4; % wartościowość modulacji
SigLen = 4000; % długość sekwencji
Fc = 4e6; % częstotliwość nośna (Hz)
SymRate = 1e6; % liczba symboli na sekundę
nSamps = 16; % Liczba próbek na symbol
Fs = SymRate * nSamps; % częstotliwość próbkowania
EbNo = 20; % stosunek energii przypadającej na jeden bit do widmowej
% gęstości mocy szumu
% Obliczanie SNR (dla funkcji awgn()), 10*log10(log2(M)/nSamps) bo
% nadpróbkowywujemy, 10*log10(2) bo N0/2
SNR = EbNo + 10*log10(log2(M)/nSamps) + 10*log10(2);
msg = randi([0, 1], SigLen, 1); % Źródło danych
% Blok modulacji (pasmo podstawowe). W przykładnie z rozdziału 22.6
% jest on zaimplementowany bez użycia funkcji modulate()
hMod = modem.pskmod('M', M, 'SymbolOrder', 'gray', 'InputType', 'bit');
% Koder modulatora
TXsig = modulate(hMod, msg); % modulacja
%% Przykład 22.2 - #2 Kształtowanie impulsów
% Parametry filtra pierwiastek z podniesionego kosinusa
FltLength = 8; % rząd filtra
RolloffFactor = 0.22; % współczynnik opadania zbocza charakterystyki
% częstotl. (rolloff factor)
filterSpec = fdesign.pulseshaping(nSamps, 'Square root raised cosine', ...
'Nsym,Beta', FltLength, RolloffFactor);
ImpShapeFlt = design(filterSpec); % tworzymy filtr (obiekt)
TXsigUpsampled = upsample(TXsig, nSamps); % nadpróbkowanie sygnału
% z kodera modulacji
TXshaped = filter(ImpShapeFlt, TXsigUpsampled); % filtracja
% Widmo
Hamming256 = spectrum.welch('Hamming', 256);
BasebandSpectrum = figure;
psd(Hamming256, TXshaped, 'Fs', Fs, 'SpectrumType', ...
'twosided','CenterDC', true);
%% Przykład 22.2 - #3 Konwersja do częstotliwości radiowej z wzoru Eulera
t = (0:1/Fs:(SigLen/(sqrt(M)*SymRate))-1/Fs).'; % czas
Carrier = sqrt(2)*exp(1i*2*pi*Fc*t); % nośna
TXup = real(TXshaped.*Carrier); % konwersja
PassbandSpectrum = figure;
psd(Hamming256, TXup, 'Fs', Fs, 'SpectrumType','twosided', ...
'CenterDC', true);
%% Przykład 22.2 - #4 Kanał AWGN
RXup = awgn(TXup, SNR, 'measured');
% Wyświetlanie widma
figure(PassbandSpectrum);
hold on;
psdProbe = psd(Hamming256, RXup, 'Fs', Fs, 'SpectrumType', ...
'twosided', 'CenterDC', true);
plot2 = plot(psdProbe);
set (plot2, 'Color','r');
legend('Sygnał na wejściu kanału', 'Sygnał na wyjściu kanału', ...
'Location', 'southwest');
hold off;
%% Przykład 22.2 - #5
% Odbiornik
% Konwerter w dół
RXdown = RXup.*conj(Carrier);
% Wyświetlanie widma
ReceivePlot = figure;
psd(Hamming256, RXdown, 'Fs', Fs, 'SpectrumType', 'twosided', ...
'CenterDC', true);
%% Przykład 22.2 - #6
% Filtr odbiorczy
ReceiveFlt = copy(ImpShapeFlt);
RXsym = filter(ReceiveFlt, RXdown);
% Wyświetlanie widma
hold on;
psdProbe = psd(Hamming256, RXsym, 'Fs', Fs, 'SpectrumType', 'twosided',...
'CenterDC', true);
plot3 = plot(psdProbe);
set (plot3, 'Color','r');
legend('Sygnał przed filtracją',...
'Sygnał po filtracji', 'Location', 'southwest')
hold off;
%% Przykład 22.2 - #7
% Detekcja
% Odrzucamy pierwsze FltLength symbole (opóźnienie filtra)
FltDelay = FltLength * nSamps;
RXdown = downsample(RXsym(FltDelay+1:end), nSamps); % decymacja
hDemod = modem.pskdemod(hMod); % demodulator (obiekt demod)
RXmsg = demodulate(hDemod, RXdown); % detekcja
liczba_bledow = biterr(msg(1:end-FltLength*log2(M)), RXmsg); % liczba błędów
Przykład 22.3
W przykładzie 22.3 (plik przyklad_22_3.m) przedstawiono implementację w języku MATLAB modelu transmisji sygnału w pasmie podstawowym w kanale AWGN.%% Przykład 22.3 -- Transmisja radiowa w paśmie podstawowym w kanale AWGN
%% Parametry poczatkowe
siglen = 100000; % dlugość sekwencji
EbNomin = 0; EbNomax = 10; % Zakres EbNo
EbNovec = EbNomin:2:EbNomax; % Wektor EbNo
numEbNos = length(EbNovec); % Liczba punktów EbNo
% Modulatory
QPSKmod = modem.pskmod('M', 4, 'SymbolOrder', 'gray'); % Modulator QPSK
QPSKdemod = modem.pskdemod(QPSKmod); % Demodulator
numerrminQPSK = 20; % Minimalna liczba błędów aby skończyć symulacje
berQPSK = zeros(1,numEbNos); % BERy
serQPSK = zeros(1,numEbNos); % SERy
intv4B = cell(1,numEbNos); % przedziały ufności QPSK
QAM16mod = modem.qammod('M', 16, 'SymbolOrder', 'gray'); % Modulator 16QAM
QAM16demod = modem.qamdemod(QAM16mod); % Demodulator 16QAM
numerrmin16QAM = 20; % Minimalna liczba błędów aby skończyć symulacje
ber16QAM = zeros(1,numEbNos); % BERy
ser16QAM = zeros(1,numEbNos); % SERy
intv16B = cell(1,numEbNos); % przedziały ufności 16QAM
% Symulacje puszczamy w pętli, każda iteracja dla danej wartości EbN0.
%% QPSK
for jj = 1:numEbNos
EbNo = EbNovec(jj);
snr = EbNo + 10*log10(log2(4));
ntrialsQPSK = 0; % Liczba iteracji petli
numerrQPSK = 0; % Liczba błędow, dla danej wartoęci EbN0 (bity)
numSymerrQPSK = 0; % Symbole
% Symuluj tak długo, aż nie wygeneruje się określona liczba błędów.
while (numerrQPSK < numerrminQPSK)
msg = randi([0, 3], siglen, 1); % Generacja danych
txsig = modulate(QPSKmod, msg); % Modulacja
rxsig = awgn(txsig, snr, 'measured'); % Dodanie szumu białego (kanał awgn)
decodmsg = demodulate(QPSKdemod, rxsig); % Demodulacja
newerrsQPSK = biterr(msg, decodmsg); % Liczba błędów w tej iteracji: bity,
newSymerrsQPSK = symerr(msg, decodmsg); % symbole
numerrQPSK = numerrQPSK + newerrsQPSK; % całkowita liczba błędów
numSymerrQPSK = numSymerrQPSK + newSymerrsQPSK;
ntrialsQPSK = ntrialsQPSK + 1; % Aktualizacja liczby iteracji pętli.
end
% Obliczanie BER
[berQPSK(jj), intv4B{jj}] = berconfint(numerrQPSK,(ntrialsQPSK * siglen * log2(4)),.95);
[serQPSK(jj), intv4S] = berconfint(newSymerrsQPSK,(ntrialsQPSK * siglen),.95);
disp(['QPSK: EbNo = ' num2str(EbNo) ' dB, ' num2str(numerrQPSK) ...
' błędów, BER = ' num2str(berQPSK(jj)) ', SER = ' num2str(serQPSK(jj))])
end
%% 16QAM
for jj = 1:numEbNos
EbNo = EbNovec(jj);
snr = EbNo + 10*log10(log2(16));
ntrials16QAM = 0;
numerr16QAM = 0;
numSymerr16QAM = 0;
while (numerr16QAM < numerrmin16QAM)
msg = randi([0, 15], siglen, 1);
txsig = modulate(QAM16mod, msg);
rxsig = awgn(txsig, snr, 'measured');
decodmsg = demodulate(QAM16demod, rxsig);
newerrs16QAM = biterr(msg, decodmsg);
newSymerrs16QAM = symerr(msg, decodmsg);
numerr16QAM = numerr16QAM + newerrs16QAM;
numSymerr16QAM = numSymerr16QAM + newSymerrs16QAM;
ntrials16QAM = ntrials16QAM + 1;
end
% Obliczanie BER
[ber16QAM(jj), intv16B{jj}] = berconfint(numerr16QAM,(ntrials16QAM * siglen * log2(16)),.95);
[ser16QAM(jj), intv16S] = berconfint(numSymerr16QAM,(ntrials16QAM * siglen),.95);
disp(['16QAM: EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr16QAM) ...
' błędów, BER = ' num2str(ber16QAM(jj)) ', SER = ' num2str(ser16QAM(jj))])
end
%% Wykres
% QPSK
plotQPSK = semilogy(EbNovec(:), berQPSK);
set(get(gca,'Xlabel'),'FontSize',10,'FontName','times');
set(get(gca,'Ylabel'),'FontSize',10,'FontName','times');
set(get(gca,'Title'),'FontSize',12,'FontName','times');
set(plotQPSK,'LineWidth',1,'LineStyle','-', 'Color','b');
set(plotQPSK, 'Marker','x', 'MarkerSize',8, 'MarkerEdgeColor', 'b');
grid on;
% 16QAM
hold on
plot16QAM = semilogy(EbNovec(:), ber16QAM);
set(plot16QAM,'LineWidth',1,'LineStyle','-.', 'Color','r');
set(plot16QAM, 'Marker','o', 'MarkerSize',8, 'MarkerEdgeColor', 'r');
hold off
% Narysuj przedziały ufności
hold on;
for jj=1:numEbNos
semilogy([EbNovec(jj) EbNovec(jj)],intv4B{jj},'g-+');
semilogy([EbNovec(jj) EbNovec(jj)],intv16B{jj},'g-+');
end
hold off;
xlabel('Eb/N0 [dB]'); ylabel('Prawdopodobieństwo błędu');
legend('QPSK','16QAM', 'Location', 'southwest');
Przykład 22.4
W przykładzie 22.4 (plik przyklad_22_4.m) przedstawiono implementację w języku MATLAB modelu matematycznego systemu radiowego pasma podstawowego wykorzystującego kanał FNFC z zanikiem płaskim.%% Przykład 22.4 Transmisja radiowa. Model kanału FNFC z zanikiem
% płaskim i addytywnym szumem białym
%% Parametry poczatkowe
siglen = 1000; % długość sekwencji
trainlen = 200; % długość ciagu uczacego
M = 4; % wartościowość modulacji
hMod = modem.pskmod('M', M, 'SymbolOrder', 'gray', 'InputType', 'bit');
hDemod = modem.pskdemod(hMod); % Demodulator
EbNomin = 0; EbNomax = 16; % Zakres EbNo
numerrmin = 10000; % Minimalna liczba błędów aby skończyć symulacje
EbNovec = EbNomin:2:EbNomax; % Wektor EbNo
numEbNos = length(EbNovec); % Liczba punktów EbNo
ber = zeros(1,numEbNos); % BERy
intv = cell(1,numEbNos); % przedziały ufności
%% Kanał
Ts = 1e-7; % czas trwania próbki (symbolu)
Fd = 14; % maksymalna częstotliwość Dopplera (Hz)
flat_ch = rayleighchan(Ts,Fd);
flat_ch.StoreHistory = 1;
%% Korektor
eqlms = lineareq(1,lms(0.03));
%% Symulacje puszczamy w petli, każda iteracja dla danej wartości EbN0.
for jj = 1:numEbNos
EbNo = EbNovec(jj);
snr = EbNo + 10*log10(log2(M));
ntrials = 0;
numerr = 0;
% Symuluj tak długo, aż nie wygeneruje się określona liczba błędów.
while (numerr < numerrmin)
msg = randi([0, 1], siglen, 1); % Generacja danych.
txsig = modulate(hMod, msg); % Modulacja.
% Przejście sygnału przez kanał FNFC z zanikiem płaskim i efektem Dopplera
rxsig = awgn(filter(flat_ch,txsig), snr);
% Korekcja
eqsig = equalize(eqlms,rxsig,txsig(1:trainlen));
decodmsg = demodulate(hDemod, eqsig);
newerrs = biterr(msg(1:trainlen), decodmsg(1:trainlen), log2(M));
numerr = numerr + newerrs;
ntrials = ntrials + 1;
end
[ber(jj), intv1] = berconfint(numerr,(ntrials * siglen),.95);
disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
' błędów, BER = ' num2str(ber(jj))])
end
%% Wykres
plotBER = semilogy(EbNovec(:), ber);
set(plotBER,'LineWidth',1,'LineStyle','-');
set(plotBER, 'Marker','x');
set(plotBER,'MarkerSize',12, 'MarkerEdgeColor', 'k')
set(plotBER, 'Color','k');
grid on;
% title('BER vs EbN0')
xlabel('SNR [dB]');
ylabel('Prawdopodobieństwo błędu');
Przykład 22.5
W przykładzie 22.5 (plik przyklad_22_5.m) przedstawiono implementację w języku MATLAB modelu matematycznego systemu radiowego pasma podstawowego wykorzystującego kanał FSFC z zanikiem selektywnym.%% Przykład 22.5 Transmisja radiowa. Model kanału FSFC z zanikiem selektywnym
% i addytywnym szumem białym
%% Parametry poczatkowe
siglen = 1000; % długość sekwencji
trainlen = 200; % długość ciagu uczacego
M = 4; % wartościowość modulacji
hMod = modem.pskmod('M', M, 'SymbolOrder', 'gray', 'InputType', 'bit');
hDemod = modem.pskdemod(hMod);
EbNomin = 3; EbNomax = 18; % Zakres EbNo
numerrmin = 30000; % Minimalna liczba błędów aby skończyć symulacje
EbNovec = EbNomin:3:EbNomax; % Wektor EbNo
numEbNos = length(EbNovec); % Liczba punktow EbNo
ber = zeros(1,numEbNos); % BERy
intv = cell(1,numEbNos); % Komórka na przedzialy ufnosci
Tber = zeros(1,numEbNos); % Teoretyczne BERy
%% Kanal
Ts = 2.6e-7; % czas trwania próbki
Fd = 10; % maksymalna częstotliwość Dopplera (Hz)
sel_ch = stdchan(Ts, Fd, '3gppTUx'); % obiekt kanału
sel_ch.NormalizePathGains = 0;
ch_delay = sel_ch.ChannelFilterDelay;
sel_ch.ResetBeforeFiltering = 0;
%% Korektor o rzędzie większym niż 1
nWeights = 30; % rząd korektora (liczba współczynników wagowych)
eqrls = lineareq(nWeights, rls(0.98,0.1)); % stworzenie obiektu korektora
eqrls.SigConst = hMod.Constellation; % ustawienie konstelacji sygnału
eqrls.ResetBeforeFiltering = 0;
eqrls.RefTap = floor(nWeights/2);
eq_delay = (eqrls.RefTap -1)/eqrls.nSampPerSym;
%% Symulacje puszczamy w petli, kazda iteracja dla danej wartości EbN0.
for jj = 1:numEbNos
EbNo = EbNovec(jj);
Tber(jj) = berfading(EbNo, 'psk', M, 10);
snr = EbNo + 10*log10(log2(M));
ntrials = 0;
numerr = 0;
% Symuluj tak długo, aż nie wygeneruje sie określona liczba błędów.
while (numerr < numerrmin)
msg = randi([0, 1], siglen, 1); % Generacja danych.
txsig = modulate(hMod, msg); % Modulacja.
% Przejście sygnału przez kanał z zanikiem płaskim i efektem Dopplera
rxsig = awgn(filter(sel_ch,txsig), snr, 'measured'); % Add noise.
% kompensacja opóźnienia kanału (selektywnego częstotliwościowo)
msg_trunc = msg(1:end-2*ch_delay); % opuść ostatnie ch_delay symbole wiadomości
txsig_trunc = txsig(1:end-ch_delay); % oraz sygnalu nadanego
rxsig_trunc = rxsig(ch_delay+1:end); % opuść pierwsze ch_delay symbole sygnału odebranego
% Korekcja
eqsig = equalize(eqrls, rxsig,txsig(1:trainlen));
% kompensacja opoznienia korektora (tez jest filtrem pamieciowym)
RXsig = eqsig(eq_delay+1:end); % opuść pierwsze eq_delay symbole odebranego sygnału
TXmsg = msg_trunc(1:end - 2*eq_delay); % opuść ostatnie eq_delay
% symbole nadanego sygnału
msglen = length(TXmsg);
decodmsg = demodulate(hDemod, RXsig); % Demodulacja.
newerrs = biterr(msg(1:trainlen), decodmsg(1:trainlen), log2(M));
numerr = numerr + newerrs;
ntrials = ntrials + 1;
nBits = ntrials * msglen;
end
% Obliczanie liczby błędów oraz BER
[ber(jj), intv1] = berconfint(numerr, nBits,.95);
disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
' błędów, BER = ' num2str(ber(jj))])
end
%% Wykres
plotBER = semilogy(EbNovec(:), ber);
set(plotBER,'LineWidth',1,'LineStyle','-');
set(plotBER, 'Marker','+');
set(plotBER,'MarkerSize',12, 'MarkerEdgeColor', 'k')
set(plotBER, 'Color','k');
grid on;
% title('BER vs EbN0')
xlabel('Eb/N0 [dB]');
ylabel('Prawdopodobieństwo błędu');
Przykład 22.6
Przykład 22.6 (plik przyklad_22_06.m) przedstawia najważniejsze operacje przetwarzania sygnałów, wykonywane podczas transmisji radiowej z wykorzystaniem modulacji QPSK (4QAM) – wysyła się w nim konkretne bity w nadajniku i dekoduje je w odbiorniku. Po drodze w kolejnych krokach: 1) przechodzi się z bitów na liczby (fazory) zespolone; 2) wstawia się zera pomiędzy te liczby i wygładza sygnał otrzymany w ten sposób, używając filtra pierwiastka z podniesionego kosinusa (funkcje firrcos() i upfirdn()); 3) moduluje się sygnał; 4) symuluje się prosty kanał AWGN, tłumiący i opóźniający; 5) demoduluje się sygnał i eliminuje wpływ kanału; 6) ponownie filtruje się sygnał filtrem pierwiastka z podniesionego kosinusa; 7) przeprowadza się synchronizację za pomocą sygnału pilota; 8) dekoduje się bity i porównuje z bitami nadanymi. Nadawany sygnał składa się ramek. Każda z nich zawiera preambułę synchronizacyjną oraz właściwą informację.
% Przykład 22.6
% Autorzy: Sebastian LIS, Jacek MAJERCZYKI (wyrazili zgodę na publikację),
% Tomasz ZIELIŃSKI (tzielin@agh.edu.pl)
clear all; close all;
% Parametry
Nsym = 100; % liczba próbek na symbol
Nmod = 10; % liczna próbek na nośną
fsym = 1; % liczba symboli na sekundę
Kramek = 3; % liczba wysłanych ramek/bloków symboli
Npsf = 1001; % liczba współczynników filtra kształtującego impuls; maks 1601
snr = 160; % stosunek sygnału do szumu
fpr = Nsym*fsym; % częstotliwość próbkowania
wzm = 0.1; % wzmocnienie kanału
faza = pi/3; % przesunięcie fazowe kanału
kanal = 1; % zniekształcenie wprowadzane przez kanał: 0/1 NIE/TAK
korektor = 1; % korektor kanału: 0/1 NIE/TAK
dane = [ 0 1 0 2 3 0 3 1 1 2 1 0 2 2 0 1 3 3 0 2 1 1 0 3 0 2 1 1 3 1 0 2 ];
synch = [ 3 0 0 1 2 1 3 0 3 2 2 1 3 0 0 1 2 1 3] ; % wzorzec synchro typu STS
% KODER: generowanie zespolnonych symboli QPSK
% nr symbolu = bity: 0=11 1=01 2=00 3=10
% odpowadające fazie: pi/4 3pi/4 7pi/4 5pi/4
% odpowiadające complex: 1+j -1+j -1-j 1-j
% koder = (1/sqrt(2))* [ 1+j, -1+j, -1-j, 1-j ];
koder = [ exp( j*pi/4 ) exp( j*3*pi/4 ) exp( j*7*pi/4 ) exp( j*5*pi/4 ) ];
psf = firrcos( Npsf-1, fsym/2, 0.35, fpr,'rollof', 'sqrt'); % filtr kształtujący
psf1 = sqrt(2)*(fpr/fsym)*psf; psf2 = 2*psf; % pulse shaping filter
plot(psf1'); xlabel('n'); title('psf[n]'); grid; pause
% NADAJNIK #########################################################################
ramka = [ synch dane ]; % pojedynczna wysłana ramka symboli
wyslane = [ 0 ramka ramka ramka ]; % wysłana sekwencja numerów stanów
Ksym = length(ramka); % liczba symboli w jednej ramce
Nramka = Ksym*Nsym; % liczba próbek sygnału na jedną ramkę symboli
IQwyslane = koder( wyslane+1 ); % numery symboli --> stany konstelacji
sygnal = upfirdn(IQwyslane,psf1,Nsym,1); % nadpróbkowanie Nsym razy z użyciem psf
Nsygnal = length(sygnal);
% Rysunki
L = 10000; % ile próbek sygnału pokazać
figure % wykres fazorowy: imag=funkcja(real)
plot(real(sygnal(1:L)),imag(sygnal(1:L)),'b-');
title('Q=f(I)'); xlabel('I[n]'); ylabel('Q[n]'); pause
figure % real, imag
subplot(311); plot(real(sygnal(1:L))); title('I[n]');
subplot(312); plot(imag(sygnal(1:L))); title('Q[n]');
subplot(313); plot( real(sygnal(1:L)).^2 + imag(sygnal(1:L)).^2 ); title('abs(IQ)');
pause
% Modulacja - przejście na wyższą częstotliwość
n = 0:Nsygnal-1;
sygnal = real(sygnal).*cos(2*pi/Nmod*n) - imag(sygnal).*sin(2*pi/Nmod*n);
% KANAŁ ############################################################################
if(kanal==1) sygnal = sygnal*wzm*exp(-j*faza); end % tłumienie i opóźnienie
sygnal = awgn( sygnal, snr ); % dodanie szumu
% ODBIORNIK ########################################################################
% Demodulacja - powrót na niską częstotliwość, eliminacja wpływu kanału
sreal = cos(2*pi/Nmod*n) .* sygnal;
simag = -sin(2*pi/Nmod*n) .* sygnal;
sygnal = sreal + j*simag;
% Korektor FEQ
if(kanal==1 & korektor==1) sygnal = sygnal*(1/wzm)*exp(j*faza); end
% Generacja sygnału do synchronizacji
IQsynch = koder( synch+1 ); % numery symboli --> stany konstelacji
sync = upfirdn( IQsynch, psf1, Nsym, 1 ); % nadpróbkowanie Nsym razy z użyciem psf
Nsync = length(sync);
% Splot z filtrem kształtującym
Nstart = (Npsf+1)/2;
sygnal= conv(sygnal,psf2); Nsygnal = length(sygnal);
sync = conv(sync,psf2); Nsync = length(sync);
Nstart = Npsf; % aktualizacja przesunięcia analizy
% Rysunki
figure % real, imag
subplot(211); plot(real(sygnal(1:L))); title('I[n]');
subplot(212); plot(imag(sygnal(1:L))); title('Q[n]'); pause
figure % wyres fazorowy: imag=f(real)
plot(sygnal(1:L),'b-'); title('Q=f(I) demod'); xlabel('I[n]'); ylabel('Q[n]'); pause
% Synchronizacja (przykład)
korel = abs( conv( fliplr(conj(sync)),sygnal) ); % korelacja sygnału i synchro.
[xmax, imax] = max( korel( 1:round(Nramka) ) ); % znalezienie pierwszego maksimum
korel = korel / max( korel ); % normalizacja
iwork = Nstart+round(Nsym);
iback = (Nsync+round(Nsym))-iwork;
nfirst = imax-iback-1; % próbka środka pierwszego symbolu
figure
n = nfirst + (0 : Ksym*Kramek-1)*Nsym;
plot( real(sygnal), imag(sygnal),'b-', real(sygnal(n)), imag(sygnal(n)),'ro');
title('Q=f(I)'); xlabel('I[n]'); ylabel('Q[n]'); pause
% Detekcja
for k = 0 : Kramek-1
k
detekt = [];
num = nfirst + k*Nramka + round((0:Ksym-1)*Nsym); % położenia symboli
I = real(sygnal(num));
Q = imag(sygnal(num));
for m=1:length(num)
if( I(m)>0 & Q(m)>0 ) detekt = [detekt 0]; end
if( I(m)<0 & Q(m)>0 ) detekt = [detekt 1]; end
if( I(m)>0 & Q(m)<0 ) detekt = [detekt 2]; end
if( I(m)<0 & Q(m)<0 ) detekt = [detekt 3]; end
if( (I(m)==0) | (Q(m)==0) ) disp('ZERO'); pause; end
end
ramka, detekt,
if( sum( abs( ramka - detekt ) ) == 0 ) disp('TO SAMO !!!'); pause;
else disp('BLĄD !!!'); pause;
end
end
Przykład 22.7
Przykład 22.7 (plik przyklad_22_07.m) to program w języku MATLAB, w którym są przedstawione podstawy cyfrowego kodowania i dekodowania audycji radia FM wraz z bitami RDS (po ich odpowiedniej interpretacji na wyświetlaczu odbiornika pojawia się wysłany tekst). Program nie jest optymalny i ma wyłącznie charakter demonstracyjny. Został on napisany na podstawie normy IEC 62106. Umożliwia tworzenie sumarycznego sygnału łącznego (kanały lewy i prawy oraz dane cyfrowe) w paśmie podstawowym (BB), który jest następnie konwertowany na wyższe częstotliwości (IF). Program daje możliwość wczytywania rzeczywistych, wąskopasmowych danych IQ z istniejących odbiorników radia SDR (np. danych z programu SDR# nagranych z częstotliwością próbkowania 250 kHz – 250 kilopróbek na sekundę) oraz ich dekodowania łącznie z bitami RDS.
% Przykład 22.7.
% Przykład cyfrowego przetwarzania sygnałów w radiu FM
% Radio FM - koder i dekoder razem
% Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl
clear all; close all;
% ########################################################
% INICJOWANIE WARTOŚCI PARAMETRÓW, ZAPROJEKTOWANIE FILTRÓW
% Program fm_param.m
fs = 250000; % częstotliwość próbkowania sygnału szerokopasmowego FM
fpilot = 19000; % częstotliwość pilota [Hz]
fsymb = fpilot/16; % częstotliwość próbkowania symboli RDS, 1187.5 Hz
fstereo = 2*fpilot; % częstotliwość nośnej sygnału L-R, 38000 Hz
frds = 3*fpilot; % częstotliwość nośnej danych RDS, 57000 Hz
Abw = 25000; % założona szerokość pasma sygnału audio
L = 500; % długość zastosowanych filtrów nierekursywnych
dt = 1/fs; % okres próbkowania
% Filtr pre-fazy i de-emfazy
f1 = 2120; tau1 = 1/(2*pi*f1); w1 = tan(1/(2*Abw*tau1));
b_de = [w1/(1+w1), w1/(1+w1)]; a_de = [1, (w1-1)/(w1+1)];
b_pre = a_de; a_pre = b_de;
% Filtr kształtowania impulsów (PSF - pulse shaping filter)
Tsymb = 1/fsymb; % czas trwania symbolu RDS
Nsymb = fs/fsymb; % liczba próbek na symbol RDS (210.5263158)
Nsymb4 = floor(4*Nsymb); % liczba próbek 4 symboli
df = fs/Nsymb4; f = 0 : df : 2/Tsymb; Nf=length(f); % widmo z normy
H = zeros(1, Nsymb4); H(1:Nf) = cos(pi*f*Tsymb/4); H(end:-1:end-(Nf-2))=H(2:Nf);
hpsf = fftshift(ifft(H)); hpsf=hpsf/max(hpsf); % odp. impulsowa
hpsf2 = firrcos( Nsymb4, fsymb, 1.0, fs,'rollof','sqrt'); % odp. impuls #2
hpsf2 = hpsf2/max(hpsf2); hpsf2 = hpsf2(1:end-1);
n=1:length(hpsf); n2=1:length(hpsf2);
plot(n,hpsf,'ro-',n2,hpsf2,'bx-'); grid; pause
phasePSF = angle( exp(-j*2*pi*fsymb/fs*[0:Nsymb4-1]) * hpsf' ); % przes. fazowe
% Filtr dolnoprzepustowy LowPass do odzyskania sygnału L+R
hLPaudio = fir1(L,(Abw/2)/(fs/2),kaiser(L+1,7));
% Filtr wąskopasmowy BandPass do odzyskania sygnału pilota (wokół 19 kHz)
fcentr = fpilot; df1 = 1000; df2 = 2000;
ff = [ 0 fcentr-df2 fcentr-df1 fcentr+df1 fcentr+df2 fs/2 ]/(fs/2);
fa = [ 0 0.01 1 1 0.01 0 ];
hBP19 = firpm(L,ff,fa);
% Filtr szerokopasmowy BandPass do odzyskania sygnału L-R (wokół 38 kHz)
fcentr = fstereo; df1 = 10000; df2 = 12500;
ff = [ 0 fcentr-df2 fcentr-df1 fcentr+df1 fcentr+df2 fs/2 ]/(fs/2);
fa = [ 0 0.01 1 1 0.01 0 ];
hBP38 = firpm(L,ff,fa);
% Filtr wąskopasmowy BandPass do odzyskania sygnału RDS (wokół 57 kHz)
fcentr = frds; df1 = 2500; df2 = 5000;
ff = [ 0 fcentr-df2 fcentr-df1 fcentr+df1 fcentr+df2 fs/2 ]/(fs/2);
fa = [ 0 0.01 1 1 0.01 0 ];
hBP57 = firpm(L,ff,fa);
% #################################
% KODER SYGNAŁU RADIA FM - przykład
% Program fm_koder_stereo_RDS.m
% Wczytaj dwa monofoniczne sygnały audio
[x1, fx1, Nbits ] = wavread('mowa44100'); x1=x1.'; % plot(x1);
[x2, fx2, Nbits ] = wavread('bach44100'); x2=x2.'; % soundsc(x1,fx1);
% Pre-emfaza, ch-ka płaska do 2.1 kHz, potem narastanie 20 dB na dekadę
x1 = filter(b_pre,a_pre,x1); x2 = filter(b_pre,a_pre,x2);
% Nadpróbkowanie do częstotliwości fs
[N1,D1] = rat(fs/fx1,1e-4); x1 = resample(x1,N1,D1);
[N2,D2] = rat(fs/fx2,1e-4); x2 = resample(x2,N2,D2);
Nx = min( length(x1),length(x2) ); n=1:Nx; x1=x1(n); x2=x2(n);
% x1 = zeros(1,Nx); % celach testowych
% Bity RDS do zakodowania
Nrds = ceil(Nx/(fs/fsymb)); rds = [];
for k=1:ceil(Nrds/10), rds = [ rds 1 1 1 1 0 0 0 1 1 0 ];
end
% rds = round(rand(1,Nrds)); % w celach testowych
% Kodowanie różnicowe bitów RDS
Nrds = length(rds); r = 0;
for k=2:Nrds, r(k)=abs(rds(k)-r(k-1));
end
Nr = ceil(fs/fsymb *length(r));
% Przejście na impulsy bi-fazowe
bip = [-1 1; 1 -1 ]; brds = bip( r+1,: ); brds=brds'; brds=brds(:);
% Nadpróbkowanie - wstawienie zer i wygładzenie RDS filtrem kształtującym
uprds = zeros(1,Nr); t = dt*(0:Nr-1);
clk = abs(diff(sign(sin(2*pi/Tsymb*t))))/2; clk = [clk 0];
indx=1;
for n=1:Nr
if(clk(n)==1) uprds(n)=brds(indx); indx=indx+1;
end
end
prds = conv( uprds, hpsf); prds = prds(1:Nx); % wygładzenie filtrem psf
clear t clk r brds uprds;
% Złożenie całkowitego sygnału radia FM
n=0:Nx-1; alpha = 2*pi*fpilot/fs*n;
x =0.9*(x1+x2)+0.1*cos(alpha)+0.9*(x1-x2).*cos(2*alpha)+0.05*prds.*cos(3*alpha);
A = max(x); x=x/A; clear prds alpha;
% Końcowa modulacja częstotliwościowa
fc = 0; % częstotliwość nośna radia FM - u nas wokół DC
kf = 75000; % głębokość modulacji częstotliwości
x = exp( j*2*pi*( fc/fs*n + kf*cumsum(x)/fs ) ); x=x.';
I = real(x); Q = imag(x);
wavwrite([I Q],fs,32,'FM_Radio_IQ_250KS.wav'); % zapisz wynik na dysku
% ###########################################################
% DEKODER SYGNAŁU RADIA FM (zsyntezowanego lub rzeczywistego)
% Dekoder RDS według standardu IEC 62106
% program fm_dekoder_stereo_RDS.m
% Określenie danych wejściowych
% y = x; clear x;
[y,fs,Nbits] = wavread('FM_Radio_IQ_250KS.wav'); % sygnał zsyntezowany
% [y,fs,Nbits] = wavread('FM_Radio_Krakow_101600kHz_IQ_250KS.wav',25*fs); % SDR#
y = y(:,1) + sqrt(-1)*y(:,2); % IQ --> complex
% ODTWORZENIE SYGNAŁU MONO
% Obliczenie częstotliwości chwilowej - demodulacja częstotliwościowa
dy = (y(2:end).*conj(y(1:end-1)));
y = atan2( imag(dy), real(dy) ); clear dy;
% Filtracja dolnoprzepustowa sygnału L+R (mono)
ym = filter( hLPaudio, 1, y );
% Pozostawienie co fs/Abw-tej próbki
ym = ym(1:fs/Abw:end);
% Odsłuchanie sygnału L+R (mono) i zapisanie go na dysk
soundsc(ym,Abw);
w = ym-mean(ym); w=w/max(abs(w)); wavwrite(w,Abw,'FM_Radio_mono.wav' ); clear w;
% ODTWORZENIE SYGNAŁU STEREO
% Odseparowanie sygnału pilota 19 kHz (filtracja pasmowoprzepustowa)
p = filter(hBP19,1,y);
% Pętla PLL z filtrem typu IIR do odtworzenia częstotliwości i fazy pilota [7]
% i na tej podstawie sygnałów nośnych: symboli c1, stereo c38 i danych RDS c57
freq = 2*pi*fpilot/fs; theta = zeros(1,length(p)+1);
alpha = 1e-2; beta = alpha^2/4;
for n = 1 : length(p)
perr = -p(n)*sin(theta(n));
theta(n+1) = theta(n) + freq + alpha*perr;
freq = freq + beta*perr;
end
c1(:,1) = cos(theta(1:end-1)/16+phasePSF); % nośna 1,1875 kHz
c38(:,1) = cos(2*theta(1:end-1)); % nośna 38 kHz
c57(:,1) = cos(3*theta(1:end-1)); clear p; clear theta; % nośna 57 kHz
% Filtracja pasmowo-przepustowa sygnału L-R wokół 38 kHZ
ys = filter(hBP38,1,y);
% Przesunięcie sygnału L-R w częstotliwości do DC (38 kHz --> 0 kHz)
ys = real(ys .* c38); clear c38;
% Filtracja dolno-przepustowa
ys = filter( hLPaudio, 1, ys );
% Pozostawienie co fs/Abw-tej próbki
ys = ys(1:fs/Abw:end);
% Synchronizacja czasowa sygnałów L+R i L-R (uwzględnienie opóźnienia L-R)
delay = (L/2)/(fs/Abw); ym = ym(1:end-delay); ys=2*ys(1+delay:end); % 2 od modul
clear ymm yss;
% Odtworzenie kanałów L i R
y1 = 0.5*( ym + ys ); y2 = 0.5*( ym - ys ); clear ym ys;
% De-emfaza, ch-ka płaska do 2.1 kHz, potem opadanie 20 dB na dekadę
y1 = filter(b_de,a_de,y1); y2 = filter(b_de,a_de,y2);
% Odsłuchanie sygnału stereo po de-emfazie
soundsc([y1' y2'],Abw); pause
% Zapisanie sygnału stereo radia FM na dysk
maxi = max( max(abs(y1)),max(abs(y2)) );
wavwrite( [ y1/maxi y2/maxi ], Abw, 'FM_Radio_stereo.wav' ); clear y1 y2;
% ZDEKODOWANIE BITÓW SYGNAŁU RDS
% Filtracja pasmowo-przepustowa wokół częstotliwości 57 kHz
y = filter(hBP57,1,y);
% Przesunięcie sygnału RDS w częstotliwości do DC (57 kHz --> 0 kHz)
y = y .* c57; clear c57;
% Filtracja LowPass dopasowana - filtr kształtujący impulsy RDS
y = filter( hpsf, 1, y );
% Synchronizacja symbolowa (timing recovery)
n1 = 5000; Nmany = 50*210;
for n = 1:2*210-1
synchro(n) = sum( ( y(n1+n-1:n1+n-1+Nmany-1).*c1(n1 : n1+Nmany-1) ) );
end
[v, indx] = max(synchro);
y = y(indx:end); c1=c1(1:end-indx+1); % zsynchronizowanie sygnału z "zegarem"
% Przesunięcie sygnału w częstotliwości do DC (1.1875 kHz --> 0 kHz)
y = y .* c1; y = y/max(y);
% Proste całkowanie - suma bieżąca ostatnich fs/Tsymb próbek
h = ones(1,round(fs/fsymb)); y = filter(h,1,y); y = y/max(y);
% Znalezienie ujemnego i dodatniego zbocza "zegara"
N = length(c1); M=5; M2=(M-1)/2; maxi1 = []; maxi2 = [];
for n = M2+1:N-M2
if( (c1(n-2)>0) & (c1(n-1)>0) & (c1(n+1)<0) & (c1(n+2)<0) ) % ujemne zbocze
maxi1 = [ maxi1 n ];
end
if( (c1(n-2)<0) & (c1(n-1)<0) & (c1(n+1)>0) & (c1(n+2)>0) ) % dodatnie zbocze
maxi2 = [ maxi2 n ];
end
end
if( std(y(maxi1)) > std(y(maxi2)) ) maxi=maxi1; else maxi=maxi2; end % DD
% Zmiana poziomów: {-1, +1} --> {0,1}
rdsbits = (-sign( y( maxi ) )+1)/2; clear c1 maxi maxi1 maxi2;
% Dekodowanie różnicowe z uwzględnieniem bi-fazowości sygnału
rdsbits = abs(rdsbits(2:end) - rdsbits(1:end-1)); rdsbits = rdsbits(2:2:end)
% Zapisanie odczytanych bitów na dysk - dla programu rds.m L. Hollevoeta
save FM_Radio_RDS.txt rdsbits /ascii