% ZZS - 9.cvičení: Spektrální analýza, R.Čmejla, 24.listopadu 2014

%% Příklad 1:    Zobrazení signálů ve frekvenční oblasti
% 
% Zobrazení ve frekvenční oblasti se nazývá spektrem nebo také periodogramem.
% Vztah mezi časovou a frekvenční oblastí je popsán pomocí Fourierovy 
% transformace (FT), při které se signál vyjádří jako lineární dekompozice 
% harmonických průběhů. Protože však pracujeme se vzorkovaným signálem, 
% používáme diskrétní Fourierovu transformaci (DFT). Zrychlení výpočtů nám 
% přinášejí efektivní algoritmy výpočtu DFT, které nazýváme rychlou 
% Fourierovou transformací (FFT). Podmínkou pro použití FFT je, že počet 
% vzorků zpracovávaného signálu N musí být mocninou dvou.

% Funkce fft.m vrací hodnoty dvoustranného komplexního spektra, 
% které musíme podělit počtem vzorků N, abychom získali správné hodnoty.
% 
% Pomocí funkce fft.m zobrazme spektrum jedné periody harmonického signálu 
% o frekvenci 500 Hz a vzorkovací frekvenci 8000 Hz.
% 
close all
figure(1)
clear
fs  = 8000;                       % vzorkovaci frekvence
f0  = 500;                        % zakladni frekvence
T   = 1/f0;                       % perioda
t   = 0:1/fs:T-1/fs;              % casova osa
sig1= 1+cos(2*pi*500*t);
subplot(211), 
stem(t,sig1), xlabel('---> cas [s]'), ylabel('amplituda')
legend('Jedna perioda harmonickeho signalu  f=500 Hz')

SIG1  = fft(sig1);                % vypocet FFT signalu
N     = length(SIG1);             % pocet vzorku signalu  osy
f_osa = fs*t/T;                   % generovani frekvencni osy
                                  % Pro ziskani spravnych hodnot 
                                  % musime spektrum podelit N
subplot(212),
stem(f_osa,abs(SIG1)/N);      % vykresleni amplitudoveho spektra signalu
xlabel('---> frekvence [Hz]'), ylabel('amplituda')
legend('Zrcadlene spektrum harmonickeho signalu v meritku')

%% Příklad 2: Frekvenční osa v periodogramu
% 
% Generujme jednu periodu obdélníkového periodického signálu o frekvenci 
% 500 Hz se stejnosměrnou složkou rovné jedné.
% Zobrazme frekvenční spektrum tohoto průběhu. Na vodorovnou osu periodogramu 
% postupně vynesme:
% a) frekvenci f v Hz na intervalu <0 ... fs),
% b) normovanou frekvenci v radiánech <0 ... 2pi),
% c) normovanou frekvenci <0 ... 2),
% d) a pořadí harmonických k.
% 
figure(2)
clear
fs  = 8000;                       % vzorkovaci frekvence
f0  = 500;                        % zakladni frekvence
T   = 1/f0;                       % perioda
t   = 0:1/fs:T-1/fs;              % casova osa

sig2=1+square(2*pi*500*t);
subplot(511), stem(t,sig2)
legend('Jedna perioda obdelnikoveho signalu')
xlabel('---> cas  [s]')

SIG2   = fft(sig2);               % vypocet FFT signalu
N      = length(SIG2);             % pocet vzorku signalu
f_osa  = fs*t/T;
subplot(512), stem(f_osa,abs(SIG2)/N);    % vykresleni amplitudoveho spektra signalu
legend('Zrcadlene spektrum obdelnika')
xlabel('---> frekvence <0...fs> [Hz]')

subplot(513), stem(2*pi*t/T,abs(SIG2)/N)
xlabel('---> uhlova frekvence <0...2\pi  [rad]>')
legend('b)')

subplot(514), stem(2*t/T,abs(SIG2)/N)
xlabel('---> normovana uhlova frekvence <0...2>  [-]')
legend('c)')

subplot(515), stem(0:N-1,abs(SIG2)/N)
xlabel('---> poradi harmonicke')
legend('d)')

%% Příklad 3: Jednostranné a dvoustranné spektrum
% 
% Velmi často se pro zobrazení signálů ve frekvenční oblasti používají 
% jednostranná spektra. Jejich vztach k dvoustranným jednoznačně vyplývá 
% z teorie Fourierových řad.
% 
% a) Zobrazme dvoustranné amplitudové spektrum se stejnosměrnou složkou 
%    uprostřed
% b) Zobrazme jednostranné amplitudove spektrum
% 
figure(3)
clear
fs  = 8000;                       % vzorkovaci frekvence
f0  = 500;                        % zakladni frekvence
T   = 1/f0;                       % perioda
t   = 0:1/fs:T-1/fs;              % casova osa

sig2=1+square(2*pi*500*t);
subplot(211), stem(t,sig2)
legend('Jedna perioda obdelnikoveho signalu')
xlabel('---> cas  [s]')

