Kodowanie sygnału mowy
Autorzy: Przemysław Dymarski
Pliki: rozdzial_12.zip
Przykład 12.1
W przykładzie 12.1 (plik przyklad_12_1.m) przedstawiono kwantyzator równomierny (funkcja kwant_rown), który przeprowadza kwantyzację pojedynczej próbki (we). Prawidłowa praca kwantyzatora odbywa się w zakresie od −xmax do xmax, liczba poziomów wynosi L. Kwantowanie przeprowadza się drogą zaokrąglenia znormalizowanej próbki do najbliższej liczby całkowitej (w programie wykorzystuje się numery poziomów kwantyzacji ze znakiem – ujemna próbka otrzymuje niedodatni numer poziomu). Proces zaokrąglania do najbliższej liczby całkowitej przebiega nieco inaczej dla parzystej liczby poziomów (tzw. kwantyzator typu mid-riser, w którym nie występuje poziom o wartości 0) i dla nieparzystej liczby poziomów (kwantyzator mid-tread, w którym występuje poziom o wartości 0). Do testowania kwantyzatora służy skrypt kwant_test, opisany w Dodatku D12.1.% Przykład_12_1: skrypt kwant_test.m - testowanie kwantyzatora równomiernego
% Dane wejściowe:
%
% Nazwa pliku dźwiękowego (bez rozszerzenia WAV)
%
% Amplituda sygnału wejściowego: wpisując 1, zapewniamy pełne wysterowanie kwantyzatora,
% gdyż przeprowadza on kwantyzację od -1 do +1; wpisując wartości np. 0,1, 0,01,
% 0,001 tłumimy frazę wejściową odpowiednio o 20, 40, 60 dB; wpisując wartości > 1
% można uzyskać efekt przesterowania;
%
% Liczba poziomów kwantyzacji L.
%
% Dane wyjściowe:
%
% Plik SYNT_WAV gotowy do odsłuchu;
%
% Wykres sygnału wejściowego, wyjściowego i błędu kwantyzacji (Figure 1);
%
% snrdb, snrsegdb – wartości SNR globalnego i segmentowego (w dB);
%
% Wykres sygnału wejściowego i wyjściowego (Figure 10);
%
% Energia sygnału w poszczególnych segmentach (w segmencie jest 80 próbek, czyli 10 ms
% sygnału) – wykres niebieski (Figure 20);
%
% SNR w dB w poszczególnych segmentach – wykres czerwony (Figure 20).
clear
close all
% nazwa pliku audio
fichier = input('plik audio ','s');
nom_fichier = [fichier '.wav'];
we=wavread(nom_fichier);
amp = input('amplituda sygnalu: ');
sig=we*amp/max(abs(we)); % tu mozna stlumic lub wzmocnic sygnal
N=length(sig);
L = input('liczba poziomow ');
% kwantyzacja
for i=1:N
[indx(i) qy] = kwant_rown(L, 1, sig(i));
end
%
% dekodowanie
for i=1:N
qsig(i) = dekod_rown(L, 1, indx(i));
end
qsig=qsig';
qerr = sig-qsig;
snr_(sig,qsig);
figure(1), hold off
plot(sig), hold on
plot(qsig, 'g')
plot(qerr, 'r')
title('we, wy i blad kwantyzacji')
wavwrite(qsig,'synt.wav')
function [ind, wy] = kwant_rown(L,xmax,we)
function [ind, wy] = kwant_rown(L,xmax,we)
%
% kwantyzator rownomierny, symetryczny:
% midtread dla nieparzystej liczby poziomów,
% midriser dla parzystej liczby poziomów.
% L: liczba poziomów
% xmax: zakres pracy "w gore" (kwantowanie w zakresie od -xmax do xmax)
% we: probka we
% wy: probka wy
% ind: nr poziomu kwantyzacji: numeracja od 1 do L/2 ze znakiem dla mid-riser
% od 0 do (L-1)/2 ze znakiem dla mid-tread
if L == 0
wy = 0;
ind = 0;
elseif fix(L/2)*2 == L % parzyste L
delta = 2*abs(xmax)/L; % odstep miedzy sasiednimi poziomami kwantyzacji
xx=delta*((L-1)/2); % najwyższy poziom kwantyzacji
if abs(we) >= abs(xmax)
tmp = xx; % probka w stanie przesterowania
ind=L/2;
else
tmp = abs(we); % probka w zakresie pracy kwantyzatora
ind=round((tmp+delta/2)/delta); % tu nastepuje kwantowanie przez zaokrąglenie do najblizszej liczby calkowitej 1,2,...,L/2
end
wy =ind*delta-delta/2;
else % nieparzyste L
delta = 2*abs(xmax)/L; % odstep miedzy sasiednimi poziomami kwantyzacji
xx=delta*((L-1)/2); % najwyższy poziom kwantyzacji
if abs(we) >= abs(xmax)
tmp = xx; % probka w stanie przesterowania
ind=(L-1)/2;
else
tmp = abs(we); % probka w zakresie pracy kwantyzatora
ind=round(tmp/delta); % tu nastepuje kwantowanie przez zaokrąglenie do najblizszej liczby calkowitej 0,1,...,(L-1)/2
end
wy =ind*delta;
end
if we ≶ 0 % odzyskanie znaku próbki
wy = -wy;
ind=-ind;
end
Przykład 12.2
W przykładzie 12.2 (plik przyklad_12_2.m) przedstawiono kwantyzator nierównomierny z kompresją logarytmiczną typu mi (kompresję (wzór 12.9) realizuje funkcja kompr_log a operację odwrotną, tzw. ekspansję, funkcja expand_log). Do testowania kwantyzatora służy skrypt qlog_test, opisany w Dodatku D12.2. Oprócz wymienionych dwóch funkcji, wywołuje on również funkcje kwant_rown (kwantyzator równomierny) i snr_ (obliczenie globalnego i segmentowego stosunku mocy sygnału do błędu kwantyzacji).% Przykład_12_2: skrypt qlog_test.m - testowanie kwantyzatora logarytmicznego typu mi
% Dane wejściowe:
%
% Nazwa pliku dźwiękowego (bez rozszerzenia WAV)
%
% Amplituda sygnału wejściowego: wpisując 1, zapewniamy pełne wysterowanie kwantyzatora,
% gdyż przeprowadza on kwantyzację od -1 do +1; wpisując wartości np. 0,1, 0,01,
% 0,001 tłumimy frazę wejściową odpowiednio o 20, 40, 60 dB; wpisując wartości > 1
% można uzyskać efekt przesterowania;
%
% Liczba poziomów kwantyzacji L.
%
% Wartość parametru mi (wzór 12.9)
%
% Dane wyjściowe:
%
% Plik SYNT_WAV gotowy do odsłuchu;
%
% Wykres sygnału wejściowego, wyjściowego i błędu kwantyzacji (Figure 1);
%
% snrdb, snrsegdb – wartości SNR globalnego i segmentowego (w dB);
%
% Wykres sygnału wejściowego i wyjściowego (Figure 10);
%
% Energia sygnału w poszczególnych segmentach (w segmencie jest 80 próbek, czyli 10 ms
% sygnału) – wykres niebieski (Figure 20);
%
% SNR w dB w poszczególnych segmentach – wykres czerwony (Figure 20).
clear
close all
% nazwa pliku audio
fichier = input('plik audio ','s');
nom_fichier = [fichier '.wav'];
we=wavread(nom_fichier);
amp = input('amplituda sygnalu: ');
sig=we*amp/max(abs(we));
N=length(sig);
L = input('liczba poziomów ');
we_nor=sig;
% kompresja mi
mu=input('wartosc mi: ');
for i=1:N
ylog(i) =kompr_log(we_nor(i),1,mu);
end
% kwantyzacja rownomierna
for i=1:N
[indx(i) qy] = kwant_rown(L, 1, ylog(i));
end
% dekodowanie rownomierne
for i=1:N
qylog(i) = dekod_rown(L, 1, indx(i));
end
% ekspander
for i=1:N
wy_nor(i)= expand_log(qylog(i),1,mu);
end
wy_nor=wy_nor';
err = wy_nor-we_nor; % blad kwantyzacji
figure(1), hold off
plot(we_nor), hold on
plot(wy_nor, 'g')
plot(err, 'r')
title('we, wy i blad kwantyzacji')
snr_(we_nor,wy_nor);
wavwrite(wy_nor,'synt.wav')
function wy = kompr_log(we, vmax, mu)
function wy = kompr_log(we, vmax, mu) % kompresja logarytmiczna typu mi % we – próbka wejściowa % wy – próbka poddana kompresji % mu – parametr krzywej kompresji, w standardzie G.711 równy 255 % vmax – zakres pracy kompresora, przy przetwarzaniu plików .wav vmax=1 % we = we/vmax; wy = vmax*sign(we)*log(1+mu*abs(we))/log(1+mu);function wy = expand_log(we,vmax, mu)
function wy = expand_log(we,vmax, mu) % ekspansja - operacja odwrotna do kompresji logarytmicznej we=we/vmax; wy=sign(we)*(vmax/mu)*((1+mu)^abs(we) -1);
Przykład 12.3
W przykładzie 12.3 (plik przyklad_12_3.m) przedstawiono kwantyzator równomierny z adaptacją momentalną. Do testowania kwantyzatora służy skrypt adpcm_test, opisany w Dodatku D12.3. Jest to ten sam skrypt, który służy do testowania układu ADPCM. Aby testować jedynie kwantyzator adaptacyjny, należy zablokować predyktor. Czyni się to się to, wpisując wartość szybkości adaptacji β = 0. Jako program kwantyzacji równomiernej jest wykorzystywana funkcja kwant_rown. Skrypt wywołuje również funkcję snr_ (obliczenie globalnego i segmentowego stosunku mocy sygnału do błędu kwantyzacji).% Przykład_12_3: skrypt adpcm_test.m - testowanie kwantyzatora adaptacyjnego
% Dane wejściowe:
%
% Nazwa pliku wejściowego (bez rozszerzenia WAV);
%
% Amplituda sygnału wejściowego – jak w kwant_test;
%
% Liczba bitów n: można wpisać 2, 3 lub 4, otrzymując kwantyzator o 4, 8 i 16 poziomach.
% Są to kwantyzatory typu mid-riser – nie mają poziomu o wartości zerowej;
%
% Liczba współczynników predykcji: wpisujemy dowolną liczbę od 1 do kilkunastu,
% predyktor i tak będzie generował zerowy sygnał predykcji;
%
% Szybkość adaptacji: 0.
%
% Dane wyjściowe – jak w programie kwant_test.
clear
close all
% nazwa pliku audio
fichier = input('plik audio ','s');
nom_fichier = [fichier '.wav'];
we=wavread(nom_fichier);
amp = input('amplituda sygnalu ≶ 1): ');
sig=we*amp/max(abs(we));
N=length(sig);
qsig=sig;
bits = input('liczba bitow (2, 3 lub 4) : ');
L=2^bits; % liczba poziomow kwantyzacji
% wspolczynniki adaptacji kwantyzatora:
if L==4
xm=[0.8,1.6];
elseif L==8
xm=[0.9,0.9,1.25,1.75];
elseif L==16
xm=[0.8,0.8,0.8,0.8,1.2,1.6,2.,2.4];
else
disp('liczba bitow = 2, 3 lub 4')
pause
end
zmax=1.;
zmin=0.001*zmax;
z=0.2*zmax; % poczatkowy zakres pracy
delta=1.e-5;
dl=1.e-10;
w=0.99;
mp = input('liczba wsp. predykcji (M > 0) :');
beta = input('szybkosc adaptacji beta = ');
ai = zeros(mp,1); % wspolczynniki predykcji
buf=zeros(mp,1);
for i=1:N
snp=buf'*ai; % predykcja
en=sig(i)-snp; % blad predykcji
[nr,wy]=kwant_rown(L,z,en); % kwantyzator rownomierny
z=z*xm(abs(nr)); %adaptacja kwantyzatora metoda "wstecz"
if z ≶ zmin
z=zmin;
end
if z > zmax
z=zmax;
end
qsig(i)=wy+snp; % syg. wyjsciowy
ai=ai+beta*wy*buf; % adaptacja predyktora
buf=[qsig(i);buf(1:mp-1)];
unstab=norm(ai); % wykrywanie niestabilnosci numerycznej
if unstab > 10^6
'niestabilnosc - zmniejsz beta!'
pause
end
end
qerr = sig-qsig;
snr_(sig(20:N),qsig(20:N));
figure(1), hold off
plot(sig), hold on
plot(qsig, 'g')
plot(qerr, 'r')
title('we, wy i blad kwantyzacji')
wavwrite(qsig,'synt.wav')
Przykład 12.4
W przykładzie 12.4 (plik przyklad_12_4.m) przedstawiono układ ADPCM z predyktorem i kwantyzatorem równomiernym z adaptacją momentalną. Do testowania układu służy skrypt adpcm_test, opisany w Dodatku D12.3. Jest to ten sam skrypt, który służy do testowania kwantyzatora adaptacyjnego. Jako program kwantyzacji równomiernej jest wykorzystywana funkcja kwant_rown. Skrypt wywołuje również funkcję snr_ (obliczenie globalnego i segmentowego stosunku mocy sygnału do błędu kwantyzacji). Aby przeprowadzić symulację układu ADPCM, postępujemy jak w Przykładzie_12_3, z tą różnicą, że wpisujemy wartości szybkości adaptacji beta > 0.Adaptacja przebiega według wzoru (12.38). Za duża szybkość adaptacji powoduje utratę stabilności numerycznej – gwałtowne powiększanie się współczynników predykcji i sygnałów, w efekcie program kończy się błędem wykonania. Optymalna szybkość adaptacji zależy też od amplitudy sygnału: ze wzoru (12.38) wynika, że wprowadzana poprawka maleje dla sygnałów o małej amplitudzie. Można to skorygować, zwiększając wartość beta dla takich sygnałów.% Przykład_12_4: skrypt adpcm_test.m - testowanie układu ADPCM
% Dane wejściowe:
%
% Nazwa pliku wejściowego (bez rozszerzenia WAV);
%
% Amplituda sygnału wejściowego – jak w kwant_test;
%
% Liczba bitów n: można wpisać 2, 3 lub 4, otrzymując kwantyzator o 4, 8 i 16 poziomach.
% Są to kwantyzatory typu mid-riser – nie mają poziomu o wartości zerowej;
%
% Liczba współczynników predykcji: typowe wartości od 1 do kilkunastu
%
% Szybkość adaptacji beta (np. 0.1). Jeśli wystąpi niestabilność numeryczna
% w procesie adaptacji, wartość tę należy zmniejszyć.
%
% Dane wyjściowe – jak w programie kwant_test.
clear
close all
% nazwa pliku audio
fichier = input('plik audio ','s');
nom_fichier = [fichier '.wav'];
we=wavread(nom_fichier);
amp = input('amplituda sygnalu ≶ 1): ');
sig=we*amp/max(abs(we));
N=length(sig);
qsig=sig;
bits = input('liczba bitow (2, 3 lub 4) : ');
L=2^bits; % liczba poziomow kwantyzacji
% wspolczynniki adaptacji kwantyzatora:
if L==4
xm=[0.8,1.6];
elseif L==8
xm=[0.9,0.9,1.25,1.75];
elseif L==16
xm=[0.8,0.8,0.8,0.8,1.2,1.6,2.,2.4];
else
disp('liczba bitow = 2, 3 lub 4')
pause
end
zmax=1.;
zmin=0.001*zmax;
z=0.2*zmax; % poczatkowy zakres pracy
delta=1.e-5;
dl=1.e-10;
w=0.99;
mp = input('liczba wsp. predykcji (M > 0) :');
beta = input('szybkosc adaptacji beta = ');
ai = zeros(mp,1); % wspolczynniki predykcji
buf=zeros(mp,1);
for i=1:N
snp=buf'*ai; % predykcja
en=sig(i)-snp; % blad predykcji
[nr,wy]=kwant_rown(L,z,en); % kwantyzator rownomierny
z=z*xm(abs(nr)); %adaptacja kwantyzatora metoda "wstecz"
if z ≶ zmin
z=zmin;
end
if z > zmax
z=zmax;
end
qsig(i)=wy+snp; % syg. wyjsciowy
ai=ai+beta*wy*buf; % adaptacja predyktora
buf=[qsig(i);buf(1:mp-1)];
unstab=norm(ai); % wykrywanie niestabilnosci numerycznej
if unstab > 10^6
'niestabilnosc - zmniejsz beta!'
pause
end
end
qerr = sig-qsig;
snr_(sig(20:N),qsig(20:N));
figure(1), hold off
plot(sig), hold on
plot(qsig, 'g')
plot(qerr, 'r')
title('we, wy i blad kwantyzacji')
wavwrite(qsig,'synt.wav')
Przykład 12.5
W przykładzie 12.5 (plik przyklad_12_5.m) przedstawiono skrypt celp_test (Dodatek D12.4), służący do uruchomienia symulatora kodera CELP, który jest zaprogramowany jako funkcja codec.m i opisany w rozdziale 12.5 książki. Skrypt wykorzystuje funkcje codec (program podstawowy), period (periodogram, szacujący widmo mocy sygnału), filt_ar (filtr o nieskończonej odpowiedzi impulsowej), filt_ma (filtr o skończonej odpowiedzi impulsowej). Testowany koder posiada uproszczoną strukturę, nie ma słownika adaptacyjnego, który jest wykorzystywany w predykcji długookresowej (p. 12.6.1). Przypomina tym koder Low-Delay CELP (standardG.728). Ma filtry adaptacyjnej preemfazy i deemfazy, umożliwiające kształtowanie widma szumu kwantyzacji (p. 12.6.3). Symulator jest elastyczny – zmieniając wartości parametrów, można go dostosować do różnych szybkości transmisji.% Przyklad_12_5: skrypt celp_test.m - testowanie układu CELP w wersji uproszczonej
% Dane wejściowe:
%
% Nazwa frazy wejściowej (bez rozszerzenia WAV);
%
% N=256 – liczba próbek tworzących ramkę, w której będą obliczane współczynniki
% predykcji (można tę wartość modyfikować);
%
% L_sf – liczba wektorów sygnału w ramce (musi być podzielnikiem N, typowa wartość
% L_sf=4). Mając dane N i L_sf program oblicza wymiar przetwarzanych wektorów
% sygnału: M_sf=N/L_sf;
%
% K – liczba wektorów słownika tworzących model sygnału według wzoru (12.46);
%
% L2 – liczbawektorów w słowniku (z L2wektorów programwybiera K,według algorytmu
% Matching Pursuit – p. 12.6.1);
%
% Rodzaj słownika: w programie wykorzystano słownik sekwencji gaussowskich
% (dic_celp = randn(M_sf, L2)). Można zamienić ten słownik na słownik impulsowy,
% zawierający L2=M_sf wektorów, z których każdy ma jedynie jedną składową
% niezerową (równą 1). Należy wówczas odblokować instrukcje:
% L2=M_sf;
% dic_celp=eye(L2);
%
% (Na podstawie tych danych jest szacowana szybkość transmisji (przepływność binarna).
% Jest to wartość orientacyjna, gdyż wiele parametrów, jak wzmocnienia i współczynniki
% predykcji, nie podlega kwantyzacji).
%
% gamma – współczynnik filtracji percepcyjnej. Podstawiono wartość 0,9 – zmiana na
% wartość 1 spowoduje wyłączenie filtracji percepcyjnej i mechanizmu kształtowania
% widma szumu kwantyzacji;
%
% P_perc – liczba współczynników predykcji wykorzystywanych w filtrze syntezy predykcyjnej
% H(z) (rys. 12.28) i w filtrach preemfazy/deemfazy.
%
% Dane wyjściowe:
%
% Plik FRAZA.WAV gotowy do odsłuchu. Jego nazwa różni się od nazwy frazy oryginalnej
% znakiem „_”;
%
% Wykres SNR [dB] w segmentach (ramkach) dla sygnału mowy i dla sygnału percepcyjnego
% (Figure 7);
%
% SNR_dB – wartość średnia SNR w segmentach dla mowy;
%
% SNR_percep_dB – wartość średnia SNR w segmentach dla sygnału percepcyjnego.
%
% Jeśli visu = 1, to są wyprowadzane co ramkę następujące wykresy:
% – sygnał wejściowy, wyjściowy i błąd kwantyzacji (Figure 1),
% – widma amplitudy [dB] ww. sygnałów (Figure 2),
% – sygnał percepcyjny (mowa po preemfazie), jego model i błąd modelowania (Figure 3),
% – widma amplitudy [dB] ww. sygnałów (Figure 4),
% – sygnał pobudzający filtr syntezy (Figure 5) i jego widmo (Figure 6).
clear
close all
visu = 0; % uaktywnienie wizualizacji
fe = 8000; % czestotliwosc probkowania
N = 256; % ramka
P_perc = 20; % rzad filru percepcyjnego
gamma = 0.9; % wspolczynnik filtracji percepcyjnej (gdy gamma=1, brak filtracji)
L2 = 128; % liczba wektorow w slowniku (mozna ja zmienic)
L_sf = input('liczba wektorów w 256-probkowej ramce (2,4,8 lub 16): ');
K = input('liczba wektorow tworzacych model sygnalu = ');
M_sf=N/L_sf; % wymiar wektora (podramki)
% slownik CELP - moze byc gaussowski
dic_celp = randn(M_sf, L2);
%L2=M_sf;
%dic_celp=eye(L2); % to jest slownik impulsowy
bitrate=40*fe/N+(3*K+ceil(log2(factorial(L2)/(factorial(K)*factorial(L2-K)))))*fe/M_sf
% plik wejsciowy
fichier = input('plik sygnalu mowy ','s');
nom_fichier = [fichier '.wav'];
speech=wavread(nom_fichier); % Fraza mowy
Nbre_ech = length(speech);
Nbre_fen = fix(Nbre_ech/N); % liczba ramek
fprintf('Przetwarzanie %3d ramek pliku %s.wav \n',Nbre_fen, fichier)
speechout=codec(speech, fe, N, K, M_sf, P_perc, gamma, dic_celp, Nbre_fen, visu);
nom_out = [fichier '_.wav'];
wavwrite(speechout,fe,nom_out); % plik wyjsciowy
function speechout=codec(speech, fe, N, K, M_sf, P_perc, gamma, dico_celp, Nbre_fen, visu)
% program uproszczonego kodera CELP
function speechout=codec(speech, fe, N, K, M_sf, P_perc, gamma, dico_celp, Nbre_fen, visu)
% wartosci poczatkowe
speechout=zeros(length(speech),1);
inits=0;
initso=0;
L2 = size(dico_celp,2); % liczba wektorow w slowniku
pn_r = zeros(N,1); % ramka sygnalu percepcyjnego
xn_r = zeros(N,1); % ramka sygnalu mowy
etat_ma_perc_r = zeros(P_perc,1); % pamieci filtrow
etat_ar_perc_r = zeros(P_perc,1);
etat_ar_mod = zeros(P_perc,1);
etat_ma_perc = zeros(P_perc,1);
etat_ar_perc = zeros(P_perc,1);
ringing = zeros(M_sf,1); % "dzwonienie" (zero input response)
visu_xn = zeros(6*N,1); % bufory wizualizacji
visu_pn = zeros(6*N,1);
visu_xn_hat = zeros(6*N,1);
visu_pn_hat = zeros(6*N,1);
visu_entr = zeros(6*N,1);
axe_freq = (0:N/2-1)'/N*fe/1000;
snr_x=zeros(Nbre_fen,1); % snr w segmentach mowy
snr_p=zeros(Nbre_fen,1); % snr w segmentach sygnalu percepcyjnego
% Petla ramek
for no_fen = 1:Nbre_fen
if visu
fprintf('%3d ', no_fen)
end
xn_r = speech(inits+1:inits+N); % odczytanie ramki sygnalu
inits=inits+N;
rk = zeros(P_perc+1, 1); % autokorelacje
for k = 0:P_perc
rk(k+1) = xn_r(1:N-k)'*xn_r(k+1:N)/N;
end
rk(1) = rk(1)*1.001+10^(-10); % zabezpieczenie przed sygnałem zerowym
ai = [1; -inv(toeplitz(rk(1:P_perc)))*rk(2:P_perc+1)]; % wspolczynniki predykcji
ai_gamma = ai.*gamma.^(0:P_perc)'; % jw, z tlumieniem
% preemfaza
[temp, etat_ma_perc_r]= filt_ma(ai, xn_r(1:N), etat_ma_perc_r);
[temp, etat_ar_perc_r]= filt_ar(1, ai_gamma, temp, etat_ar_perc_r);
pn_r(1:N) = temp; % sygnal percepcyjny
% Modelowanie CELP: en => en_hat
en = pn_r(1:N);
dico_f = zeros(size(dico_celp));
energies = zeros(1,L2);
for l = 1:L2 % filtracja slownika CELP
dico_f(:,l) = filter(1, ai_gamma, dico_celp(:,l));
energies(l) = dico_f(:,l)'*dico_f(:,l);
end
en_hat = zeros(N,1);
entry = zeros(N,1); % bufor dla wizualizacji sygnalu pobudzajacego
Nbre_sfen = N/M_sf; % petla "podramek" (wektorow)
for no_sous_fen = 1:Nbre_sfen
entree=zeros(M_sf,1);
n1 = (no_sous_fen-1)*M_sf + 1;
n2 = no_sous_fen*M_sf;
target = en(n1:n2) - ringing; % sygnal docelowy
targ=target;
for iter=1:K % modelowanie w K krokach
temp = ((targ'*dico_f).^2)./energies;
[val_max, index_max] = max(temp); % wybor wektora ze slownika
gain = targ'*dico_f(:,index_max)/energies(index_max); % wzmocnienie
targ=targ-gain*dico_f(:,index_max);
entree = entree+gain*dico_celp(:,index_max);
end
[temp, etat_ar_mod] = filt_ar(1, ai_gamma, entree, etat_ar_mod);
en_hat(n1:n2) = temp; % model sygnalu percepcyjnego
entry(n1:n2) = entree;
ringing = filt_ar(1, ai_gamma, zeros(M_sf,1), etat_ar_mod); % "dzwonienie"
end
pn_hat = en_hat;
%%% pn_hat=pn_r; % test: odblokowanie tej instrukcji prowadzi do bezblednego modelowania
% Filtracja pn_hat => xn_hat (deemfaza)
temp = pn_hat(1:N);
[temp, etat_ma_perc]= filt_ma(ai_gamma, temp, etat_ma_perc);
[temp2, etat_ar_perc]= filt_ar(1, ai, temp, etat_ar_perc);
xn_hat = temp2; % model sygnalu mowy (sygnal wyjsciowy dekodera CELP)
speechout(initso+1:initso+N)=xn_hat; % zapis ramki sygnalu wyjsciowego
initso=initso+N;
% snr
temp = norm(pn_r(1:N))/(norm(pn_r(1:N)-pn_hat)+10^(-10));
rsb_pn = max([20*log10(temp+0.000001),0]);
snr_p(no_fen)= rsb_pn; % snr dla ramki sygnalu percepcyjnego
temp = norm(xn_r(1:N))/(norm(xn_r(1:N)-xn_hat)+10^(-10));
rsb_xn = max([20*log10(temp+0.000001),0]);
snr_x(no_fen)= rsb_xn; % snr dla ramki mowy
fprintf(' snr_percep = %4.1f snr_speech = %4.1f \n', rsb_pn, rsb_xn)
% wizualizacje
if visu
n5 = length(visu_xn);
n1 = n5 - N + 1;
visu_xn(1:n1-1) = visu_xn(N+1:n5);
visu_xn(n1:n5) = xn_r(1:N);
visu_xn_hat(1:n1-1) = visu_xn_hat(N+1:n5);
visu_xn_hat(n1:n5) = xn_hat;
visu_err_mod = visu_xn(1:n5) - visu_xn_hat(1:n5);
Vmax_xn = max(abs(visu_xn))+10^(-10);
Xk = period(visu_xn(n1:n5));
Xk_hat = period(visu_xn_hat(n1:n5));
Ek_mod = period(visu_err_mod(n1:n5));
visu_pn(1:n1-1) = visu_pn(N+1:n5);
visu_pn(n1:n5) = pn_r(1:N);
visu_pn_hat(1:n1-1) = visu_pn_hat(N+1:n5);
visu_pn_hat(n1:n5) = pn_hat;
visu_err_pmod = visu_pn(1:n5) - visu_pn_hat(1:n5);
Vmax_pn = max(abs(visu_pn))+10^(-10);
Pk = period(visu_pn(n1:n5));
Pk_hat = period(visu_pn_hat(n1:n5));
Ek_pmod = period(visu_err_pmod(n1:n5));
visu_entr(1:n1-1) = visu_entr(N+1:n5);
visu_entr(n1:n5) = entry(1:N);
Vmax_en = max(abs(visu_entr))+10^(-10);
Enk = period(visu_entr(n1:n5));
figure(1), hold off
plot(visu_xn), hold on
plot(visu_xn_hat, 'r')
plot(visu_err_mod, 'g')
plot([n1 n1], [-1.1*Vmax_xn 1.1*Vmax_xn], ':')
plot([n1 n1], [-1.1*Vmax_xn 1.1*Vmax_xn])
plot([n5 n5], [-1.1*Vmax_xn 1.1*Vmax_xn], ':')
axis([1 length(visu_xn) -1.1*Vmax_xn 1.1*Vmax_xn])
title('sygnal we, wy, blad')
figure(2), hold off
plot(axe_freq, Xk), hold on
plot(axe_freq, Xk_hat, 'r')
plot(axe_freq, Ek_mod, 'g')
% axis([0 16 0 100])
title('widma s. we, wy, bledu')
figure(3), hold off
plot(visu_pn), hold on
plot(visu_pn_hat, 'r')
plot(visu_err_pmod, 'g')
plot([n1 n1], [-1.1*Vmax_pn 1.1*Vmax_pn], ':')
plot([n1 n1], [-1.1*Vmax_pn 1.1*Vmax_pn])
plot([n5 n5], [-1.1*Vmax_pn 1.1*Vmax_pn], ':')
axis([1 length(visu_pn) -1.1*Vmax_pn 1.1*Vmax_pn])
title('sygnal percepcyjny: we, wy, blad')
figure(4), hold off
plot(axe_freq, Pk), hold on
plot(axe_freq, Pk_hat, 'r')
plot(axe_freq, Ek_pmod, 'g')
title('widma s. percep.: we, wy, bledu')
figure(5), hold off
plot(visu_entr), hold on
plot([n1 n1], [-1.1*Vmax_en 1.1*Vmax_en], ':')
plot([n1 n1], [-1.1*Vmax_en 1.1*Vmax_en])
plot([n5 n5], [-1.1*Vmax_en 1.1*Vmax_en], ':')
axis([1 length(visu_entr) -1.1*Vmax_en 1.1*Vmax_en])
title('sygnal pobudzajacy filtr syntezy')
figure(6), hold off
plot(axe_freq, Enk), hold on
title('widmo sygnalu pobudzajacego')
drawnow
fprintf('pause \n'); pause
end
end
fprintf('\n')
SNR_dB=mean(snr_x)
SNR_percep_dB=mean(snr_p)
figure(7), hold off
plot(snr_x), hold on
plot(snr_p,'r'), hold on
title('SNR(dB) w segmentach: syg. mowy (nieb.) i percepcyjny (czerw.)')