Generování signálu v Matlabu

Základy zpracování signálu

Jiří Krejsa, 2012

Generování signálu v Matlabu Jak získám signál ? - měření Jak vyrobit signál v Matlabu ? Řada funkcí, základní jsou v klasickém Matlabu, další v SPT (Signal Processing Toolbox)

1. Deterministické signály Klasické harmonické - funkce sin a cos 2 1.5 1 0.5 amplituda

t = 0:0.001:1; amplituda = 2; skutecna_frekvence = 10; uhlova_frekvence = 2*pi*skutecna_frekvence; uhlova_faze = 0; x = amplituda*sin(uhlova_frekvence*t + uhlova_faze); plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

0 -0.5 -1 -1.5 -2

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1

sin / cos jsou v základním matlabu, na všechno další je již nutný SPT - Signal Processing Toolbox


Jiří Krejsa, 2012

Harmonický signál s proměnným kmitočtem - příkaz chirp 2 1.5 1 0.5 amplituda

koncovy_cas = 1; t = 0:0.001:koncovy_cas; amplituda = 2; pocatecni_frekvence = 10; koncova_frekvence = 10; x = amplituda*chirp(t,pocatecni_frekvence,koncovy_cas,koncova_frekvence); plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

0 -0.5

v čem je to jiné ?

-1 -1.5 -2

0.1

0.2

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1

2 1.5 1 0.5 amplituda

koncovy_cas = 1; t = 0:0.001:koncovy_cas; amplituda = 2; pocatecni_frekvence = 1; koncova_frekvence = 10; x = amplituda*chirp(t,pocatecni_frekvence,koncovy_cas,koncova_frekvence); plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

0

0 -0.5 -1 -1.5 -2

0

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1


Jiří Krejsa, 2012

rozšířená verze (metoda, fáze) - fáze nefachá !

další změny kmitočtu: linear, quadratic, logarithmic

1.5 1 0.5 amplituda

koncovy_cas = 1; t = 0:0.001:koncovy_cas; amplituda = 2; pocatecni_frekvence = 1; koncova_frekvence = 10; pocatecni_faze = pi/8; x= amplituda*chirp(t,pocatecni_frekvence,koncovy_cas,koncova_frekvence,'quadratic',pocatecni_faze); plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

2

0 -0.5 -1 -1.5 -2

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1


Jiří Krejsa, 2012

Harmonický oscilátor s frekvencí závislou na vektorové veličině - příkaz vco (voltage controlled oscillator) 0.5

ridici signal

amplituda = 2; Fs = 100; % vzorkovaci frekvence t = 0:1/Fs:2; % cas v rozsahu 0-2 sec Fc = Fs/10; % referencni frekvence ridici_signal = zeros(1,length(t)); % rozsah hodnot je -1 az 1, -1 odpovida nulove frekvenci 0 odpovida Fc, % 1 odpovida 2*Fc prvniinterval = floor(length(t)/3); ridici_signal(1:prvniinterval) = -0.5; ridici_signal(2*prvniinterval:length(t)) = 0.5;

0

-0.5

x = amplituda*vco(ridici_signal,Fc,Fs);

subplot(2,1,2) plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

0.2

0.4

0.6

0.8

1 1.2 cas t [sec]

1.4

1.6

1.8

2

0

0.2

0.4

0.6

0.8

1 1.2 cas t [sec]

1.4

1.6

1.8

2

2 1 amplituda

subplot(2,1,1) plot(t,ridici_signal); xlabel('cas t [sec]'); ylabel('ridici signal');

0

0 -1 -2

Co když bude řídící signál harmonický ? Nebo lineární ?


Jiří Krejsa, 2012

amplituda = 2; Fs = 100; % vzorkovaci frekvence t = 0:1/Fs:2; % cas v rozsahu 0-2 sec Fc = Fs/10; % referencni frekvence

ridici signal

1 0.5 0 -0.5

ridici_signal = sin(5*t); -1



0.2

0.4

0.6

0.8

1 1.2 cas t [sec]

1.4

1.6

1.8

2

0

0.2

0.4

0.6

0.8

1 1.2 cas t [sec]

1.4

1.6

1.8

2