SIG2   = fft(sig2);               % vypocet FFT signalu
N      = length(SIG2);             % pocet vzorku signalu
f_osa  = fs*t/T;

%%%%%%%%%%%%%%%%
ShiftSIG2   = fftshift(SIG2);           % vypocet FFT signalu
ind         = floor(N/2)+1;             % index Nyquist.f.
shift_f_osa = fs*([1:N]-ind)/N;         % generovani frekvencni osy
subplot(211)
stem(shift_f_osa,abs(ShiftSIG2)/N); % vykr. ampl.spektra signalu
title('Dvoustranne spektrum obdelnika')
xlabel('---> f [Hz]')
subplot(212)
stem(f_osa(1:ind),[SIG2(1) 2*abs(SIG2(2:ind))]/N);
title('Amplitudove spektrum obdelnika v meritku')
xlabel('---> f [Hz]')

%% Příklad 4: Porovnání oken
% 
% Spektrální rozlišení můžeme ilustrovat na rozlišení dvou blízkých sinusovek.
% Nejmenší možné rozlišení je limitováno šířkou propustného pásma okna
% (šířkou hlavního laloku). Nejlepší rozlišení má pravoúhlé okno. Zlepšení 
% spektrálního rozlišení u ostatních oken se dosáhne zvětšením N (nikoliv 
% však doplněním nul).
% Velikost prvního postranního laloku určuje omezení dynamického rozsahu. 
% Největší omezení má pravoúhlé okno, kdy malý odstup prvního postranního 
% laloku může maskovat signály s malou amplitudou.
% 
% okno               Amin [dB]     šířka propust.pásma
% pravoúhlé             13,3             2pi/N
% von Hannovo           31,5             4pi/N
% Hammingovo            42,5             4pi/N
% Blackmanovo           58,2             6pi/N
% Dolphovo-Čebyševovo   48,6             4pi/N
% Dolphovo-Čebyševovo   76,1             6pi/N
% 	
% Amin odstup prvního postranního laloku [dB]
% N    délka okna
% 
% Vykresleme časové průběhy oken a jejich normovaná amplitudová spektra pro:
% a) pravoúhlé okno
% b) Hannovo okno
% c) Hammingovo okno
% d) Blackmanovo okno
% e) Dolphovo-Čebyševovo okno pro potlačení 48,6 dB
% f) Dolphovo-Čebyševovo okno pro potlačení 76,1 dB
% 
% Porovnejme potlačení z hlediska nepropustného pásma a spektrálního rozlišení
% 
% Řešení:
% a) Pravoúhlé okno má nejmenší potlačení nepropustného pásma 
%    a velmi úzké propustné pásmo. Z toho výplývá, že spektrální 
%    rozlišení je při váhování pravoúhlým oknem velmi dobré.
% b) Potlačení Hannova okna roste. Je to však za cenu dvojnásobné
%    šířky nepropustného pásma (oproti pravoúhlému) a tudíž i 
%    horšího spektrálního rozlišení.
% c) Hammingovo okno má téměř konstantní potlačení v celém 
%    nepropustném pásmu.
% d) Blackmanovo okno má největší potlačení nepropustného pásma, 
%    avšak nehorší spektrální rozlišení.
% e,f) Dolphovo-Čebyševova okna mají rovnoměrnou aproximaci 
%      nepropustného pásma. Oproti tradičním oknům není zde šířka 
%      hlavního laloku určena délkou okna N. 
% 


%% Příklad 5: Potlačení prosakování váhováním signálu
% a) Zobrazme dvakrát opakovanou posloupnost x1[n]=sin(k pi/8) 
%    pro k=0 až 63 a amplitudové spektrum |X1(k)|=FFT{x1[n]}.
% b) Vynásobme posloupnost x1[n] Hammingovým oknem a zobrazme 
%    amplitudové spektrum |X1H(k)| posloupnosti x1H[n].
% c) Zobrazme dvakrát opakovanou posloupnost x2[n]=sin(k pi/10) 
%    pro k=0 až 63 a amplitudové spektrum |X2(k)|=FFT{x2[n]}.
% d) Vynásobme posloupnost x2[n] Hammingovým oknem a zobrazme 
%    amplitudové spektrum |X2H(k)| posloupnosti x2H[n].
% 
% Řešení:
% a) Posloupnost x1[n] bez váhování:
%    V tomto případě je délka signálu celistvým násobkem počtu 
%    period. Zobrazené periodické opakování (dvě posloupnosti 
%    x1[n] za sebou) ukazuje plynulý průběh bez nespojitostí. 
%    K jevu, který se nazývá prosakováním ve spektru zde nedochází. 
%    Zobrazené spektrum |X1(k)| odpovídá harmonickému signálu x1[n].
% 
% b) Posloupnost x1H[n] je x1[n] s váhováním:
%    Posloupnost z předcházejícího případu je násobená (váhovaná) 
%    Hammingovým oknem. To odpovídá konvoluci spektrálních obrazů 
%    harmonického průběhu a Hammingova okna. Do spektra pronikají 
%    dvě falešné spektrální čáry, což je způsobeno šířkou hlavního 
%    laloku spektrálního obrazu Hammingova okna.
% 
% c) Posloupnost x2[n] bez váhování:
%    Posloupnost x2[n] má necelistvý počet period. Na obrázku je 
%    zobrazeno periodické opakování této posloupnosti. Vidíme zde 
%    nespojitost. Protože DFT předpokládá periodické opakování 
%    posloupnosti i mimo analyzované okno, uvedená nespojitost 
%    způsobuje jev, který se nazývá prosakování do spektra (leakage).
% 
% d) Posloupnost x2H[n] je x2[n] s váhováním:
%    Bude-li posloupnost váhovaná, pak nespojitosti na okrajích 
%    budou potlačeny a prosakování bude rovněž potlačeno. Je to 
%    však za cenu horšího spektrálního rozlišení.
figure(5)


