Difference between revisions of "Компресия на звук"

From Ilianko
 
(73 intermediate revisions by the same user not shown)
Line 2: Line 2:
  
 
Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 44000 пъти в sec) се вземат отчети (samples) за стойността на електричския сигнал. Всяко измерване се съхранява като число с фиксирана точност (примерно 8, 16 бита).  
 
Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 44000 пъти в sec) се вземат отчети (samples) за стойността на електричския сигнал. Всяко измерване се съхранява като число с фиксирана точност (примерно 8, 16 бита).  
 +
 +
 +
==импулсно кодова модулация==
 +
 +
*линейно кодиране
 +
*логаритмично кодиране
 +
 +
 +
 +
 +
 +
Дигитализирането по този начин и директното съхраняването на отчетите се нарича линейна импулсно кодова модулация. В такъв формат за записани CD-тата и wav файловете.
  
 
Компютрите излъчват звуков сигнал, като съхранети отчети са подават към устройство генериращо електрически сигнал, който се подава към тонколоните.
 
Компютрите излъчват звуков сигнал, като съхранети отчети са подават към устройство генериращо електрически сигнал, който се подава към тонколоните.
  
[[Image:sampling.png|thumb|none|650px|Самплиране на сигнал]] 
 
 
  
 +
[[Image:sampling.png|thumb|none|650px|аналогов сигнал/ дискретизиран сигнал/ квантуван сигнал - самплиране/ дигитализиране на сигнал]] 
 +
 +
 +
==Цифрова обработка на сигналите==
 +
 +
[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));
 +
 +
[[Image:sound.png]]
 +
 +
stem(ll(250:300,1));
 +
 +
 +
[[Image:sample.png]]
 +
 +
 +
Филтър:
 +
 +
[[Image:AC_.gif|none|frame]]
 +
[[Image:DSP_.png|none|thumb|300px]]
 +
 +
 +
* signal = [0 0 0 0 0 2 2 2 2 0 0 0 0 0 0 ];
 +
[[Image:signal.png]]
 +
* filter = [-1 1]
 +
[[Image:filter.png]]
 +
* conv(signal, filter)
 +
[[Image:Convolution.png]]
 +
 +
===Филтриране IIR ===
 +
 +
*s1 = sinetone(500,44100,1,1);
 +
[[image:sine500.png]]
 +
*s2 = sinetone(1000,44100,1,1);
 +
[[image:sine1000.png]]
 +
*s = s1+s2
 +
[[image:sine500_1000.png]]
 +
 +
*Филтър
 +
[[Image:DSP.png|300px|thumb]]
 +
 +
*НЧФ (Ниско Честотен Филтър)
 +
*Fmax = 500;
 +
*Fs = 44100;
 +
*size = 5;
 +
*[b a] = butter(size, 2* Fmax / Fs);
 +
[[Image:bcoef.png]]
 +
[[Image:acoef.png]]
 +
*Входен сигнал
 +
[[Image:soef.png]]
 +
*Филтриран сигнал
 +
*sf = filter(b,a,s);
 +
[[Image:sfilt.png]]
 +
 +
===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);
 +
 +
== Цифрово композиране на музика ==
 +
 +
<math>f = 2^{n/12} \times 440 \,\text{Hz}\,</math>
 +
 +
<code><pre>
 +
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);
  
*mp3 player
+
[H,W] = freqz(b, a, F, Fs);
*GSM телефония
+
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');
 +
</pre></code>
  
 
==bit rate==
 
==bit rate==
Line 19: Line 132:
  
 
Сравнително големият размер на аудио файловете в CD формат, ни дава основание да търсим по ефективни методи за съхранение на звук.
 
Сравнително големият размер на аудио файловете в CD формат, ни дава основание да търсим по ефективни методи за съхранение на звук.
 +
 +
== Психоакустичен модел ==
 +
 +
*Честотен обхват 20 Hz - 20 KHz
 +
 +
[[Image:f_A.png]]
 +
 +
*Динамичен диапазон
 +
 +
[[Image:shadow.png]]
  
 
== Компресия без загуби ==
 
== Компресия без загуби ==
 +
 +
flac
 +
WavePack
  
 
== Компресия със загуби ==
 
== Компресия със загуби ==
Оригиналният сигнал се замества от линиейна комбинация на косиносови функции.
+
 
 +
 
 +
*mp3
 +
*GSM телефония
 +
 
 +
 
 +
===Принцип на работа===
 +
Заместване на (голям) набор от данни с друг (по-малък) набор от моделиращи коефициенти, които заместват данните чрез минимизиране на разликите между модела и данните.
 +
 
  
 
[[Image:line.png|right|thumb|333px|...]]
 
[[Image:line.png|right|thumb|333px|...]]
 
[[Image:cubic.png|right|thumb|333px|...]]
 
[[Image:cubic.png|right|thumb|333px|...]]
 +
 
Задача. 1. Да се намери y = f(x) по зададени точки
 
Задача. 1. Да се намери y = f(x) по зададени точки
 
  x = [1 2 4 5];
 
  x = [1 2 4 5];
Line 47: Line 182:
 
  plot(x1,y1,'-b',x,y,'*r')
 
  plot(x1,y1,'-b',x,y,'*r')
 
  grid on
 
  grid on
 +
  
 
*Намиране на апроксимираща функция от трети ред у = а + b*x + c*x^2 +d*x^3
 
*Намиране на апроксимираща функция от трети ред у = а + b*x + c*x^2 +d*x^3
Line 63: Line 199:
 
  grid on
 
  grid on
  
Заместване на (голям) набор от данни с друг (по-малък) набор от моделиращи коефициенти, които заместват данните чрез минимизиране на разликите между модела и данните.
+
 
 +
===Заместващи функции===
 +
Оригиналният сигнал се замества от линиейна комбинация на косиносови функции.
 +
 
 +
'''Задача 2.'''
 +
Да разгледаме функцията '' '''f(t) = cos(t) + 5 cos(2t) + cos(3t) + 2 cos(4t)''' '' в интервала 0 < t < 2pi.
 +
В този интервал може да заместим функцията с равномерно взети отчети '' '''s''' ''за стойността на функцията.
 +
 +
[[Image:cosine.png|thumb|330px|right|...]]
 +
%Разделяме периода 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 = c<sub>0</sub> + c<sub>1</sub>*cos(ω*t) + c<sub>2</sub>cos(2*ω*t) + · · · + c<sub>n−1</sub>*cos((n-1)*ω*t)
 +
 
 +
Като максималната честота  (n-1)*ω според теоерамата на [http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem Котелников-Ш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;
 +
 
 +
Примерна компресия
 +
<code><pre>
 +
% Зареждане на аудио файл
 +
[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);
 +
</pre></code>
 +
 
 +
== Литература / References ==
 +
[http://ilianko.com/files/lecture14.pdf Lothar Reichel, Digital Audio Compression, http://www.math.kent.edu/~reichel/courses/intr.num.comp.2/lecture14/lecture14.pdf]
 +
 
 +
[[Category:ММС]]

Latest revision as of 21:30, 6 November 2015

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

Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 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));

Sound.png

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


Sample.png


Филтър:

AC .gif
DSP .png


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

Signal.png

  • filter = [-1 1]

Filter.png

  • conv(signal, filter)

Convolution.png

Филтриране IIR

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

Sine500.png

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

Sine1000.png

  • s = s1+s2

Sine500 1000.png

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

Bcoef.png Acoef.png

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

Soef.png

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

Sfilt.png

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

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

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

F A.png

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

Shadow.png

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

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