2 1 amplituda


0

0 -1 -2

amplituda = 2; Fs = 100; % vzorkovaci frekvence t = 0:1/Fs:2; % cas v rozsahu 0-2 sec Fc = Fs/10; % referencni frekvence ridici_signal = linspace(-1, 1, length(t));

ridici signal

1 0.5 0 -0.5 -1



0.2

0.4

0.6

0.8

1 1.2 cas t [sec]

1.4

1.6

1.8

2

0

0.2

0.4

0.6

0.8

1 1.2 cas t [sec]

1.4

1.6

1.8

2

2 1 amplituda


0

0 -1 -2


Jiří Krejsa, 2012

Obdélník a pila - příkazy square a sawtooth square - Příkaz umožňuje vytvářet sled impulzů o obdélníkovém průběhu, defaultně je stejně dlouhý v kladné i záporné oblasti, změnit to jde parametrem (v procentech) 1.5

amplituda = 1.2; frekvence_signalu = 10; %Hz t = 0:0.001:0.5; % cas v rozsahu 0-0.5 sec

1

x = amplituda * square(2*pi*frekvence_signalu*t); amplituda

plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

0.5

0

-0.5

-1

-1.5

0.05

0.1

0.15

0.2

0.25 0.3 cas t [sec]

0.35

0.4

0.45

0.5

0

0.05

0.1

0.15

0.2

0.25 0.3 cas t [sec]

0.35

0.4

0.45

0.5

1.5

amplituda = 1.2; frekvence_signalu = 10; %Hz t = 0:0.001:0.5; % cas v rozsahu 0-0.5 sec

1

x = amplituda * square(2*pi*frekvence_signalu*t,25); % 25 procent je delka kladne

0.5 amplituda

plot(t,x) xlabel('cas t [sec]'); ylabel('amplituda');

0

0

-0.5

-1

-1.5


Jiří Krejsa, 2012

sawtooth - Příkaz umožňuje vytvářet sled impulzů o pilovitém průběhu, tvar pily se dá měnit parametrem amplituda = 1.2; frekvence_signalu = 10; %Hz t = 0:0.001:0.75; % cas v rozsahu 0-0.75 sec

subplot(2,2,2) plot(t,x2) xlabel('cas t [sec]'); ylabel('amplituda'); title('parametr 0.2') subplot(2,2,3) plot(t,x3) xlabel('cas t [sec]'); ylabel('amplituda'); title('parametr 0.5') subplot(2,2,4) plot(t,x4) xlabel('cas t [sec]'); ylabel('amplituda'); title('parametr 1')

amplituda

amplituda

-1 0

0.2

2

0.4 0.6 cas t [sec] parametr 0.5

-1 0

0.2

0.4 0.6 cas t [sec] parametr 1

0.8

0

0.2

0.4 0.6 cas t [sec]

0.8

2

1

1

0 -1 -2

0

-2

0.8

amplituda

subplot(2,2,1) plot(t,x1) xlabel('cas t [sec]'); ylabel('amplituda'); title('parametr 0')

1

0

-2

parametr 0.2

2

1

amplituda

parametr = 0; x1 = amplituda * sawtooth(2*pi*frekvence_signalu*t,parametr); parametr = 0.2; x2 = amplituda * sawtooth(2*pi*frekvence_signalu*t,parametr); parametr = 0.5; x3 = amplituda * sawtooth(2*pi*frekvence_signalu*t,parametr); parametr = 1; x4 = amplituda * sawtooth(2*pi*frekvence_signalu*t,parametr);

parametr 0

2

0

0.2

0.4 0.6 cas t [sec]

0.8

0 -1 -2


Jiří Krejsa, 2012

Jediný impuls o různých průbězích vyrábí příkazy rectpuls, tripuls a gauspuls. Impulsy jsou centrované kolem nuly a délka je na intervalu otevřeném zprava (na to pozor!) amplituda = 1.2; t = -4:0.01:4; % cas v rozsahu -5 5 sec x1 = amplituda * rectpuls(t); x2 = amplituda * rectpuls(t,3); x3 = amplituda * tripuls(t,3,-1); x4 = amplituda * tripuls(t,2,1);

