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

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

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

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

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

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

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

bit rate

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

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

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

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

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

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

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

...

...

Задача. 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] = waveread('audio.wav');
% s - стойност на отчет
% Fs - стойност семплиращата честота

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

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

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

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');
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);

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