Ocena jakości sygnałów fonicznych
Autorzy: Jakub Rachwalski, Maciej Bartkowiak
Pliki: rozdzial_16.zip
Przykład 16.1
Przykład (plik przyklad_16_1.m) ilustruje implementację funkcji rozpraszania sygnałów subpasmowych w dziedzinie częstotliwości. Argumentem wejściowym funkcji są sygnały subpasmowe x_re i x_im w postaci macierzy: każdy wiersz reprezentuje kolejną ramkę, a każda kolumna reprezentuje podpasmo. Dodatkowo, funkcja korzysta z wektora częstotliwości środkowych filtrów fc[i], i=1..40.% przykład 16.1:
function [y_re,y_im] = Sprdng(fc,x_re,x_im)
% Rozpraszanie częstotliwościowe dla modelu A
%
% [x_re,x_im] = Sprdng(x_re,x_im)
% f0 - poziomy wektor częstotliwości środkowych podpasm
% x_re, x_im - próbki podpasm w kolejnych ramkach (kolumnami)
N = size(x_re,1); % Liczba ramek
dist = 0.1.^((z(fc(40))-z(fc(1)))/(39*20)); % Odległość w Barkach pomiędzy pasmami
a = exp(-32/(48000*0.1)); % Współczynnik filtru wygładzającego
b = 1-a;
Ex = x_re.^2 + x_im.^2; % Energia chwilowa
% Górne zbocze funkcji poszerzającej
Su = max(4,(24+230./repmat(f0,N,1)-2*log10(Ex)));
Cu = filter(a,[1, -b],dist.^Su); % Wygładzenie w czasie
% Poszerzanie w kierunku wyższych częstotliwości
for n = 1:N % Operacja osobna dla każdej ramki
A_re(n,:) = conv(Cu(n,:),x_re(n,:));
A_im(n,:) = conv(Cu(n,:),x_im(n,:));
end
% Usunięcie zbędnej części splotu i odwrócenie
A_re = fliplr(A_re(:,1:40));
A_im = fliplr(A_im(:,1:40));
% Poszerzanie w kierunku niższych częstotliwości
Cd = dist^31; % Dolne zbocze funkcji poszerzającej
for n = 1:N % Operacja osobna dla każdej ramki
B_re(n,:) = conv(Cd(n,:),A_re(n,:));
B_im(n,:) = conv(Cd(n,:),A_im(n,:));
end
% Ponowne usunięcie zbędnej części splotu i odwrócenie
y_re = fliplr(B_re(:,1:40));
y_im = fliplr(B_im(:,1:40));
end
function bark = z(f)
% Konwersja częstotliwości na Barki
bark = 7*asin(f/650);
Przykład 16.2
Przykład (plik przyklad_16_2.m) pokazuje, jak przeprowadzić kalibrację sygnału, tak, aby średnia moc sygnału była równa wartości wcześniej zadanej (oczekiwnanej). Zmienna ‘dane’ zawiera próbki sygnału, a stała ‘OCZEKIWANA_MOC’ zawiera oczekiwaną wartość mocy.% Przykład 16.2
moc = sum(dane.^ 2) / length(dane); % obliczenie średniej mocy badanego sygnału
faktor= sqrt(OCZEKIWANA_MOC/ moc); % obliczenie stosunku mocy oczekiwanej do średniej
% mocy sygnału
dane = dane * faktor; % przemnożenie badanego sygnału przez otrzymany
% stosunek
Przykład 16.3
Przykład (plik przyklad_16_3.m) pokazuje, jak w algorytmie PESQ przeprowadzane jest filtrowanie IRS badanego sygnału. Zmienna ‘dane’ zawiera próbki sygnału, a zmienna ‘filtr’ zawiera pary: (częstotliwość, wartości tłumienia filtru).% Przykład 16.3
n = length (dane); %obliczenie liczby próbek
pow_of_2= 2^ (ceil( log2( n))); %obliczenie liczby próbek potrzebnych do ...
%... przeprowadzenia operacji filtrowania
rozdzielczosc= Fs/ pow_of_2; % obliczenie jaka, dla takiej liczby próbek...
%...będzie rozdzielczość częstotliwościowa
x= zeros( 1, pow_of_2); %uzupełnienie próbek zerami na końcu...
x( 1: n)= dane; %...do obliczonej wcześniej długości
dane_fft = fft( x, pow_of_2); % przejście na dziedzinę częstotliwościową
faktorDb( 1: pow_of_2/2+ 1)= interp1( filtr( :, 1), ... %zwiększenie rozdzielczości zdefiniowanego..
filtr( :, 2), (0: pow_of_2/2)* rozdzielczosc); %...filtra do obliczonej ‘rozdzielczosc’
faktor= 10.^ (faktorDb/ 20); %zmiana skali filtra na liniową
faktor= [faktor, fliplr( faktor( 2: pow_of_2/2))]; %doklejenie lustrzanego odbicia filtra
dane_fft = dane_fft.* faktor; % przemnożenie spektrum sygnału przez filtr
dane = ifft(dane_fft, pow_of_2); % powrót do dziedziny czasowej
dane_y = dane_y ( 1: n); %skrocenie przefiltrowanego sygnalu do
% pierwotnej dlugosci
Przykład 16.4
Przykład (plik przyklad_16_4.m) pokazuje, jak przeprowadzić dopasowanie czasowe dwóch sygnałów (dane_1 , dane_2) za pomocą korelacji. W algorytmie PESQ jest ono przeprowadzone najpierw na całych sygnałach, a później, ta sama operacja jest wykonywana na poszczególnych fragmentach sygnałów. Zmienna ‘dane_1’ zawiera próbki sygnału odniesienia, a zmienna ‘dane_2’ zawiera próbki ocenianego sygnału.% Przykład 16.4
sygnal_1= fliplr(dane_1); % Odwrócenie sygnal_1 w czasie
Y = conv( dane_2, dane_1); % Obliczenie splotu pomiędzy sygnal_1 (odwróconym), a
% syngal_2
[max_Y, przesuniecie ] = max(Y); % wartość przesunięcia to punkt, dla którego splot przyjmuje
% wartość największą
Przykład 16.5
Przykład (plik przyklad_16_5.m) pokazuje, jak powstaje macierz gęstości mocy poprzez zamianę skali częstotliwościowej z Hz na Bark. Zmienna ‘dane_fft’ zawiera sygnał w dziedzinie częstotliwości.% Przykład 16.5
for ramka = 1: ostatnia_ramka % kolejne kroki są powtarzane dla każdej ramki z
% osobna
for bark = 1: liczba_pasm % kolejne kroki są powtarzane dla każdego
% pasma krytycznego
suma = 0; % zmienna ‘suma’ przechowuje moc
% poszczególnych pasm częstotliwościowych
f=PRAZKI_PER_BARK(bark,:); % zmienna f przechowuje informacje które
% prążki należą do obecnie rozważanego pasma
for k = f(1) : f(2) % kolejne kroki są powtarzane dla każdego
% pasma obecnie rozważanej wstęgi z osobna
suma = suma + ...
dane_fft( ramka, k )^2; % sumujemy moce sygnału ‘dane_fft’ należące
% do rozważanego pasma
end
gestosc_mocy (ramka, bark) = suma; % gęstość mocy dla rozważanego pasma
end
end
Przykład 16.6
Przykład (plik przyklad_16_6.m) pokazuje, jak obliczana jest macierz gęstości głośności.% Przykład 16.6
for ramka = 1: ostatnia_ramka % kolejne kroki są powtarzane dla
% każdej ramki z osobna
for z = 1: liczba_pasm % kolejne kroki są powtarzane dla
% każdego pasma krytycznego
h = Zwicker ( gestosc_mocy (ramka, z) ); % oblicz głośność z mocy zgodnie z
% prawem Zwicker’a ...
gestosc_glosnosci(ramka, z) = h; % ...i zapisz do macierzy ‘gestosc_gloscnosci’
end
end
Przykład 16.7
Przykład (plik przyklad_16_7.m) pokazuje, jak obliczana jest gęstość zaburzeń, biorąc pod uwagę efekt asymetrii.% Przykład 16.7
for ramka = 1: ostatnia_ramka % kolejne kroki są powtarzane dla
% każdej ramki z osobna
for z = 1: liczba_pasm % kolejne kroki są powtarzane dla
% każdego pasma z osobna
ratio = (gestosc_mocy_oceniany (ramka, z) + 50)... % oblicz stosunek gęstości mocy...
/ (gestosc_mocy_odniesienia (ramka, z) + 50); % sygnału ocenianego i odniesienia
h = ratio^ 1.2; % podnieś stosunek do potęgi 1.2
if (h > 12) % stosunek nie może być większy...
h = 12; % niż 12
elseif (h < 3) % jeśli stosunek ma wartość...
h = 0.0; %.mniejszą niż 3, wyzeruj go
end
gestosc_zaburzen_asym(ramka, bark) = ... % przemnóż gęstość zaburzeń przez...
gestosc_zaburzen(ramka, bark) * h; %... obliczony stosunek mocy
end
end
Przykład 16.8
Przykład (plik przyklad_16_8.m) pokazuje, jak obliczany jest współczynnik zaburzeń (tak normalnych jak i asymetrycznych). Sygnał rozważany jest na poziomie fonemów(ramek) i sylab, zgodnie z Rys 14.% Przykład 16.8
wspolczynnik_zaburzen= 0;
for sylaba = 1 : liczba_sylab % rozważnie kolejnych sylab
zaburzenie = 0;
for fonem = 1 : 20 % rozważnie kolejnych fonemów
h = gestosc_zaburzen(fonem); % gęstość zaburzeń w rozważanym fonemie
zaburzenie = zaburzenie + (h^ 6); % sumowanie zaburzeń w fonemach
end
zaburzenie = (zaburzenie / 20)^ (1/6) ;
wspolczynnik= wspolczynnik+ zaburzenie^ 2; % sumowanie zaburzeń kolejnych sylab
end
wspolczynnik = wspolczynnik/ liczba_sylab;
wspolczynnik= wspolczynnik^ (1/ 2);