subplot(2,2,2) plot(t,x2) xlabel('cas t [sec]'); ylabel('amplituda'); title('rectpuls(t,3)')

1

1

0.8 0.6

0.8 0.6

0.4

0.4

0.2

0.2 -2

0 cas t [sec]

2

tripuls(t,3,-1)

0.5

-2

0 cas t [sec]

-2

2

4

0 cas t [sec]

2

4

2

4

tripuls(t,2,1)

1.5

1

0 -4

0 -4

4

amplituda

subplot(2,2,4) plot(t,x4) xlabel('cas t [sec]'); ylabel('amplituda'); title('tripuls(t,2,1)')

1.2

1.5

amplituda

subplot(2,2,3) plot(t,x3) xlabel('cas t [sec]'); ylabel('amplituda'); title('tripuls(t,3,-1)')

1.2

0 -4

rectpuls(t,3)

1.4

amplituda

amplituda

subplot(2,2,1) plot(t,x1) xlabel('cas t [sec]'); ylabel('amplituda'); title('rectpuls')

rectpuls

1.4

1

0.5

0 -4

-2

0 cas t [sec]


Jiří Krejsa, 2012

gauspuls vyrábí harmonický signál, jehož „tvar“ (amplituda) se mění podle Gaussovky. Má více parametrů: gauspuls(t,fc,bw), kde t je čas, fc je frekvence harmonického signálu a bw je šířka pásma (bw = band width). Šířka pásma je ve zlomku (má rozsah 0 - 1). 1.5

t = -1:0.01:1; % cas v rozsahu -1 1 sec

subplot(1,2,1) plot(t,x1) xlabel('cas t [sec]'); ylabel('amplituda'); title('sirka pasma 0.1') subplot(1,2,2) plot(t,x2) xlabel('cas t [sec]'); ylabel('amplituda'); title('sirka pasma 0.5')

amplituda

x1 = amplituda * gauspuls(t, frekvence, sirka_pasma); sirka_pasma = 0.5; x2 = amplituda * gauspuls(t, frekvence, sirka_pasma);

sirka pasma 0.1

1.5

1

1

0.5

0.5 amplituda

amplituda = 1.2; frekvence = 10; sirka_pasma = 0.1

0

0

-0.5

-0.5

-1

-1

-1.5 -1

-0.5

0 0.5 cas t [sec]

Poznámka: tenhle způsob modulace signálu máte každý u sebe (GSM)

1

sirka pasma 0.5

-1.5 -1

-0.5

0 0.5 cas t [sec]

1


Jiří Krejsa, 2012

Někdy je fajn zjistit, kdy už bude „mrnění“ signálu příliš malé. tc = gauspuls('cutoff',fc,bw,bwr,tpe) kde: o tc je doba (cutoff time), po které obálka opisující generovaný impuls klesne o tpe (ten poslední parametr ve funkci) decibelů vůči maximu obálky. o bwr je referenční úroveň šířky pásma tpe i bwr musí být menší než nula !

% zjistim dobu (v sekundach) tc1 = gauspuls('cutoff',frekvence,sirka_pasma,[],limit)

1

0

-1 -4

-3

-2

-1

1.5

0 cas t [sec] limit -60dB

1

2

3

4 x 10

-5

x 10

-5

1 amplituda

subplot(2,1,1); plot(t1,x1) xlabel('cas t [sec]'); ylabel('amplituda'); title('limit -40dB') subplot(2,1,2); plot(t2,x2) xlabel('cas t [sec]'); ylabel('amplituda'); title('limit -60dB')

0.5

-0.5

% cas si nastavim jen odtud potud t1 = -tc1 : vz_per : tc1; x1 = amplituda * gauspuls(t1,frekvence,sirka_pasma); % jeste jednou totez pro krutejsi limit sirka_pasma = 0.6; % 60 % limit = -90; % az bude signal o -90dB mensi nez spicka tc2 = gauspuls('cutoff',frekvence,sirka_pasma,[],limit) t2 = -tc2 : vz_per : tc2; x2 = amplituda * gauspuls(t2,frekvence,sirka_pasma);

limit -40dB

1.5

amplituda

