Компресия на звук
Звукът най-често е в следствие на движение на тяло във някаква среда (въздух,...). Движението предизвиква промяна на налягането, което се разпространява както вълна във водата. Тъпънчето на ухото преобразува промяната на налягането в сигнал, който мозъка ни възприема като звук.
Компютрите използват микрофон вместо тъпанче за преобразуването на звуковото налягане в електрически сигнал. След това на определен интервал от време (примерно - 44000 пъти в sec) се вземат отчети (samples) за стойността на електричския сигнал. Всяко измерване се съхранява като число с фиксирана точност (примерно 8, 16 бита).
Дигитализирането по този начин и директното съхраняването на отчетите се нарича линейна импулсно кодова модулация. В такъв формат за записани CD-тата и wav файловете.
Компютрите излъчват звуков сигнал, като съхранети отчети са подават към устройство генериращо електрически сигнал, който се подава към тонколоните.
- mp3 player
- GSM телефония
Contents
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] = 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;
Примерна компресия
% Load an audio sample data set [s, Fs] = wavread (abc.wav); N = length(s); % Compute the interpolation model coefficients x = dct(s); w = sqrt(2/N); f = linspace(0,Fs/2,N)’; % Lets look at the weighted coefficients and pick a cut-off value plot (f,w*c) % Pick a cut-off value and split the coefficients into low- and high-precision sets: cutoff = 0.00015 mask = (abs(w*x)<cutoff); low=mask.*x; high=(1-mask).*x; % This plot nicely illustrates the cut-off region: plot(f,w*high,’-R’,f,w*low,’-b’) % Now pick a precision (in bits) for the low precision data set: lowbits=8 % We wont quantize the high-precision set of coefficients (high), only the % low precision part (requires quantize.m): low = quantize(low, lowbits); % Finally, let’s reconstruct our compressed audio sample and listen to it! y=idct(low+high); sound (y,Fs);
