Компресия на звук

From Ilianko

Звукът най-често е в следствие на движение на тяло във някаква среда (въздух,...). Движението предизвиква промяна на налягането, което се разпространява както вълна във водата. Тъпънчето на ухото преобразува промяната на налягането в сигнал, който мозъка ни възприема като звук.

Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 44000 пъти в sec) се вземат отчети (samples) за стойността на електричския сигнал. Всяко измерване се съхранява като число с фиксирана точност (примерно 8, 16 бита).

Дигитализирането по този начин и директното съхраняването на отчетите се нарича линейна импулсно кодова модулация. В такъв формат за записани CD-тата и wav файловете.

Компютрите излъчват звуков сигнал, като съхранети отчети са подават към устройство генериращо електрически сигнал, който се подава към тонколоните.


Сeмплиране на сигнал


bit rate

Количество битове необходимо за запис на 1 секунда звук.

CD-to има фиксирана семплираща честота =>

44100samples/sec * 16bits/samples = 705600bits/sec

Сравнително големият размер на аудио файловете в CD формат, ни дава основание да търсим по ефективни методи за съхранение на звук.

Компресия без загуби

flac WavePack

Компресия със загуби

  • mp3
  • GSM телефония


Принцип на работа

Заместване на (голям) набор от данни с друг (по-малък) набор от моделиращи коефициенти, които заместват данните чрез минимизиране на разликите между модела и данните.


...
...

Задача. 1. Да се намери y = f(x) по зададени точки

x = [1 2 4 5];
y = [1 2 2 3];


  • Намиране на най-близката функция от първи ред - y = a + b*x
A = [ [1 1 1 1]' , x'];
z = A\y';
%коефиценти на функцията
a = z(1)
b = z(2) 

%генериране на стойности за x
x1 =linspace(0,6,20);
%стойност на функцията за тези стойности
y1 = a + b*x1;
%изчертаване на дадените точки и най-близката права
plot(x1,y1,'-b',x,y,'*r')
grid on


  • Намиране на апроксимираща функция от трети ред у = а + b*x + c*x^2 +d*x^3
A = [ [1 1 1 1]' , x', x'.^2, x'.^3];
z = A\y';
%коефиценти на функцията
a = z(1)
b = z(2)
c = z(3)
d = z(4)
x1 =linspace(0,6,20);
%стойност на функцията за тези стойности
y1 = a + b*x1+c*x1.^2+d*x1.^3;
%изчертаване на дадените точки и най-близката права
plot(x1,y1,'-b',x,y,'*r')
grid on


Заместващи функции

Оригиналният сигнал се замества от линиейна комбинация на косиносови функции.

Задача 2. Да разгледаме функцията f(t) = cos(t) + 5 cos(2t) + cos(3t) + 2 cos(4t) в интервала 0 < t < 2pi. В този интервал може да заместим функцията с равномерно взети отчети s за стойността на функцията.

...
%Разделяме периода 2pi на броя отчети които се ползват
t = linspace (0,2*pi,50)'; % t = 0, pi/50, 2pi/50, 3pi/50 ... 50pi/50

%За всяка стойност на s = f(t)
s = cos(t) + 5*cos(2*t) + cos(3*t) + 2*cos(4*t); %(1)
  
%Обратно генериране на коефициентите
%Създаваме линейна система уравнения
A = [cos(0*t), cos(t), cos(2*t), cos(3*t), cos(4*t)];
z = A\s

%За решения се получават същите коефициенти като в %(1)
plot(t,s);

Заместваща функция: Коефициентите пред cos(0*t),cos(t) и cos(3*t) са малки, затова ги игнорираме.

A = [cos(2*t), cos(4*t)];
z = A\s
s = z(1)*cos(2*t)  + z(2)*cos(4*t);
hold;
plot(t,s,'r');

Обработка на звук в MATLAB

1. Звуков файл в ЛИКМ формат.

 http://ilianko.com/audio/audio.wav