amplituda = 1.2; frekvence = 50000; % 50kHz vzorkovani = 1000000; % MHz vz_per = 1/vzorkovani; % vzorkovaci perioda sirka_pasma = 0.6; % 60 % limit = -40; % az bude signal o -40dB mensi nez spicka

0.5 0 -0.5 -1 -6

-4

-2

0 cas t [sec]

2

4

6


Jiří Krejsa, 2012

K čemu jsou dobré jednotlivé pulsy ? Dá se z nich sestavit celý signál (pulse train) pomocí funkce pulstran: y = pulstran(t,d,'func'), kde t je vektor časové osy d určuje kolikrát se budou impulsy opakovat ‘func’ označuje typ signálu, konkrétně 'gauspuls' pro generování gaussovsky modulovaného harmonického signálu 'rectpuls' pro obdélníkový signál 'tripuls' pro pilovitý signál

sirka_obdelniku = 10*(1/vzorkovani); % desetinasobek periody vzorkovani y1 = amplituda*pulstran(t,d, 'rectpuls',sirka_obdelniku);

rectpuls

1.5

amplituda

amplituda = 1.2; frekvence = 5; % 50Hz vzorkovani = 1000; % 1 kHz t = 0:1/vzorkovani:1; d = 0:1/frekvence:1;

sirka_pily = (1/frekvence)/3; % tretina frekvence pulsu tvar_pily = -1; y2 = amplituda*pulstran(t,d, 'tripuls',sirka_pily,tvar_pily);

1 0.5 0

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec] tripuls

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1

subplot(2,1,2); plot(t,y2) xlabel('cas t [sec]'); ylabel('amplituda'); title('tripuls')

amplituda

1.5

subplot(2,1,1); plot(t,y1) xlabel('cas t [sec]'); ylabel('amplituda'); title('rectpuls')

1 0.5 0


Jiří Krejsa, 2012

Parametr d může být dvousloupcová matice. V tom případě je první sloupec frekvence impulsů a druhý sloupec je váha (gain). Příklad s váhou o průběhu půlky sinusovky:

d= [d1' d2']

rectpuls

2 1 amplituda

amplituda = 1.2; frekvence = 10; % 50Hz vzorkovani = 1000; % 1 kHz t = 0:1/vzorkovani:1; d1 = 0:1/frekvence:1; % odstupy d2 = sin(2*pi*d1); % gain

sirka_obdelniku = 10*(1/vzorkovani); % desetinasobek periody vzorkovani y1 = amplituda*pulstran(t,d, 'rectpuls',sirka_obdelniku);

0 -1 -2

sirka_pily = (1/frekvence)/3; % tretina frekvence pulsu tvar_pily = -1; y2 = amplituda*pulstran(t,d, 'tripuls',sirka_pily,tvar_pily);

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec] tripuls

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1

2

amplituda

1 0 -1 -2


Jiří Krejsa, 2012

Jako pulsu se dá použít i nějaký vlastní tvar, přičemž se použije interpolace z tohoto tvaru: pulstran(t, d, muj_puls, vzorkovani_meho_pulsu) přičemž muj_puls pokrývá [0,(length(muj_puls)-1)/ vzorkovani_meho_pulsu] a mimo rozsah se bere roven nule. amplituda = 1.2; frekvence = 10; % Hz vzorkovani = 10000; % Hz

0.8 hodnota

muj_puls = hamming(25); % moj vlastni prubeh pulsu t = 0:1/vzorkovani:1; d = 0:1/frekvence:1; % odstupy

muj puls

1

0.6 0.4 0.2

y = amplituda*pulstran(t,d,muj_puls,400);

0

subplot(2,1,1); plot(muj_puls) xlabel('vzorek'); ylabel('hodnota'); title('muj puls')

0

5

10

15

20

25

vzorek sled mych vlastnich pulsu

1.5

amplituda

1

subplot(2,1,2); plot(t,y) xlabel('cas t [sec]'); ylabel('amplituda'); title('sled mych vlastnich pulsu')

0.5 0 -0.5

0

0.1

0.2

0.3

0.4

0.5 0.6 cas t [sec]

0.7

0.8

0.9

1