%% Příklad 6: Spektrogramy – rozlišení 
% Příklady, které jsme doposud vyzkoušeli, byly jednoduché, neboť 
% frekvenční obsah signálu se neměnil s časem. Většina reálných 
% signálů však bývá komplikovanější. K zobrazení frekvenční informace 
% signálu s časově měnícím se frekvenčním obsahem můžeme použít 
% spektrogram. Při tomto druhu kreslení je na horizontální ose čas 
% a na vertikální ose je frekvence. Barvy popisují úroveň signálu 
% pro všechny frekvenční složky a pro všechny časové okamžiky.
% 
% V této úloze se seznámíme s principem neurčitosti při zpracování 
% číslicových signálů. Uvidíme, že nastavení parametrů spektrální 
% analýzy nám dovolí buď přesné nebo frekvenční rozlišení. V praxi 
% často volíme kompromis.
% 
% Gerujme harmonický signál o frekvenci 1920 Hz a časovém trvání 0,6 s 
% a vzorkovací frekvenci 800 Hz. K němu přičteme signál o frekvenci 
% 2000 Hz, trvající 0,2 s a začínající v čase 0.2 s. Sledujme vztahy 
% mezi časovým a frekvenčním rozlišením. 


%% Příklad 7:    Spektrogramy řečových signálů
% Pro lepší rozlišení změn v čase se při analýze řeči používá kratšího 
% okna a zobrazuje se tzv. širokopásmový spektrogram. Chceme-li dosáhnout 
% lepšího frekvenčního rozlišení, pak použijeme delšího okna a zobrazíme 
% úzkopásmový spektrogram. „Kratší“ a „delší“ volíme vzhledem k délce 
% pitch periody.
% 
% a) Zobrazme širokopásmový spektrogram řečového signálu
% b) Zobrazme úzkopásmový spektrogram řečového signálu
% 
% Řešení:
% a) V širokopásmovém spektrogramu vidíme hrubou informaci o rozložení 
%    frekvencí (rozložení hlavních formantů) a přesnou informaci o časovém 
%    vývoji (přesné určení hranic fonémů).
% 
% b) V úzkopásmovém spektrogramu vidíme velmi přesnou informaci o 
%    spektrálním složení (včetně základní frekvence a jejich násobků) 
%    za cenu nepřesné časové informace. 
% 

%% Příklad 8: 3d zobrazení
% 
% Jiným způsob zobrazení časové a frekvenční závislosti je 3 dimenzionální 
% kreslení, kdy na souřadnici x vynášíme čas, souřadnici y frekvenci a na 
% souřadnici z úroveň signálu.
% 
figure(8)
clear

%% Příklad 9: Fázový vokodér
% Pomocí fázového vokodéru:
% a) Proveďte přímou resyntézu písně why11.wav (pouze analýza - syntéza)
%    Oba signály poslechněte a zobrazte v časové a frekvenční oblasti
% 
% b) Prodlužujte signál o 25 až 250 % 
% c) Zkraťte signál na jednu polovinu
% d) Posuňte skladbu frekvenčně výš
% e) Posuňte skladbu frekvenčně níž
% f) Realizujte robot
% g) Realizujte whisper
% 
% 
clear
figure(9)

[signal,fs]   = wavread('why11.wav');
N             = length(sig);              
window_length = 
window_shift  = 
R             =  

% Analýza
     X = []; k=1;
     for start = 0:window_shift:N-window_length,     
        frame = sig(start+1:start+window_length)...
	        .*hamming(window_length);
        X(:,k)=fft(frame);     
        k=k+1;
     end

% Transformace
     Y=X

% Syntéza
     k=1; N_new=size(Y,2)*window_shift+window_length;
     signal_y=zeros(N_new,1);   
     for start = 0:window_shift:(size(Y,2)-1)*window_shift,     
         segment=real(ifft(Y(:,k),window_length))...
                 .*hamming(window_length);
         signal_y((start+1):(start+window_length))...            
                =signal_y((start+1):(start+window_length))+segment;
         k=k+1;
     end