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

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

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

импулсно кодова модулация

• линейно кодиране
• логаритмично кодиране

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

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

аналогов сигнал/ дискретизиран сигнал/ квантуван сигнал - самплиране/ дигитализиране на сигнал

Цифрова обработка на сигналите

[ll Fs] = wavread('LL.wav');

ll =
  0.07413   0.04541
  0.06058   0.02930
  0.05316   0.00235
  0.04892  -0.02866
  0.04269  -0.05762
  0.03540  -0.08215
  0.02930  -0.10242
  0.02271  -0.11783
  0.01526  -0.12515
  0.00735  -0.12766
 -0.00253  -0.13477
 -0.01025  -0.14554
 -0.01047  -0.15143

plot(ll(:,1));

stem(ll(250:300,1));

Филтър:

• signal = [0 0 0 0 0 2 2 2 2 0 0 0 0 0 0 ];

• filter = [-1 1]

• conv(signal, filter)

Филтриране IIR

• s1 = sinetone(500,44100,1,1);

• s2 = sinetone(1000,44100,1,1);

• s = s1+s2

• Филтър

• НЧФ (Ниско Честотен Филтър)
• Fmax = 500;
• Fs = 44100;
• size = 5;
• [b a] = butter(size, 2* Fmax / Fs);

• Входен сигнал

• Филтриран сигнал
• sf = filter(b,a,s);

efekti

• echo = 1/5*Fs;
• LL_echo = zeros( length(ll) + echo);
• LL_echo = [zeros(echo,1)', ll(:,1)'] + [ll(:,1)', zeros(echo,1)'];
• sound(LL_echo, Fs);

Цифрово композиране на музика

$f=2^{n/12}\times 440\text{Hz}$

A = 110
Fs = 44100
x = zeros(Fs*4, 1);
F = linspace(1/Fs, 1000, 2^12);
fret  = 4;
delay = round(Fs/(A*2^(fret/12)));

b  = firls(42, [0 1/delay 2/delay 1], [0 0 1 1]);
a  = [1 zeros(1, delay) -0.5 -0.5];

zi = rand(max(length(b),length(a))-1,1);
note = filter(b, a, x, zi);

note = note-mean(note);
note = note/max(abs(note));

sound(note,Fs);

[H,W] = freqz(b, a, F, Fs);
hold on
plot(W, 20*log10(abs(H)), 'r');
title('Harmonics of the A string');
legend('Open A string', 'A string on the 4th fret');

bit rate

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

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

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

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

Психоакустичен модел

• Честотен обхват 20 Hz - 20 KHz

• Динамичен диапазон

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

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