2. Modulované signály Modulace - změna charakteru vhodného nosného signálu pomocí modulujícího (řídícího)signálu Způsoby modulace: moc moc moc, ukážeme ty základní, je nutný SPT Příkaz modulate Syntaxe: modulate(x, fc, fs, ‘metoda’, argumenty) x - řídící signál fc - frekvence nosné fs - vzorkovací frekvence metoda - string, je jich spousta argumenty - většina metod umožňuje rozsáhlou parametrizaci základní modulace: PPM - pulzně polohová modulace (pulse position modulation) PWM - pulzně šířková modulace (pulse width modulation) AM - amplitudová modulace FM - frekvenční modulace

Jiří Krejsa, 2012


Jiří Krejsa, 2012

PPM - pulzně polohová modulace frekvence1 = 40; % Hz vzorkovani = 1000; % Hz

ridici signal 1

1 0.5

% ridici signal musi byt v rozmezi 0-1, ale pozor na prekryti pulsu! ridici_signal1 = [0 0 0 0 0 0 0 0 0]; ridici_signal2 = [0 0 0 0.5 0.5 0.5 0.9 0.9 0.9];

0 -0.5 -1

1

2

3

4

sirka_pulsu = 0.1;

subplot(4,1,1); stem(ridici_signal1) title('ridici signal 1') subplot(4,1,2); stem(t1,y) xlabel('cas [sec]'); title('modulovany signal 1') subplot(4,1,3); stem(ridici_signal2) title('ridici signal 2') subplot(4,1,4); stem(t2,y2) xlabel('cas [sec]'); title('modulovany signal 2')

6

7

8

9

modulovany signal 1

1

% modulate umi i vratit casovou osu, staci pridat parametr na vystup [y,t1] = modulate(ridici_signal1,frekvence1,vzorkovani,'ppm',sirka_pulsu); [y2,t2] = modulate(ridici_signal2,frekvence1,vzorkovani,'ppm',sirka_pulsu);

5

0.5

0

0

0.05

0.1

0.15

0.2

0.25

cas [sec] ridici signal 2

1

0.5

0

1

2

3

4

5

6

7

8

9

modulovany signal 2

1

0.5

0

0

0.05

0.1

0.15 cas [sec]

0.2

0.25


Jiří Krejsa, 2012

PWM - pulsně šířková modulace frekvence1 = 40; % Hz frekvence2 = 80; vzorkovani = 1000; % Hz ridici_signal1 = [0 0.2 0 0.5 0.1 0.5 0.9 0.7 0.9]; % modulate umi i vratit casovou osu, staci pridat parametr na vystup [y,t1] = modulate(ridici_signal1,frekvence1,vzorkovani,'pwm'); [y2,t2] = modulate(ridici_signal1,frekvence2,vzorkovani,'pwm'); subplot(3,1,1); stem(ridici_signal1) title('ridici signal 1') subplot(3,1,2); stem(t1,y) xlabel('cas [sec]'); title('modulovany signal - frekvence 40 Hz') subplot(3,1,3); stem(t2,y2) xlabel('cas [sec]'); title('modulovany signal - frekvence 80 Hz')

ridici signal 1

1

0.5

0

1

2

3

4

5

6

7

8

9

modulovany signal - frekvence 40 Hz

1

0.5

0

0

0.05

0.1

0.15

0.2

0.25

cas [sec] modulovany signal - frekvence 80 Hz

1

0.5

0

0

0.02

0.04

0.06 cas [sec]

0.08

0.1

0.12


Jiří Krejsa, 2012

Amplitudová modulace

ridici_signal = cos(2*pi*t);

0.5 hodnota

t = 0:1/vzorkovani:2;

-1

subplot(3,1,3); plot(t,y2) xlabel('cas t [sec]'); ylabel('amplituda'); title('modulovany signal (20 Hz)')

0

0.2

0.4

0.6

0

0.2

0.4

0.6

0

0.2

0.4

0.6

1

0.8

1 1.2 1.4 cas [sec] modulovany signal (10 Hz)

1.6

1.8

2

0.8

1.6

1.8

2

1.6

1.8

2

amplituda

0.5 0 -0.5 -1

1

1 1.2 1.4 cas t [sec] modulovany signal (20 Hz)