2. Прочитане на звуковия файл.

[s, Fs] = wavread('audio.wav');
% s - стойност на отчет
% Fs - стойност семплиращата честота

3. Възпроизвеждане на звук

sound(s, Fs);
plot(s, (0:length(s))/Fs)

Семплираният звук изглежда по-сложен от разглежданите по-горе примери. Въпреки това данните биха могли да се да се апроксимират по подобен начин. За базовa функция ще се използва косиносова функция. Моделиращата функция би изглеждала така:

y = c0 + c1*cos(ω*t) + c2cos(2*ω*t) + · · · + cn−1*cos((n-1)*ω*t)

Като максималната честота (n-1)*ω според теоерамата на Котелников-Шeнон-Найкуист, трябва да е два пъти по-голяма от честотата на семплирания сигнал.

Изчисляване на модел с ДПФ(DCT)

Нека s съдържа една секунда семплиран звук, с честота на семплиране Fs . В този случай s има Fs наброй стойности.

То моделът би трябвало да се намери по този начин:

t = linspace(0,1,Fs); % време на отчета
А = [cos(0*2*pi*t)), cos(1*2*pi*t), cos(2*2*pi*t), cos(3*2*pi*t), ..., cos((Fs/2-1)*2*pi*t)];
x = A\b;

...(44100 x 22050)

x = dct(s);

Fs = 44100;
t = linspace (0,1,Fs)';
s = cos(2*pi*t) + 5*cos(2*2*pi*t) + cos(3*2*pi*t) + 2*cos(4*2*pi*t);
x = dct(s);
w = sqrt(2/Fs);
f = linspace(0, Fs/2, Fs)';
plot (f(1:10),w*x(1:10),'x');

реконструкция на оригиналния сигнал

y = idct(x);

Цифров филтър

[s,Fs] = wavread ('abc.wav');
s = s/max(s);
N = length(s);
x = dct(s); % изчисляване на апроксимиращия модел
w = sqrt(2/N);
f = linspace(0,Fs/2,N)';
plot (f,w*x); % визуализира коефициентите на съществуващите честоти
hold on
m = (f<3000); % генериране на маска за честотите на 3000 Hz
plot (f,w*m.*x,'r');

y = idct(m.*x); % обратна трансформация, без филтрираните честоти
sound(y,Fs);

Задача: Експеримментирайте с няколко стойности на режащата честота

Задача: Създайте маска, която да режи честотите между 200 и 5000 Херца

задача: С колко се намаля бит райта в горния пример


Идея на mp3

В мп3 вместо отрязване на честотната лента, честотите с по малка значимост се предствавят с по-малка прецизност. Честоти с по малка значимост са тези, чиито коефициенти са с относително малка стойност.

Примерно коефициенти с висока прецизност ще се съхраняват с 16 бита, а тези с малка с 8


Функция за квантуване

function y = quantize (x, bits)

m = max(abs(x));
y = x/m;
y = floor((2^bits - 1)*y/2);
y = 2*y/(2^bits -1);
y = m*y;

Примерна компресия

 % Зареждане на аудио файл
[s, Fs] = wavread ('audio.wav');
% Извличане на 10 сек.
 s = s(44100*20:44100*30,1);
N = length(s);
% Преминаване в честототмна област
x = dct(s);
w = sqrt(2/N);
f = linspace(0,Fs/2,N)';
% Коефициенти
plot (f,w*x)
pause;
% прагова стойност
cutoff = 1
mask = (abs(w*x)<cutoff);
low=mask.*x;
high=(1-mask).*x;
% Визуализация прагови стойности
plot(f,w*high,'r',f,w*low,'b')
% Кванизация
lowbits=8
low = quantize(low, lowbits);
% Реконструиране на сигнала!
y=idct(low+high);
sound (y,Fs);

Литература / References

Lothar Reichel, Digital Audio Compression, http://www.math.kent.edu/~reichel/courses/intr.num.comp.2/lecture14/lecture14.pdf