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

From Ilianko
 
(53 intermediate revisions by the same user not shown)
Line 2: Line 2:
  
 
Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 44000 пъти в sec) се вземат отчети (samples) за стойността на електричския сигнал. Всяко измерване се съхранява като число с фиксирана точност (примерно 8, 16 бита).  
 
Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 44000 пъти в sec) се вземат отчети (samples) за стойността на електричския сигнал. Всяко измерване се съхранява като число с фиксирана точност (примерно 8, 16 бита).  
 +
 +
 +
==импулсно кодова модулация==
 +
 +
*линейно кодиране
 +
*логаритмично кодиране
 +
 +
 +
 +
  
 
Дигитализирането по този начин и директното съхраняването на отчетите се нарича линейна импулсно кодова модулация. В такъв формат за записани CD-тата и wav файловете.  
 
Дигитализирането по този начин и директното съхраняването на отчетите се нарича линейна импулсно кодова модулация. В такъв формат за записани CD-тата и wav файловете.  
Line 8: Line 18:
  
  
[[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 23: Line 132:
  
 
Сравнително големият размер на аудио файловете в CD формат, ни дава основание да търсим по ефективни методи за съхранение на звук.
 
Сравнително големият размер на аудио файловете в CD формат, ни дава основание да търсим по ефективни методи за съхранение на звук.
 +
 +
== Психоакустичен модел ==
 +
 +
*Честотен обхват 20 Hz - 20 KHz
 +
 +
[[Image:f_A.png]]
 +
 +
*Динамичен диапазон
 +
 +
[[Image:shadow.png]]
  
 
== Компресия без загуби ==
 
== Компресия без загуби ==
 +
 +
flac
 +
WavePack
  
 
== Компресия със загуби ==
 
== Компресия със загуби ==
 +
 +
 +
*mp3
 +
*GSM телефония
 +
  
 
===Принцип на работа===
 
===Принцип на работа===
Line 108: Line 235:
  
 
2. Прочитане на звуковия файл.
 
2. Прочитане на звуковия файл.
  [s, Fs] = waveread('audio.wav');
+
  [s, Fs] = wavread('audio.wav');
 
  % s - стойност на отчет
 
  % s - стойност на отчет
 
  % Fs - стойност семплиращата честота
 
  % Fs - стойност семплиращата честота
Line 116: Line 243:
 
  plot(s, (0:length(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)
 
  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)
  
Line 134: Line 261:
 
   
 
   
 
  Fs = 44100;
 
  Fs = 44100;
  t = linspace (0,1,Fs);
+
  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);
 
  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);
 
  x = dct(s);
 
  w = sqrt(2/Fs);
 
  w = sqrt(2/Fs);
  f = linspace(0, Fs/2, Fs);
+
  f = linspace(0, Fs/2, Fs)';
  plot (f(1:10),w*x(1:10),’x’);
+
  plot (f(1:10),w*x(1:10),'x');
  
 
реконструкция на оригиналния сигнал
 
реконструкция на оригиналния сигнал
Line 145: Line 272:
 
  y = idct(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