0.5 amplituda

subplot(3,1,2); plot(t,y) xlabel('cas t [sec]'); ylabel('amplituda'); title('modulovany signal (10 Hz)')

0 -0.5

y = modulate(ridici_signal,frekvence1,vzorkovani,'am'); y2 = modulate(ridici_signal,frekvence2,vzorkovani,'am'); subplot(3,1,1); plot(t,ridici_signal) xlabel('cas [sec]'); ylabel('hodnota'); title('ridici signal')

ridici signal

1

frekvence1 = 10; % Hz frekvence2 = 20; % Hz vzorkovani = 10000; % Hz

0 -0.5 -1

0.8

1 1.2 cas t [sec]

1.4


Jiří Krejsa, 2012

Frekvenční modulace ridici signal

1

t = 0:1/vzorkovani:2;

0.5 hodnota

frekvence1 = 10; % Hz frekvence2 = 20; % Hz vzorkovani = 10000; % Hz

-0.5

ridici_signal = cos(2*pi*t);

-1

y = modulate(ridici_signal,frekvence1,vzorkovani,'fm'); y2 = modulate(ridici_signal,frekvence2,vzorkovani,'fm');

0

0.2

0.4

0.6

0

0.2

0.4

0.6

0

0.2

0.4

0.6

1

0.8

1 1.2 1.4 cas [sec] modulovany signal (10 Hz)

1.6

1.8

2

0.8

1.6

1.8

2

1.6

1.8

2

amplituda

0.5 0 -0.5 -1

1

1 1.2 1.4 cas t [sec] modulovany signal (20 Hz)

0.5 amplituda

poznámka - není problém využít již probraného příkazu vco a dosáhnout téhož výsledku, vyzkoušejte

0

0 -0.5 -1

0.8

1 1.2 cas t [sec]

1.4


Jiří Krejsa, 2012

Jak z toho zase dostanu zpátky původní signál, když mám k dispozici ten modulovaný ? Příkaz demod frekvence1 = 200; % Hz vzorkovani = 10000; % Hz t = 0:1/vzorkovani:1;

y = modulate(ridici_signal,frekvence1,vzorkovani,'am');

pred po

0.5 hodnota

ridici_signal = sin(2*pi*t);

ridici signal pred a po demodulaci

1

0 -0.5

y2 = demod(y,frekvence1,vzorkovani,'am'); -1

subplot(2,1,2); plot(t,y) xlabel('cas t [sec]'); ylabel('amplituda'); title('modulovany signal')

0

0.1

0.2

0.3

0.5 0.6 cas [sec] modulovany signal

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.7

0.8

0.9

1

1

0.4

0.5 amplituda

subplot(2,1,1); plot(t,ridici_signal,'r',t,y2,'b') xlabel('cas [sec]'); ylabel('hodnota'); title('ridici signal pred a po demodulaci') legend('pred','po');

0 -0.5 -1

0.5 0.6 cas t [sec]

Není to nějaké podivné ? Vždyť po demodulaci je to přesně půlka původního signálu ! Jak to ? Demod to dělá blbě!


Jiří Krejsa, 2012

Generování zašumělých signálů Aditivní šum - šum je se signálem sečtený

x [ n ] = s [ n ] + sum [ n ] Jak vyrobit signál, který je zašumělý aditivním šumem o stanoveném poměru SNR (signal to noise ratio) ? 2

⎛ Psignál ⎞ ⎛ Asignál ⎞ ⎛ Asignál ⎞ SNRdB = 10 log10 ⎜ ⎟ = 10 log10 ⎜ ⎟ = 20 log10 ⎜ ⎟ P A A ⎝ šum ⎠ ⎝ šum ⎠ ⎝ šum ⎠

Potřebuji tedy vypočítat efektivní hodnoty signálu a šumu. Když má šum nulovou střední hodnotu, je jeho efektivní hodnota rovna směrodatné odchylce. Tuto směrodatnou odchylku potřebuji určit:

σ=

Asignál 10

SNR 20

1 P = Připomínám, že výkon signálu je suma čtverců podělená počtem vzorků: N

N −1

∑ s [ n] 2

n=0

A efektivní hodnota je druhá odmocnina z výkonu (pro harmonický signál je to amplituda/ 2 ).


Jiří Krejsa, 2012

Příklad - chceme vyrobit sinusovku o amplitudě 2, která bude zatížena aditivním šumem tak, že SNR je 10 a 20 dB. amplituda = 2; vzorkovani = 1000; %Hz frekvence = 2; %Hz % casova osa t = 0:1/vzorkovani:2; % 0 - 2 sec % cisty signal signal = amplituda * sin(2*pi*frekvence*t); % vypocet sumu SNR = 10; %dB A = amplituda / sqrt(2); sigma = A / (10^(SNR/20)); sum = sigma * randn(1,length(t)); %pozor na rand/randn ! SNR2 = 20; %dB sigma2 = A / (10^(SNR2/20)); sum2 = sigma2 * randn(1,length(t)); %celkovy signal signal_sesumem = signal + sum - sigma/2; signal_sesumem2 = signal + sum2 - sigma2/2; subplot(2,1,1) plot(t,signal_sesumem,'b',t,signal,'r') title('signal s aditivnim sumem o SNR = 10 dB'); subplot(2,1,2) plot(t,signal_sesumem2,'b',t,signal,'r') title('signal s aditivnim sumem o SNR = 20 dB');

signal s aditivnim sumem o SNR = 10 dB

3 2 1 0 -1 -2 -3 -4

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.6

1.8

2

signal s aditivnim sumem o SNR = 20 dB

3 2 1 0 -1 -2 -3

0

0.2

0.4

0.6

0.8

1

1.2

1.4


Jiří Krejsa, 2012

Příklad toho, co se stane s nebohým zvukem, když ho zašumíme. Na webu sosněte zvuk „rec.wav“, nebo namluvte vlastní. % nasosnuti zvuku, jeho zasumeni clear all [vzorky,fs,nbits]=wavread('rec.wav'); % zasumime SNR = 10; %dB amplituda = (max(vzorky) - min(vzorky))/2; A = amplituda / sqrt(2); sigma = A / (10^(SNR/20)); sum = sigma * randn(1,length(vzorky)); sum = sum';

0.2 original se sumem

0.15 0.1 0.05 0

vzorky2 = vzorky + sum - sigma/2; % ulozime wavwrite(vzorky2,fs,nbits,'recsum.wav');

-0.05 -0.1

% zobrazime od = 24200; do = od+500; plot(1:do-od+1, vzorky(od:do),'b',1:do-od+1,vzorky2(od:do),'r'); legend('original','se sumem')

-0.15 -0.2

0

100

200

300

400

500

600


Jiří Krejsa, 2012

Šum se přidá spíše při přenosu signálu. Takže nejprve modulujeme, pak zašumíme, demodulujeme a koukáme… % nasosnuti zvuku, jeho modulace, zasumeni a nasledna demodulace clear all [vzorky,fs,nbits]=wavread('rec.wav'); % zmodulujeme AM, nosna bude dvojnasobek vzorkovacky originalnich dat frekvence = fs * 2; vzorkovani = frekvence * 10; radio = modulate(vzorky,frekvence,vzorkovani,'am'); % zasumime modulovany signal SNR = 20; %dB amplituda = (max(radio) - min(radio))/2; A = amplituda / sqrt(2); sigma = A / (10^(SNR/20)); sum = sigma * randn(1,length(radio)); sum = sum'; radio2 = radio + sum - sigma/2; % domodulujeme vzorky2 = demod(radio2,frekvence,vzorkovani,'am'); vzorky2 = vzorky2 * 2; % % ulozime wavwrite(vzorky2,fs,nbits,'recsummod.wav'); % zobrazime od = 24200; do = od+100; plot(1:do-od+1, vzorky(od:do),'b',1:do-od+1,radio(od:do),'r',1:dood+1,radio2(od:do),'k',... 1:do-od+1,vzorky2(od:do),'m'); legend('original','po modulaci','modulece se sumem','demodulovany')

0.08 0.06 0.04 0.02 0 -0.02 -0.04

original po modulaci modulece se sumem demodulovany

-0.06 -0.08

0

20

40

60

80

100

120

Generování signálu v Matlabu

Recommend Documents